微分方程数值解

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

微分方程数值解

4.1

当常微分方程能解析求解时,可利用Matlab符号工具箱中的功能找到精确解. 见下例

求解方程,,,. 键入: yyy,,,20

syms x y %定义符号变量diff_equ= ‘D2y+2*Dy-y=0’; %D2y表示,,,Dy= y,yy=dsolve (diff_equ, ‘x’) %定义x为自变量 y=cl*exp ((2^(1/2)-

1)*x+c2*exp (-(2^(1/2)+1)*x)

%表达式中含c1与c2,表示通解.

%初始条件为y (0)=0,,y(0)=1时,按如下方式调用 y=dsolve (diff_equ,

‘y (0)=0’, ‘Dy (0)=1’, ‘x’) y=1/4*2^(1/2)*exp ((2^(1/2)-1)*x)—

1/4*2^(1/2)*exp (-(2^(1/2)+1)*x)

%画出函数y=y (x)的图形

ezplot (y,[-2,2])

图形具体形式请上机试之.

在方程无法获得解析解的情况下,可方便地获得数值解. 下面的例子说明用Matlab求

数值解的方法及应注意的问题.

[例1] 求解范德堡(vander pol)方程

2dxdx2,,,,,(1)0xx 2dtdt

求解高阶方程,必须等价地变换为一阶微分方程组,对本例,通过定义两个新的变量,

实现这一变换

yxydxdt1,2/,, 则令 dydty1/2,

2dydtyyy2/(11)*21,,,,

编写求解程序分为两部分,第一部分为待求解的方程,存盘的文件名为,待求

解方程

的函数名.m,,第二部分为求解主程序,本例中取名为main1.m.

首先编写待求解方程的文件. 文件存盘名为“vdpol.m”. M,

function yprime=vdpol(,)ty

yprime (1)=y (2);

; (1(1)^2)*(2)(1),,yyy

mu=2 yprime=[yprime (1);yprime (2)]; yprime (2)=mu*说明函数

yprime=vdpol中. 定义为自变量,的形式取决于求解方程的阶数,本(,)tyyt 例中,,为解向量,为导数向量. yprime, y(2)yyyy,[(1),(2)],(1)(1)(1),y yprime,,,函数返回vander pol方程的导数列向量. 因为所求结果为方程数

值解,(2)(1),y

所以各向量维数只有在主程序求解时定下精度后才能确定.

主程序定名为main1.m,你可用你所喜欢的其它名子,但vdpol.m除外. clear functions

%调试程序时,放置这一语句是必要的. 它清除前边已编译的存在于内存中的

废弃程序[]=ode23 (‘vdpol’,[0,30],[1,0]); ty,

y1=y (:,1); %解曲线.

y2=y (:,2); %解曲线的导数.

polt ( ‘_ _’) tyty,1,,2,

说明龙格_库塔的2阶与4阶改进型求解公式的实现,其指令分别为:

[]=ode23 (‘f’,tsx,0,options) tx,

[]=ode45 (‘f’,tsx,0,options) tx,

其中可由系统依据精度要求自动设定,亦可由使用者依据实际需要自己确定,分别说明之. ts

(1)若令tstttf,[0,1,,],则输出在指定时刻tttf0,1,,给出,当tstktf,0::时,

输出在区间[0,]ttf的等分点上给出,为步长. k

(2)若tsttft,[0,],0为自变量初值,tf为终值,此时,options决定自变

量的维数,t中的时间点不是等间隔的,这是为了保证所需的相对精度,积分算法改变了步长. 用于t

,3,6设定误差限的参数options可缺省,此时系统设定相对误差为,绝对误差为,若1010

自行设定误差限,可用如下语句:

options=odeset (‘reltol’,, ‘abstol’,) rtat

这里的与分别为设定的相对与绝对误差. rtat

须注意的是无论用哪种方法确定ttf0,的取值方式,必须由使用者确定且应与

相匹配. x0t

,y0[1,0],ts,[0,30]y(0)1,y(0)0,为初始条件,本例中,因为,这意味着解曲线,,x0

一般说,当解nnn个未知函数的方程组时,为维向量,共含有个初始条件. x0

两个输出参数是列向量xx与矩阵,它们具有相同的行数,而矩阵的列数等于

方程t

组的个数,本例中y(:,1)y(:,2)的列数为2,其中,为自变量上各点函数值,为上各ytt

点导数值.

最后,提请读者注意的是:ode45也不总是比ode23好,在很多时候,低阶算法更有

效,有关微分方程数值解法的更进一步信息,请参考数值分析方面的书籍. 有些参考书提

供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法. 4.2 -

设有一阶方程与初始条件

,yfxy,(,), (4.1) ,yxy(),00,

其中适当光滑,关于满足Lipschitz条件,即存在使 fxy(,)Ly

fxyfxyLyy(,)(,),,,1212

则(4.1)式的解存在且惟一.

关于yyx,()的解析解一般难以求到或根本无解析解,因而,实际问题中,通常,采

用差分的方法. 在一系列离散点xxx,,,,yyy,,,,上寻求其数值近似解. 12n12n

相关文档
最新文档