微分方程数值解.docx

1、微分方程数值解微分方程数值解4.1 当常微分方程能解析求解时，可利用Matlab符号工具箱中的功能找到精确解. 见下例求解方程,. 键入： yyy，,20syms 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

2、=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, 编写求解

3、程序分为两部分，第一部分为待求解的方程，存盘的文件名为,待求解方程的函数名.m,，第二部分为求解主程序，本例中取名为main1.m. 首先编写待求解方程的文件. 文件存盘名为“vdpol.m”. M,function yprime=vdpol(,)ty yprime (1)=y (2); ; (1(1)2)*(2)(1),yyymu=2 yprime=yprime (1);yprime (2); yprime (2)=mu*说明 函数yprime=vdpol中. 定义为自变量，的形式取决于求解方程的阶数，本(,)tyyt例中，,为解向量，为导数向量. yprime， y(2)yyyy,(1),

4、(2),(1)(1)(1),yyprime,，函数返回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阶改进型求解公式的

5、实现，其指令分别为： =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可缺省，此时系统设定相对

6、误差为，绝对误差为，若1010自行设定误差限，可用如下语句： options=odeset (reltol, abstol,) rtat这里的与分别为设定的相对与绝对误差. rtat须注意的是无论用哪种方法确定ttf0,的取值方式，必须由使用者确定且应与相匹配. x0t,y01,0,ts,0,30y(0)1,y(0)0,为初始条件，本例中，因为，这意味着解曲线，x0一般说，当解nnn个未知函数的方程组时，为维向量，共含有个初始条件. x0两个输出参数是列向量xx与矩阵，它们具有相同的行数，而矩阵的列数等于方程t组的个数，本例中y(:,1)y(:,2)的列数为2，其中，为自变量上各点函数值，为上

7、各ytt点导数值. 最后，提请读者注意的是：ode45也不总是比ode23好，在很多时候，低阶算法更有效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法. 4.2 - 设有一阶方程与初始条件 ,yfxy,(,), （4.1） ,yxy(),00,其中适当光滑，关于满足Lipschitz条件，即存在使 fxy(,)LyfxyfxyLyy(,)(,),1212则（4.1）式的解存在且惟一. 关于yyx,()的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采用差分的方法. 在一系列

8、离散点xxx,yyy,上寻求其数值近似解. 12n12n相邻两个节点间的间距xxnhn,，,1,2,hxx,称为步长，一般地取等步长，则 hn0nnn，11、欧拉方法 在区间,xx上用差商 nn，1yxyx()(),nn，1 h代替（4.1）式中,xxxxxy，对fxy(,)中在上取值还是，而形成向前欧拉公式nn，1nn，1与向后欧拉公式. （1）向前欧拉公式 xfxy(,)取左端点，得如下公式 nyxyxhfxyx()()(,(),， （4.2） nnnn，1从yxy(),x点出发，由初值代入（4.2）求得 000yyhfxy,，(,) （4.3） 1000反复利用（4.2），有 yyhfx

9、yn,，,(,) 0,1,2, （4.4） nnnn，1其几何意义如图4.1所示. y 图中yyx,()为方程（4.1）的精确 P P43P 2解曲线，其上任意点(,)xy处切线斜率为 误差 P 1yyx,() 32Pxy(,)fxy(,). 从初值点出发，用该 P000 0y 0点斜率fxy(,)xx,作一直线段，在 001yyx,() yx() 3处得到Pxy(,)y，由（4.2）式确定， 1111y 3再从Pxy(,)fxy(,)出发，以为斜率 11111作直线段，在xx,Pxy(,)处得到， 2222xxxxx x O03124 PPP， 012作为积分曲线yx()的近似，用 图4.1

10、 yyx,()n，1这一过程继续下去，形成折线表示在xy处的精确值，为解的近似值，不难得到 n，1n，12h32,()()()()yxyyxOhOh,，, nnn，112P，1这一误差称为局部截断误差. 若一种算法局部截断误差为Oh()，则称该算法具有阶P精度，所以向前欧拉公式具有1阶精度. （2）向后欧拉公式 若,xxxx中取中的，则有如下公式： fxy(,)nn，1n，1yyhfxyn,，,(,) 0,1,2, （4.5） nnnn，111称式（4.5）为向后欧拉公式，因为此式中y未知，故称其为隐式公式，无法用其直n，1接计算y，一般用向前欧拉公式产生初值. n，1(0)yyhfxyn,，

11、,(,) 0,1,2, 11nnnn，再按下式迭代 (1)()kk，yyhfxykn,，,(,),0,1,0,1, nnnn，111其误差估计如下 2h32,()()()()yxyyxohoh,，, nnn，112精度亦为1阶，将向前欧拉公式（4.4）与向后欧拉公式（4.5）及它们的误差的几何说明作一对比，是十分有益的，见图4.2. y 为讨论局部截断误差，在图4.2中设点 APxy(,)落在积分曲线yyx,()上，按式 nnnyyx,() （4.4）及式（4.5）分别得 ,P点为与， ABn，1 B且P AB,yyx,()点一定在积分曲线上相应 n点的上、下两边,所以将式（4.4）与（4.5

12、） ,平均之，一定能得到更好的结果. xxx （3）梯形公式 nn，1 将向前与向后欧拉公式加以平均得到所 图4.2 谓梯形公式 hyyfxyfxyn,，,(,)(,) 0,1,2, （4.6） nnnnnn，11123其局部截断误差为Oh()，具有2阶精度. （4）改进的欧拉公式 为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步 yyhfxy,，(,)nnnn，1h （4.7） yyfxyfxyn,，,(,)(,), 0,1,2,nnnnn，11n，12或写为 h,yykk,，()nn，112,2,1nn （4.8） n,0,1,2,kfxy,(,),211nn，kfxyhk,，(

13、,),最后指出，上述欧拉方法可推广至微分方程组，如 ,yfxyz,(,),zgxyz,(,), ,yxy(),00,zxz(),00向前欧拉公式为 yyhfxyz,，(,),nnnnn，1 n,0,1,2, ,zzhgxyz,，(,),nnnnn，12、龙格_库塔方法 由微分中值定理 ,()()/(),01yxyxhyxh,，, nnn，1又因为,yxhfxhyxh()(,()，,，,yfxy,(,)，所以 nnn从而有 yxyxhfxhyxh()()(),(),，, （4.9） nnnn，1令,xx，称其为区间上的平均斜率，由（4.9）可知，给kfxhyxh,，(,(),nn，1nn出一种平

14、均斜率，可相应导出一种算法. 向前欧拉公式中，精度低. 改进欧kfxy,(,)nn1拉公式中取,xxkfxyfxy,，(,)(,)，精度提高，下面，我们在区间内nn，1nnn，1n，12多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体给出，只开列具体结果. （1）2阶龙格_库塔公式 yyhkk,，(),nn，11122,kfxy,(,) （4.10） 1nn,21nnkfxahyhka,，,(,),0,1,1,其中,，,1,1a，由于4个未知数只有3个方程，所以解不惟一，若令1222a1，即得改进的欧拉公式，具有2阶精度. ,，,1a122（3）4阶龙格_库塔公式

15、 只给出精典格式中最常用的一种. h,yykkkk,，(22)nn，11234,6,kfxy,(,),1nn,hh, （4.11） kfxyk,，(,),21nn22,hh,kfxyk,，(,)32nn,22,kfxhyhk,，(,),43nn,其计算精度为4阶 4.3 1、模型与问题 例2 单摆运动 图4.3中一根长的细线，一端固定，另一端悬挂质量为 lm的小球，在重力作用下，小球处于竖直的平衡位置. 现使 小球偏离平衡位置一个小的角度，然后使其自由运动，在不 ,考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动. 为平衡位置，在小球摆动过程中，当与平衡位置夹 ,0角为,mgsin,时，小球所