下载文档原格式
常微分方程的数值解法在自然科学的许多领域中，都会遇到常微分方程的求解问题。
然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。
在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。
这些方法可以给出解的近似表达式，通常称为近似解析方法。
还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。
利用计算机解微分方程主要使用数值方法。
我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(yx y y x f dx dy在区间[a, b]上的解，其中f (x, y )为x, y 的已知函数，y 0为给定的初始值，将上述问题的精确解记为y (x )。
数值方法的基本思想是：在解的存在区间上取n + 1个节点b x x x x a n =<<<<= 210这里差i i i x x h -=+1，i = 0,1, …, n 称为由x i 到x i +1的步长。
这些h i 可以不相等，但一般取成相等的，这时na b h -=。
在这些节点上采用离散化方法，（通常用数值积分、微分。
泰勒展开等）将上述初值问题化成关于离散变量的相应问题。
把这个相应问题的解y n 作为y (x n )的近似值。
这样求得的y n 就是上述初值问题在节点x n 上的数值解。
一般说来，不同的离散化导致不同的方法。
§1 欧拉法与改进欧拉法 1.欧拉法1．对常微分方程初始问题(9.2))((9.1) ),(00⎪⎩⎪⎨⎧==y x y y x f dx dy用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从（9.2）式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。
第四章 微分方程数值解4.1 算 法当常微分方程能解析求解时，可利用Matlab 符号工具箱中的功能找到精确解. 见下例求解方程20y y y '''++=. 键入：syms x y %定义符号变量 diff_equ= ‘D2y+2*Dy -y=0’; %D2y 表示,y ''Dy=y 'y=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 ）方程222(1)0d x dx x x dt dtμ--+=求解高阶方程，必须等价地变换为一阶微分方程组，对本例，通过定义两个新的变量，实现这一变换1,2/y x y dx dt ==则令 1/2dy dt y =22/(11)*21dy dt y y y μ=--编写求解程序分为两部分，第一部分为待求解的方程，存盘的文件名为<待求解方程的函数名.m >，第二部分为求解主程序，本例中取名为main1.m.首先编写待求解方程的M -文件. 文件存盘名为“vdpol.m ”. function yprime=vdpol (,)t y yprime (1)=y (2);mu=2yprime (2)=mu*(1(1)^2)*(2)(1)y y y --;yprime=[yprime (1);yprime (2)];说明 函数yprime=vdpol (,)t y 中. 定义t 为自变量，y 的形式取决于求解方程的阶数，本例中，[(1),(2)],(1)y y y y =为解向量，(2)y 为导数向量. yprime (1)(1)y '=， yprime (2)(1)y ''=，函数返回vander pol 方程的导数列向量. 因为所求结果为方程数值解，所以各向量维数只有在主程序求解时定下精度后才能确定.主程序定名为main1.m ，你可用你所喜欢的其它名子，但vdpol.m 除外. clear functions%调试程序时，放置这一语句是必要的. 它清除前边已编译的存在于内存中的废弃程序 [,t y ]=ode23 (‘vdpol’,[0,30],[1,0]); y1=y (:,1); %解曲线. y2=y (:,2); %解曲线的导数. polt (,1,,2,t y t y ‘_ _’)说明 龙格_库塔的2阶与4阶改进型求解公式的实现，其指令分别为： [,t x ]=ode23 (‘f ’,,0,ts x options) [,t x ]=ode45 (‘f ’,,0,ts x options)其中ts 可由系统依据精度要求自动设定，亦可由使用者依据实际需要自己确定，分别说明之. （1）若令[0,1,,]ts t t tf =，则输出在指定时刻0,1,,t t tf 给出，当0::ts t k tf =时，输出在区间[0,]t tf 的等分点上给出，k 为步长.（2）若[0,],0ts t tf t =为自变量初值，tf 为终值，此时，options 决定自变量t 的维数，t 中的时间点不是等间隔的，这是为了保证所需的相对精度，积分算法改变了步长. 用于设定误差限的参数options 可缺省，此时系统设定相对误差为310-，绝对误差为610-，若自行设定误差限，可用如下语句： options=odeset (‘reltol’,rt , ‘abstol’,at ) 这里的rt 与at 分别为设定的相对与绝对误差.须注意的是无论用哪种方法确定t 的取值方式，0,t tf 必须由使用者确定且应与0x 相匹配. 0x 为初始条件，本例中0[1,0]y =，因为[0,30]ts =，这意味着解曲线(0)1y =，(0)0y '=，一般说，当解n 个未知函数的方程组时，0x 为n 维向量，共含有n 个初始条件.两个输出参数是列向量t 与矩阵x ，它们具有相同的行数，而矩阵x 的列数等于方程组的个数，本例中y 的列数为2，其中，(:,1)y 为自变量t 上各点函数值，(:,2)y 为t 上各点导数值.最后，提请读者注意的是：ode45也不总是比ode23好，在很多时候，低阶算法更有效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法.4.2 欧拉与龙格-库塔方法设有一阶方程与初始条件0(,)()y f x y y x y '=⎧⎨=⎩ （4.1） 其中(,)f x y 适当光滑，关于y 满足Lipschitz 条件，即存在L 使1212(,)(,)f x y f x y L y y -≤-则（4.1）式的解存在且惟一.关于()y y x =的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采用差分的方法. 在一系列离散点12n x x x <<<<上寻求其数值近似解12,,,,n y y y .相邻两个节点间的间距1n n n h x x +=-称为步长，一般地取等步长h ，则0,1,2,n x x nh n =+=1、欧拉方法在区间1[,]n n x x +上用差商1()()n n y x y x h+-代替（4.1）式中y '，对(,)f x y 中x 在1[,]n n x x +上取值n x 还是1n x +，而形成向前欧拉公式与向后欧拉公式.（1）向前欧拉公式(,)f x y 取左端点n x ，得如下公式1()()(,())n n n n y x y x hf x y x +≈+ （4.2）从0x 点出发，由初值00()y x y =代入（4.2）求得1000(,)y y hf x y =+ （4.3）反复利用（4.2），有1(,) 0,1,2,n n n n y y hf x y n +=+= （4.4）其几何意义如图4.1所示.图中()y y x =为方程（4.1）的精确 解曲线，其上任意点(,)x y 处切线斜率为(,)f x y . 从初值000(,)P x y 点出发，用该点斜率00(,)f x y 作一直线段，在1x x = 处得到111(,)P x y ，1y 由（4.2）式确定， 再从111(,)P x y 出发，以11(,)f x y 为斜率 作直线段，在2x x =处得到222(,)P x y ，y0yO0x 误差4x1x 2x 3x x()y y x =3()y x32()y y x - 0P1P2P3P 4P3y这一过程继续下去，形成折线012P PP ，作为积分曲线()y y x =的近似，用1()n y x + 图4.1 表示在1n x +处的精确值，1n y +为解的近似值，不难得到23211()()()()2n n n h y x y y x O h O h ++''-=+=这一误差称为局部截断误差. 若一种算法局部截断误差为1()P O h +，则称该算法具有P 阶精度，所以向前欧拉公式具有1阶精度.（2）向后欧拉公式若(,)f x y 中x 取1[,]n n x x +中的1n x +，则有如下公式：111(,) 0,1,2,n n n n y y hf x y n +++=+= （4.5）称式（4.5）为向后欧拉公式，因为此式中1n y +未知，故称其为隐式公式，无法用其直接计算1n y +，一般用向前欧拉公式产生初值.(0)11(,) 0,1,2,n n n n y y hf x y n ++=+=再按下式迭代(1)()111(,),0,1,,0,1,k k n n n n y y hf x y k n ++++=+==其误差估计如下23211()()()()2n n n h y x y y x o h o h ++''-=+=精度亦为1阶，将向前欧拉公式（4.4）与向后欧拉公式（4.5）及它们的误差的几何说明作一对比，是十分有益的，见图4.2. 为讨论局部截断误差，在图4.2中设点 (,)n n n P x y 落在积分曲线()y y x =上，按式 （4.4）及式（4.5）分别得1n P +点为A 与B ， 且,A B 点一定在积分曲线()y y x =上相应 点θ的上、下两边,所以将式（4.4）与（4.5） 平均之，一定能得到更好的结果. （3）梯形公式 将向前与向后欧拉公式加以平均得到所 图4.2谓梯形公式111[(,)(,)] 0,1,2,2n n n n n n hy y f x y f x y n +++=++= （4.6）其局部截断误差为3()O h ，具有2阶精度.（4）改进的欧拉公式为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步1(,)n n n n y y hf x y +=+x 1n x +n x y()y y x =n PAθB111[(,)(,)], 0,1,2,2n n n n n n hy y f x y f x y n +++=++= （4.7）或写为1121211()2(,)(,)n n n n n n h y y k k k f x y k f x y hk ++⎧=++⎪⎪=⎨⎪=+⎪⎩0,1,2,n= （4.8）最后指出，上述欧拉方法可推广至微分方程组，如0000(,,)(,,)()()y f x y z z g x y z y x y z x z '=⎧⎪'=⎪⎨=⎪⎪=⎩ 向前欧拉公式为11(,,)(,,)n n n n n n n n n n y y hf x y z z z hg x y z ++=+⎧⎨=+⎩ 0,1,2,n=2、龙格_库塔方法 由微分中值定理1[()()]/(),01n n n y x y x h y x h θθ+'-=+<<又因为(,)y f x y '=，所以()(,())n n n y x h f x h y x h θθθ'+=++从而有1()()(),()n n n n y x y x hf x h y x h θθ+=+++ （4.9）令(,())n n k f x h y x h θθ=++，称其为区间1[,]n n x x +上的平均斜率，由（4.9）可知，给出一种平均斜率，可相应导出一种算法. 向前欧拉公式中(,)n n k f x y =，精度低. 改进欧拉公式中取111((,)(,))2n n n n k f x y f x y ++=+，精度提高，下面，我们在区间1[,]n n x x +内多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体给出，只开列具体结果. （1）2阶龙格_库塔公式11122121()(,)(,),0,1n n n n n n y y h k k k f x y k f x ah y hk a λλββ+=++⎧⎪=⎨⎪=++<<⎩ （4.10） 其中12211,,12a aβλλλ+===，由于4个未知数只有3个方程，所以解不惟一，若令121,12a λλβ+===，即得改进的欧拉公式，具有2阶精度.（3）4阶龙格_库塔公式 只给出精典格式中最常用的一种.112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y k k k k k f x y h h k f x y k h h k f x y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ （4.11） 其计算精度为4阶4.3 模型与实验1、模型与问题 [例2] 单摆运动图4.3中一根长l 的细线，一端固定，另一端悬挂质量为m 的小球，在重力作用下，小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度θ，然后使其自由运动，在不考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动. 0θ=为平衡位置，在小球摆动过程中，当与平衡位置夹 角为θ时，小球所受重力在其动运轨迹的分量为sin mg θ- （负号表示力的方向使θ减少），由牛顿第二定律可得微分方 程()sin ml t mg θθ''=- （4.12） 设小球初始偏离角度为0θ，且初速为0，式（4.12）的初 始条件为0(0),(0)0θθθ'== （4.13）当0θ不大时，sin θθ≈，式（4.12）化为线性常系数微分方程 图4.30g lθθ''+= （4.14）解得θlθmg0()t θθ= （4.15）简谐运动的周期为2T =. 现在的问题是：当0θ较大时，仍用θ近似sin θ，误差太大，式（4.12）又无解析解，试用数值方法在0030,10θθ==两种情况下求解，画出()t θ的图形，与近似解（4.15）比较，这里设25cm l =. [例3] 捕食与被捕食当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻t ，小鱼的数量为()x t ，鲨鱼的数量为()y t ，当甲独立生存时它的（相对）增长率与种群数量成正比，即有()()x t rx t '=，r 为增长率，而乙的存在使甲的增长率r 减少，设减少率与乙的数量成正比，而得微分方程()()(())x t x t r ay t rx axy '=-=- （4.16）比例系数a 反映捕食者掠取食饵的能力.乙离开甲无法生存，设乙独自存在时死亡率为d ，()()y t dy t '=-，甲为乙提供食物，使乙的死亡率d 降低，而促其数量增长，这一作用与甲的数量成正比，于是()y t 满足()(())y t y d bx t dy bxy '=-+=-+ （4.17）比例系数b 反映甲对乙的供养能力，设若甲，乙的初始数量分别为00(0);(0)x x y y == （4.18） 则微分方程（4.16），（4.17）及初始条件（4.18）确定了甲，乙数量()x t 、()y t 随时间变化而演变的过程，但该方程无解析解，试用数值解讨论以下问题：（1）设001,0.5,0.1,0.02,25,2r d a b x y ======，求方程（4.16），（4.17）在条件（4.18）下的数值解，画出(),()x t y t 的图形及相图(,)x y ，观察解(),()x t y t 的周期变化，近似确定解的周期和,x y 的最大、小值，近似计算,x y 在一个周期内的平均值. （2）从式（4.16）和（4.17）消去dt 得到()()dx x r ay dy y d bx -=-+ （4.19） 解方程（4.19），得到的解即为相轨线，说明这是封闭曲线，即解确为周期函数. （3）将方程（4.17）改写为1()[]y x t d b y'=+ （4.20）在一个周期内积分，得到()x t 一周期内的平均值，类似可得()y t 一周期内的平均值，将近似计算的结果与理论值比较. 2、实验（1）方程（单摆问题）0sin (0),(0)0ml mg θθθθθ''=-⎧⎨'==⎩无解析解，为求其数值解，先将其化为等价的一阶方程组. 令12,x x θθ'==，方程化为1221102sin (0),(0)0x x g x x l x x θ'=⎧⎪⎪'=-⎨⎪==⎪⎩ 其中，09.8,25,g l θ==为100.1745≈（弧度）与300.5236≈（弧度）两种情况，具体编程如下：先以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)=x '(1)=θ' 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])%只计算0100.1745θ==（弧度）的情形. %对近似解0()cos ,t wt w θθ==2T = w=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）食饵_捕食者 方程()x t rx axy '=- ()y t dy bxy '=-+可化为如下形式00x r ay x d bx y y ⎡⎤-⎡⎤⎡⎤⎢⎥=⎢⎥⎢⎥⎢⎥-+⎣⎦⎣⎦⎣⎦初始条件00(0),(0)x x y y ==表示为00(0)(0)x x y y ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦以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;%x=(1)(2)x x x y ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦,xdot=x 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 1(t),x 2(t)放至指定点 plot (x (:,1),x (:,2)) %以x (1)与x (2)为坐标点画相图xlabel (‘x’),ylabel (‘y’)可以猜测，1()x t 与2()x t 是周期函数，相图12(,)x x 是封闭曲线，从数值解可近似定出周期为110.7,x 的最大和最小值分别为99.3与22.0,x 的最大和最小值分别为28.4和2.0，为求1x 与2x 在一个周期的平均值，只需键入： y1=x (1:108,1); %x 1周期为10.7. x1p=trapz (y1)*0.1/10.7, %trapz (y1)返回按 y2=x (1:108,2); %梯形法对y1的积 x2p=trapz (y2)*0.1/10.7, %分值10711()1(1)2i i y i y i x =++∆∑ 可得1225,10x x == 对方程()()dx x r ay dy y d bx -=-+化为d bx r aydx dy x y-+-=积分得解()()d bx r ay x e y e c --= （4.21）即为原方程组的相轨线，其中c 由初始条件确定. 为说明上述相轨线是封闭的，令();()d bx r ay f x x e g y y e --==设其最大值分别为,m m f g ，若00,x y 满足00(),()m m f x f f y g ==则有00,d rx y b a==（令0,0f g ''==，可解出00,x y ），又当m m f g c ⋅≥时，相轨线（4.21）有意义. 当m m f g c =时，相轨线退化为一个点00(,)P x y ，对0m m c f g <<时，相轨线如图 4.4（c ），而图（a ），（b ）分别为()f x 与()g y 的图像，下面讨论如何由（a ），（b ）画（c ）.（a ） （b ） （c ）图4.4设(0)m m c pg p f =<<，若令0y y =，则有()f x p =，由（a ）知，12,x x ∃，使12()()f x f x p ==，且102x x x <<于是相轨线应过110(,)q x y 与220(,)q x y ，对12(,)x x x ∈，()f x p >，由()()()m m f x g x pg g y g =⇒<，令()g y q =，又由（b ）知，存在12,y y 使12()()g y g y q ==，且102y y y <<，于是相轨线又过31(,)q x y 与42(,)q x y 两点，所以对12,q q 间每个x ，轨线总要通过纵坐标为12102,()y y y y y <<的两点，于是相轨线是一条封闭曲线.相轨线封闭等价于(),()x t y t 是周期函数，记周期为T ，为求其在一个周期内的平均值(,)x y ，由1()()yx t d b y=+两边在一个周期内积分有((0)())y y T =：011ln ()ln (0)()T y T y dT d x x t dt T T b b b-⎡⎤==+=⎢⎥⎣⎦⎰ 同理()f xmfP1x 0x 2xxyy 1y ()()n m n mf x fg y g ==2yyqm g()g vx1y 0y 2y1xx 2x1q 3q2q00(,)P x yr y a=从而 00,x x y y ==即(),()x t y t 的平均值恰为相轨线中心点p 的坐标，提请读者注意的是：这里的0x 与0y 与初始条件中的00,x y 不是一件事.注意到,,,r d a b 在生态学上的意义，上述结果表明，捕食者的数量与食饵的增长率成正比，与它掠取食饵的能力成反比，食饵的数量与捕食者死亡率成正比，与它供养捕食者的能力成反比.3、练习内容（1）编写改进欧拉公式求微分方程数值解的程序，并用其与ode23求下列微分方程数值解，对二者作出比较.a ）22,(0)0y x y y '=-=或(0)1y =.b)222()0x y xy x n y '''++-=2()2,()22y y πππ'==-（Bessel 方程，这里令12n =，其精确解为sin y =. c)cos 0,(0)1,(0)0y y x y y '''+===. （2）倒圆锥形容器，上底面直径为1.2m ，容器的高亦为1.2m ，在锥尖的地方开有一直径为3cm 的小孔，容器装满水后，下方小孔开启，由水利学知识可知当水面高度为h 时，水从小孔中流出的速度为v g =为重力加速度，若孔口收缩系数为0.6（即若一个面积单位的小孔向外出水时，水柱截面积为0.6），问水从小孔中流完需多少时间？2分钟时，水面高度是多少？（3）一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸上B 点，已知河水流速1v 与船在静水中的速度2v 之比为k .（a ）建立小船航线的方程，求其解析解.（b ）设12100m,1m/s,2m/s d v v ===，用数值解法求渡河所需时间，任意时刻小船的位置及航行曲线，作图并与解析解比较.（c ）若流速1v 为0,0.5,2 (m/s)，结果将如何？（4）研究种群竞争模型. 当甲、乙两个种群各自生存时，数量演变服从下面规律1212()(1),()(1)x y x t r x y t r y n n =-=- 其中，(),()x t y t 分别为t 时刻甲，乙两个种群的数量，12,r r 为其固有增长率，12,n n 为它们的最大容量，而当这两个种群在同一环境中生存时，由于乙消耗有限资源而对甲的增长产生影响，将甲的方程修改为1112(1)x y x r x s n n =-- （4.22） 这里1s 的含意是：对于供养甲的资源而言，单位数量乙（相对2n ）的消耗率为单位数量甲（相对1n ）消耗的1s 倍，类似地，甲的存在亦影响乙的增长，乙的方程应改为2212()(1)x y y t r y s n n =-- （4.23） 给定种群的初始值为 00(0),(0)x x y y == （4.24）及参数121212,,,,,r r s s n n 后，方程（4.22）与（4.23）确定了两种群的变化规律，因其解析解不存在，试用数值解法研究以下问题：（a ）设121212001,100,0.5,2,10r r n n s s x y ========，计算(),()x t y t ，画出它们的图形及相图(,)x y ，说明时间t 充分大以后(),()x t y t 的变化趋势（人们今天看到的已经是自然界长期演变的结果）.（b ）改变1200,,,r r x y ，但1s 与2s 不变（保持121,1s s <>），分析所得结果，若121.5(1),0.7(1)s s =>=<，再分析结果，由此你得到什么结论，请用各参数生态学上的含义作出解释.（c ）试验当120.8(1),0.7(1)s s =<=<时会有什么结果；当121.5(1), 1.7(1)s s =>=>时又会出现什么结果，能解释这些结果吗？。
