微分方程数值解

微分方程数值解 4.1 当常微分方程能解析求解时，可利用Matlab符号工具箱中的功能找到精确解. 见下例 求解方程,,,. 键入： yyy，，,20 syms x y %定义符号变量 diff_equ= ‘D2y+2*Dy-y=0’; %D2y 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示,,,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 相邻两个节点间的间距xxnhn,，,,1,2,hxx,,称为步长，一般地取等步长，则 hn0nnn，1 1、欧拉方法 在区间[,]xx上用差商 nn，1 yxyx()(),nn，1 h代替（4.1）式中,[,]xxxxxy，对fxy(,)中在上取值还是，而形成向前欧拉公式nn，1nn，1 与向后欧拉公式. （1）向前欧拉公式 xfxy(,)取左端点，得如下公式 n yxyxhfxyx()()(,()),， （4.2） nnnn，1 从yxy(),x点出发，由初值代入（4.2）求得 000 yyhfxy,，(,) （4.3） 1000 反复利用（4.2），有 yyhfxyn,，,(,) 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 yyx,()n，1这一过程继续下去，形成折线表示在xy处的精确值，为解的近似值，不难得到 n，1n，1 2h32,,()()()()yxyyxOhOh,,，, nnn，，112 P，1这一误差称为局部截断误差. 若一种算法局部截断误差为Oh()，则称该算法具有阶P 精度，所以向前欧拉公式具有1阶精度. （2）向后欧拉公式 若[,]xxxx中取中的，则有如下公式： fxy(,)nn，1n，1 yyhfxyn,，,(,) 0,1,2, （4.5） nnnn，，，111 称式（4.5）为向后欧拉公式，因为此式中y未知，故称其为隐式公式，无法用其直n，1接计算y，一般用向前欧拉公式产生初值. n，1 (0)yyhfxyn,，,(,) 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,()上，按式 nnn yyx,() （4.4）及式（4.5）分别得 ,P点为与， ABn，1 B且P AB,yyx,()点一定在积分曲线上相应 n点的上、下两边,所以将式（4.4）与（4.5） , 平均之，一定能得到更好的结果. xxx （3）梯形公式 nn，1 将向前与向后欧拉公式加以平均得到所 图4.2 谓梯形公式 hyyfxyfxyn,，，,[(,)(,)] 0,1,2, （4.6） nnnnnn，，，1112 3其局部截断误差为Oh()，具有2阶精度. （4）改进的欧拉公式 为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步 yyhfxy,，(,)nnnn，1 h （4.7） yyfxyfxyn,，，,[(,)(,)], 0,1,2,nnnnn，，11n，12或写为 h,yykk,，，()nn，112,2,1nn （4.8） n,0,1,2,kfxy,(,), ,211nn，kfxyhk,，(,), , 最后指出，上述欧拉方法可推广至微分方程组，如 ,yfxyz,(,,), ,,zgxyz,(,,), ,yxy(),00, ,zxz(),,00向前欧拉公式为 yyhfxyz,，(,,),nnnnn，1 n,0,1,2, ,zzhgxyz,，(,,),nnnnn，1 2、龙格_库塔方法 由微分中值定理 ,[()()]/(),01yxyxhyxh,,，,,,, nnn，1 又因为,,yxhfxhyxh()(,())，,，，,,,yfxy,(,)，所以 nnn从而有 yxyxhfxhyxh()()(),(),，，，,, （4.9） nnnn，1令[,]xx，称其为区间上的平均斜率，由（4.9）可知，给kfxhyxh,，，(,()),,nn，1nn 出一种平均斜率，可相应导出一种算法. 向前欧拉公式中，精度低. 改进欧kfxy,(,)nn 1拉公式中取[,]xxkfxyfxy,，((,)(,))，精度提高，下面，我们在区间内nn，1nnn，1n，12 多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体 给出，只开列具体结果. （1）2阶龙格_库塔公式 yyhkk,，，(),,,nn，11122,kfxy,(,) （4.10） 1nn, ,21nnkfxahyhka,，，,,(,),0,1,,, 1,其中,,,，,,,1,,1a，由于4个未知数只有3个方程，所以解不惟一，若令1222a 1，即得改进的欧拉公式，具有2阶精度. ,,,，,,,,1a122 （3）4阶龙格_库塔公式 只给出精典格式中最常用的一种. 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中一根长的细线，一端固定，另一端悬挂质量为 l m的小球，在重力作用下，小球处于竖直的平衡位置. 现使 小球偏离平衡位置一个小的角度，然后使其自由运动，在不 , 考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动. 为平衡位置，在小球摆动过程中，当与平衡位置夹 ,,0 角为,mgsin,时，小球所受重力在其动运轨迹的分量为 , , l（负号表示力的方向使减少），由牛顿第二定律可得微分方 , 程 ,,mltmg,,()sin,, （4.12） 设小球初始偏离角度为,，且初速为0，式（4.12）的初 0 ,始条件为 ,,,,(0),(0)0,, （4.13） 0 mg 当,不大时，，式（4.12）化为线性常系数微 sin,,,0 分方程 图4.3 g,,,,，,0 （4.14） l 解得 g,,()costt, （4.15） 0l l简谐运动的周期为. T,2,g 现在的问题是：当,较大时，仍用近似，误差太大，式（4.12）又无解析解，,sin,0 试用数值方法在,,30,10,,两种情况下求解，画出的图形，与近似解（4.15）,()t00 比较，这里设. l,25cm [例3] 捕食与被捕食 当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻，小鱼的数量为，鲨鱼的数量为，xt()yt()t 当甲独立生存时它的（相对）增长率与种群数量成正比，即有,r，为增长率，xtrxt()(),而乙的存在使甲的增长率r减少，设减少率与乙的数量成正比，而得微分方程 ,xtxtraytrxaxy()()(()),,,, （4.16） 比例系数a反映捕食者掠取食饵的能力. 乙离开甲无法生存，设乙独自存在时死亡率为,，ytdyt()(),,，甲为乙提供食物，d 使乙的死亡率降低，而促其数量增长，这一作用与甲的数量成正比，于是yt()满足 d ,ytydbxtdybxy()(()),,，,,， （4.17） 比例系数反映甲对乙的供养能力，设若甲，乙的初始数量分别为 b xxyy(0);(0),, （4.18） 00 则微分方程（4.16），（4.17）及初始条件（4.18）确定了甲，乙数量xt()、yt()随时间变化而演变的过程，但该方程无解析解，试用数值解讨论以下问题： （1）设rdabxy,,,,,,1,0.5,0.1,0.02,25,2，求方程（4.16），（4.17）在00 条件（4.18）下的数值解，画出xtyt(),()的图形及相图(,)xy，观察解xtyt(),()的周期变化，近似确定解的周期和的最大、小值，近似计算在一个周期内的平均值. xy,xy, （2）从式（4.16）和（4.17）消去得到 dt dxxray(),, （4.19） dyydbx(),， 解方程（4.19），得到的解即为相轨线，说明这是封闭曲线，即解确为周期函数. （3）将方程（4.17）改写为 ,1yxtd()[],， （4.20） by 在一个周期内积分，得到xt()yt()一周期内的平均值，类似可得一周期内的平均值，将近似计算的结果与理论值比较. 2、实验 （1）方程（单摆问题） ,,mlmg,,,,sin, ,,,,,(0),(0)0,,0, 无解析解，为求其数值解，先将其化为等价的一阶方程组. 令,xx,,,,,，方程化为 12 ,,xx,12,g,,21,,xxsin ,l,102,,xx(0),(0)0,,, 其中，gl,,9.8,25,,为（弧度）与（弧度）两种情况，具100.1745,300.5236,0 体编程如下： 先以danbai.m为文件名存放待解方程. 键入： function xdot =danbai (t,x) %x=[x (1),x (2)], g=9.8;1=25; %x (1)为解向量, x (2)是导数 xdot (1)=x (2); %xdot (1)=,,(1)= x,xdot (2)=-g/1*sin (x (1)); %xdot (2)=,, , xdot=[xdot (1);xdot (2)]; %必须将导数向量以列向量形式给出. 再以主程序(求数值解)调用待求方程，主程序用main2.m为文件名存盘，其代码如下： clear functions [t,x]=ode23 (‘danbai’,[0,10],[0.1745,0]) %只计算,,100.1745,（弧度）的情形.0 1g,,()cos,twtw,,%对近似解，周期 T2,,01gw=sqrt (9.8/25); y=0.1745*cos (w*t); [t,x (:,1),y] %显示数据，无分号. hold on %欲在同一幅图中画近似解. plot (t,x (:,1), ‘b’) %画数值解，绿色 plot (t,y, ‘g*’) %用*号，红色画近似解. Hold off （2）食饵_捕食者 方程 ,xtrxaxy(),, ,ytdybxy(),,， 可化为如下形式 ，，xrayx,0，，，，,,, ,,,,,,0,，dbxy，，，，y，， 初始条件xxyy(0),(0),,表示为 00 xx(0)，，，，0 ,,,,,yy(0)，，，，0以shier.m存盘如下代码 function xdot=shier (t,x) r=1;d=0.5;a=0.1;b=0.02; xdot=diag ([r-a*x (2),-d+b*x (1)])*x; ，，xxx(1)，，，，%x=,,,xdot= ,,,,,,,xy(2)，，，，y，， 以main3.m存盘如下代码. clear functions ts=0:0.1:15; x0=[25,2]; [t,x]=ode45 (‘shier’,ts,x0); [t,x] %显示数据t,x,y plot (t,x) gird %加网格线 gtext (‘x(t)’),gtext (‘y(t)’), %用点鼠标方式 pause,figure (2) %将x (t),x(t)放至指定点 12plot (x (:,1),x (:,2)) %以x (1)与x (2)为坐标点画相图 xlabel (‘x’),ylabel (‘y’) 可以猜测，xt()xt()(,)xx与是周期函数，相图是封闭曲线，从数值解可近似定出2121 周期为10.7,x2.0,x的最大和最小值分别为与的最大和最小值分别为和，99.328.42.012为求xx与在一个周期的平均值，只需键入： 12 y1=x (1:108,1); %x周期为10.7. 1 x1p=trapz (y1)*0.1/10.7, %trapz (y1)返回按 y2=x (1:108,2); %梯形法对y1的积 107yiyi1()1(1)，，x2p=trapz (y2)*0.1/10.7, %分值,x ,i12i,可得xx,,25,10 12 对方程 dxxray(),,，,dbxray,dxdy,化为 dyydbx(),，xy 积分得解 dbxray,,()()xeyec, （4.21） 即为原方程组的相轨线，其中c由初始条件确定. 为说明上述相轨线是封闭的，令 dbxray,,fxxegyye();(),, 设其最大值分别为fg,xy,，若满足 mm00 fxffyg(),(),, 00mm dr则有,,xy,fgc,,（令，可解出），又当时，相轨线（4.21）fg,,0,0xy,,,00mm00ba 有意义. 当fgc,0,,cfgPxy(,)时，相轨线退化为一个点，对时，相轨线如mmmm00 图4.4（c），而图（a），（b）分别为fx()与gy()的图像，下面讨论如何由（a），（b）画（c）. fxf(),nm y gyg(),gv()nm fx() y2 gm fm Pxy(,)qq 00 1 yP q02 q3 y1x y x x xxxy xyxy012001221 （a） （b） （c） 图4.4 设cpgpf,,,(0)yy,，xx,fxp(),，若令，则有，由（a）知，，使mm012 qxy(,)xxx,,qxy(,)fxfxp()(),,，且于是相轨线应过与，对11010222012 fxgxpggyg()()(),,,xxx,(,)fxp(),gyq(),，，由，令，又由（b）知，mm12 存在yyy,,qxy(,)yy,gygyq()(),,qxy(,)使，且，于是相轨线又过与10231121242 yyyyy,(),,xqq,两点，所以对间每个，轨线总要通过纵坐标为的两点，于是1210212 相轨线是一条封闭曲线. 相轨线封闭等价于xtyt(),()是周期函数，记周期为，为求其在一个周期内的平均T 值(,)xy，由 1yxtd()(),， by 两边在一个周期内积分有((0)())yyT,： T11ln()ln(0)yTydTd,，，xxtdt,,，,() ,,,0TTbbb，， 同理 ry, a 从而 xxyy,,,00 即xy的平均值恰为相轨线中心点的坐标，提请读者注意的是：这里的与与xtyt(),()p00 初始条件中的xy,不是一件事. 00 注意到在生态学上的意义，上述结果表明，捕食者的数量与食饵的增长率成rdab,,, 正比，与它掠取食饵的能力成反比，食饵的数量与捕食者死亡率成正比，与它供养捕食者 的能力成反比. 3、练习内容 （1）编写改进欧拉公式求微分方程数值解的程序，并用其与ode23求下列微分方程数值解，对二者作出比较. a）22,yxyy,,,,(0)0或y(0)1,. 222 b),,,xyxyxny，，,,()0 2,,,21 yxsin,n,,（Bessel方程，这里令，其精确解为）. yy()2,(),,,x222, c),,,yyxyy，,,,cos0,(0)1,(0)0. （2）倒圆锥形容器，上底面直径为1.2m，容器的高亦为1.2m，在锥尖的地方开有一直径为3cm的小孔，容器装满水后，下方小孔开启，由水利学知识可知当水面高度为时，h 水从小孔中流出的速度为为重力加速度，若孔口收缩系数为0.6（即若一个面vghg,2, 积单位的小孔向外出水时，水柱截面积为0.6），问水从小孔中流完需多少时间？2分钟时，水面高度是多少？ （3）一只小船渡过宽为的河流，目标是起点正对着的另一岸上点，已知河水ABd 流速vv与船在静水中的速度之比为. k12 （a）建立小船航线的方程，求其解析解. （b）设dvv,,,100m,1m/s,2m/s，用数值解法求渡河所需时间，任意时刻小12 船的位置及航行曲线，作图并与解析解比较. （c）若流速v为0,0.5,2 (m/s)，结果将如何？ 1 （4）研究种群竞争模型. 当甲、乙两个种群各自生存时，数量演变服从下面规律 xyxtrxytry()(1),()(1),,,, 12nn12 其中，rr,nn,xtyt(),()分别为时刻甲，乙两个种群的数量，为其固有增长率，为它t1212 们的最大容量，而当这两个种群在同一环境中生存时，由于乙消耗有限资源而对甲的增长 产生影响，将甲的方程修改为 xyxrxs,,,(1) （4.22） 11nn12 这里s的含意是：对于供养甲的资源而言，单位数量乙（相对n）的消耗率为单位数12 量甲（相对n）消耗的s倍，类似地，甲的存在亦影响乙的增长，乙的方程应改为 11 xyytrys()(1),,, （4.23） 22nn12给定种群的初始值为 xxyy(0),(0),, （4.24） 00 及参数rrssnn,,,,,后，方程（4.22）与（4.23）确定了两种群的变化规律，因其解析121212 解不存在，试用数值解法研究以下问题： （a）设rrnnssxy,,,,,,,,1,100,0.5,2,10，计算，画出xtyt(),()12121200 它们的图形及相图(,)xy，说明时间充分大以后xtyt(),()的变化趋势（人们今天看到的t 已经是自然界长期演变的结果）. （b）改变rrxy,,,，但s与s不变（保持ss,,1,1），分析所得结果，若12001212ss,,,,1.5(1),0.7(1)，再分析结果，由此你得到什么结论，请用各参数生态学上的12 含义作出解释. （c）试验当ss,,,,0.8(1),0.7(1)ss,,,,1.5(1),1.7(1)时会有什么结果；当1212时又会出现什么结果，能解释这些结果吗？