常微分方程模型及其数值解课件.ppt

1、0、导言 在许多实际问题中，例如物理中的速率问题，人口的增长问题，放射性衰变问题，经济学中的边际问题等，常常涉及到两个变量之间的变化规律。微分方程是研究上述问题的一种机理分析方法，它在科技、工程、生态、环境、人口以及经济管理等领域中有着十分广泛的应用。在应用微分方程解决实际问题时，必须经过两个阶段。一是微分方程的建立，建立一个微分方程的实质就是构建函数、自变量以及函数对自变量的导数之间的一种平衡关系。而正确地构建这种平衡关系，需要对实际问题的深入浅出的刻画，根据物理的和非物理的原理、定律或定理，作出合理的假设和简化并将它升华成数学问题。另一个是方程的求解和结果分析。对一些常系数的或特殊函数形式

2、的微分方程，往往能得到解析解，这对实际问题的分析和应用都是有利的,但是大多数变系数的、非线性函数形式的微分方程都是求不出解析解的,此时就需要应用求解微分方程的另一个重要方法数值解法。本章简要介绍有关微分方程模型的概念，微分方程的数值解法和图解法，主要介绍若干建模实例，通过它们展示微分方程模型的建模步骤及解决实际问题的全过程。1、实例及其数学模型、实例及其数学模型 例例1 海上缉私海上缉私 问题问题 海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船正以速度a向正北方向行驶，缉私艇立即以最大速度b前往拦截。用雷达进行跟踪时，可保持缉私艇的速度方向始终指向走私船。建立任意时刻缉私艇的位置和缉私

3、艇航线的数学模型，讨论缉私艇能够追上走私船的条件，求出追上的时间。建立直角坐标系如图，设在t=0时刻缉私艇发现走私船，此时缉私艇的位置在(0, 0)，走私船的位置在(c, 0)。走私船以速度a平行于y轴正向行驶，缉私艇以速度b按指向走私船的方向行驶。在任意时刻t缉私艇位于P(x, y)点，而走私船到达Q(c, at)点，直线PQ与缉私艇航线（图中曲线）相切，切线与x轴正向夹角为。Q(c,at)P(x,y)R(c,y )0yxc缉私艇在x, y方向的速度分别为 ，由直角三角形PQR写出sin 和cos 的表达式，得到微分方程： （1） 初始条件为 （2） 这就是缉私艇位置(x(t), y(t)的

4、数学模型。但是由方程（1）无法得到x(t), y(t)的解析解，需要用数值算法求解。我们将在后面继续讨论这个问题。 sin,cosbdtdybdtdx2222)()()()()()(yatxcyatbdtdyyatxcxcbdtdx0)0(, 0)0(yx例例2 弱肉强食弱肉强食问题问题 自然界中在同一环境下的两个种群之间存在着几种不同的生存方式，比如相互竞争，即争夺同样的食物资源，造成一个种群趋于灭绝，而另一个趋向环境资源容许的最大容量；或者相互依存，即彼此提供部分食物资源，二者和平共处，趋于一种平衡状态；再有一种关系可称之为弱肉强食，即某个种群甲靠丰富的自然资源生存，而另一种群乙靠捕食种群

5、甲为生，种群甲称为食饵（Prey），种群乙为捕食者（Predator），二者组成食饵-捕食者系统。海洋中的食用鱼和软骨鱼（鲨鱼等）、美洲兔和山猫、落叶松和蚜虫等都是这种生存方式的典型。这样两个种群的数量是如何演变的呢？近百年来许多数学家和生态学家对这一系统进行了深入的研究，建立了一系列数学模型，本节介绍的是最初的、最简单的一个模型，它是意大利数学家Volterra在上个世纪20年代建立的。模型模型 用x(t)表示时刻t食饵（如食用鱼）的密度，即一定区域内的数量，y(t)表示捕食者（如鲨鱼）的密度。假设食饵独立生存时的（相对）增长率为常数r0，即 ，而捕食者的存在使食饵的增长率减小，设减小量与捕

6、食者密度成正比，比例系数为a0，则 。 捕食者离开食饵无法生存，设它独自存在时死亡率为常数d0，即 ，而食饵的存在为捕食者提供了食物，使捕食者的死亡率减小，设减小量与食饵密度成正比，比例系数为b0，则 ，实际上，当bxd时捕食者密度将增长。 给定食饵和捕食者密度的初始值x0, y0，由上可知x(t), y(t)满足以下方程： （3）（3）的解x(t), y(t)描述了食饵和捕食者密度随时间的演变过程。但是我们同样得不到x(t), y(t)的解析解，需要用数值算法求解。我们将在3继续讨论这个问题 00)0(,)0()()(yyxxbxydyybxdyaxyrxxayrxrxx/ ayrxx/dy

7、y/)(/bxdyy2 欧拉方法和龙格欧拉方法和龙格库塔方法库塔方法一阶常微分方程初值问题的一般形式为 y=(x,y) ,axb(4)y(a)=其中(x,y)是已知函数,为给定的初值. 如果函数(x,y)在区域axb,-y0为LipschitzLipschitz常数常数,则初值问题(4)有唯一解.yxyyLyxfyxf,),(),( 所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值. a=x0 x1x2xnxN=b其中剖分节点xn=a+nh,n=0,1,N, h称为剖分步长.数值解法就是求精确解y(x)在剖分节点xn上的近似值yny(xn), n=1,2,n

8、. 假设初值问题(4)的解y=y(x)唯一存在且足够光滑.对求解区域a,b做剖分 我们采用数值积分方法来建立差分公式. 2.1 2.1 构造数值解法的基本思想构造数值解法的基本思想 在区间xn,xn+1上对方程(4)做积分,则有对右边的积分应用左矩形公式，则有)5()(,()()(11nnxxnndxxyxfxyxy梯形公式oxyab左矩形公式y=(x)babfafabdxxf)()(2)(baafabdxxf)()()(右矩形公式babfabdxxf)()()(中矩形公式babafabdxxf)2()()(对右边的积分应用左矩形公式，则有)6()(,()()(11nnxxnndxxyxfxy

9、xy因此，建立节点处近似值yn满足的差分公式称之为EulerEuler公式公式. 称为梯形公式梯形公式. )(,()()(1nnnnxyxhfxyxy),(1nnnnyxhfyy1,2 , 1 , 0,0Nny 若对(6)式右边的积分应用梯形求积公式，则可导出差分公式1,2 , 1 , 0,0Nny),(),(2111nnnnnnyxfyxfhyy 利用Euler方法求初值问题 解解 此时的Euler公式为称为EulerEuler中点公式中点公式或称双步双步EulerEuler公式公式. . 若在区间xn-1,xn+1上对方程(4)做积分,则有11)(,()()(11nnxxnndxxyxfx

10、yxy对右边的积分应用中矩形求积公式，则得差分公式),(211nnnnyxhfyy1,2 , 1 , 0,0Nny例例320 ,21122xyxy 0)0(y的数值解.此问题的精确解是y(x)=x/(1+x2).分别取步长h=0.2 ,0.1 ,0.05,计算结果如下)211(221nnnnyxhyy2 , 1 , 0,00nyhxnyny(xn)y(xn)-ynh=0.20.000.400.801.201.602.000.000000.376310.542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00

11、000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.201.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.49

12、1800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.00227Euler中点公式则不然, 计算yn+1时需用到前两步的值yn , yn-1 ,称其为两步方法两步方法,两步以上的方法统称为多步法多步法. 在Euler公式和梯形公式中,为求得yn+1,只需用到前一步的值yn,这种差分方法称为单步法单步法,这是一种自开始方法. 隐式公式中,每次计算yn+1都需解方程,要比显式公式需要更多的计算量,但其计算稳定性较好. 在Euler公式和Euler中点公式中,需要计算的yn+1已被显式表示出来,称这类差分公式为显式公式显式公式, ,而梯

13、形公式中,需要计算的yn+1隐含在等式两侧,称其为隐式公式隐式公式. 从数值积分的角度来看,梯形公式计算数值解的精度要比Euler公式好,但它属于隐式公式,不便于计算. 实际上,常将Euler公式与梯形公式结合使用: 2.2 改进的改进的Euler方法方法),(),(2111nnnnnnyxfyxfhyy1,2 , 1 , 0,0Nny),(01nnnnyxhfyy),(),(21111knnnnnknyxfyxfhyy1,2 , 1 , 0,0Nny 由迭代法收敛的角度看，当 (是给定的精度要求)时, 取 就可以保证迭代公式收敛, 而当h很小时, 收敛是很快的. 而且, 只要|111knkn

14、yy.111knnyy, 12Lyfh),(1nnnnyxhfyy),(),(2111nnnnnnyxfyxfhyy1,2 , 1 , 0,0Nny 通常采用只迭代一次的算法:称之为改进的改进的Euler方法方法. 这是一种单步显式方法. 改进的Euler方法也可以写成)(2211KKhyynn),(1nnyxfK 1,2 , 1 , 0,0Nny y=y-2x/y , 0 x1的数值解, 取步长h=0.1 . 精确解为y(x)=(1+2x)1/2.),(12hKyhxfKnn例例4 求初值问题 y(0)=1 解解 (1) 利用Euler方法nnnnyxyy/2 . 01 . 119 ,2 ,

15、 1 , 0,10ny)(05. 0211KKyynnnnnyxyK/219 ,2 , 1 , 0,10ny计算结果如下:1121 . 0) 1 . 0(21 . 0KyxKyKnnn (2) 利用改进Euler方法nxnEuler方法yn改进Euler法yn精确解y(xn)01234567891000.10.20.30.40.50.60.70.80.9111.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.7177791.78477011.0959091.1840961.2662011.3433601.4164021.48

16、59561.5525151.6164761.6781681.73786911.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.732051 在节点xn+1的误差y(xn+1)-yn+1 ,不仅与yn+1这一步计算有关,而且与前n步计算值yn,yn-1,y1都有关. 为了简化误差的分析,着重研究进行一步计算时产生的误差.即假设yn=y(xn),求误差y(xn+1)-yn+1,这时的误差称为局部截断误差局部截断误差,它可以反映出差分公式的精度.2.3 差分公式的误差分析差分公式的误差分析 如果单步差分公

17、式的局部截断误差为O(hp+1),则称该公式为p p阶方法阶方法.这里p为非负整数.显然,阶数越高,方法的精度越高. 研究差分公式阶的重要手段是Taylor展开式,一元函数和二元函数的Taylor展开式为:另外,在yn=y(xn)的条件下,考虑到y(x)=(x,y(x),则有 321! 3)(! 2)()()()()(hxyhxyhxyxyhxyxynnnnnn2222222),(),(),(! 21),(),(),(),(kyyxfhkyxyxfhxyxfkyyxfhxyxfyxfkyhxfnnnnnnnnnnnnnn y(xn)=(xn,y(xn)=(xn,yn)=n y(xn)=(xn,

18、y(xn)=x(xn,yn)+y(xn,yn)(xn,yn)nnnfyfxfnnnnnnnnnnfyfyfxffyffyxfxfxy2222222)(2)( yn+1=yn+h(xn,yn) 对Euler方法，有 21! 2)()()()()(hxyhxyxyhxyxynnnnn =yn+(xn,yn)h+O(h2)从而有: y(xn+1)-yn+1=O(h2)所以Euler方法是一阶方法.再看改进Euler方法, 因为),(12hKyhxfKnn1hKyfhxffnnn)(21321222122222hOKhyfKhyxfhxfnnn可得所以, 改进的Euler方法是二阶方法.而nnnnnn

19、fyxfhhfyy221)(44222223hOfyffyxfxfhnnnnn)(! 3)(2)()()()(4321hOhxyhxyhxyxyxynnnnn nnnnnfyfxfhhfy22nnnnnnnnnfyfyfxffyffyxfxfh22222223)(26从而有: y(xn+1)-yn+1=O(h3)2.4 Taylor展开方法展开方法 设y(x)是初值问题（4)的精确解, 利用Taylor展开式可得称之为p阶Taylor展开方法. 1)1()(21)!1()(!)(! 2)()()()( pppnpnnnnhpyhPxyhxyhxyxyxy)()(,(!)(,(! 2)(,()(

20、1)1()1(2pnnppnnnnnhOxyxfPhxyxfhxyxhfxy因此，可建立节点处近似值yn满足的差分公式),(!),(! 2),()1()1(21nnppnnnnnnyxfPhyxfhyxhfyy1,2 , 1 , 0,0Nny),(),(),(),()1(yxfyyxfxyxfyxffyfyfxffyffyxfxfyxf2222222)2()(2),(其中所以,此差分公式是p阶方法. 由于Taylor展开方法涉及很多复合函数(x,y(x)的导数的计算,比较繁琐,因而很少直接使用,经常用它为多步方法提供初始值.然而, Taylor展开方法给出了一种构造单步显式高阶方法的途径. E

21、uler方法可写为 可见,公式的局部截断误差为: y(xn+1)-yn+1=O(hp+1).2.5 Runge-Kutta方法方法hKyynn1),(nnyxfK 构造差分公式 改进的Euler方法可写为)(2211KKhyynn),(1nnyxfK ),(12hKyhxfKnn)(22111ppnnKKKhyy),(1nnyxfK ),(12122KhyhxfKnn),(11piipinPnPKhyhxfK其中i,i,ij为待定参数. 若此公式的局部截断误差为由于 yn+1=yn+h1n+h2(n+hxn+hn yn)+O(h3)O(h3),称此公式为p p阶阶Runge-kuttaRung

22、e-kutta方法方法,简称p p阶阶R-KR-K方法方法. 对于p=2的情形, 应有)(22111KKhyynn),(1nnyxfK ),()7(12hKyhxfKnn =yn+h(1+2)n+h22(xn+n yn)+O(h3)(2)(321hOfffhhfyxyynnxnnnn所以,只要令 1+2=1, 2=1/2, 2=1/2 (8) 一般地, 参数由(8)确定的一族差分公式(7)统称为二二阶阶R-KR-K方法方法. .称之为中点公式中点公式,或可写为若取=1,则得1=2=1/2,=1,此时公式(7)就是改进的Euler公式; 若取1=0,则得2=1,=1/2,公式(7)为21hKyy

23、nn),(1nnyxfK ),(121212hKyhxfKnn),(,(21211nnnnnnyxhfyhxhfyy 高阶R-K公式可类似推导. 下面列出常用的三阶、四阶R-K公式. 四阶标准四阶标准R-KR-K公式公式 三阶三阶R-KR-K公式公式)4(63211KKKhyynn),(1nnyxfK )2,(213hKhKyhxfKnn)22(643211KKKKhyynn),(1nnyxfK ),(121212hKyhxfKnn),(34hkyhxfKnn),(121212hKyhxfKnn),(221213hKyhxfKnn 解解 四阶标准R-K公式为)22(4321611KKKKhyy

24、nnnnnyxyK/21例例3 用四阶标准R-K方法求初值问题 y=y-2x/y , 0 x1 y(0)=1的数值解, 取步长h=0.2 .)/()2(2212213hKyhxhKyKnnn)/()2(1211212hKyhxhKyKnnn)/()(2334hKyhxhKyKnnn计算结果如下:nxnyny(xn)nxnyny(xn)0120.00.20.41.001.18321.34171.001.18321.34163450.60.81.01.48331.61251.73211.48321.61251.7321 也可以构造隐式R-K方法,其一般形式为prrrnnKhyy11prKhyhxf

25、Kpiirinrnr, 2 , 1,),(1称之为p p级隐式级隐式R-KR-K方法方法,同显式R-K方法一样确定参数.如)(21211KKhyynn),(1nnyxfK ),(2211212hKhKyhxfKnn是二级二阶隐式R-K方法,也就是梯形公式.但是p级隐式R-K方法的阶可以大于p,例如,一级隐式中点公式为11hKyynn),(121211hKyhxfKnn或写为)(,(121211nnnnnyyhxhfyy它是二阶方法.2.6 变步长变步长Runge-Kutta方法方法 以p阶R-K方法为例讨论.设从xn以步长h计算y(xn+1)的近似值为)(1hny ,局部截断误差为1)(11)

26、(phnnChyxy其中,C是与h无关的常数. 如果将步长减半,取h/2为步长, 从xn经两步计算得到y(xn+1)的近似值记为 ,其局部截断误差为于是有从而,得到事后误差估计)(12hny11)(1121)2(2)(2pppnnChhCyxyh可见,当phnnnnyxyyxyh21)()()(11)(112)(121)()(1)(1)(1122hnnpnnyyyxyhh|)(1)(12hnnyyh 成立时,可取)(112)(hnnyxy .否则,应将步长再次减半进行计算. 求解初值问题的单步显式方法可一统一写为如下形式 yn+1=yn+h(xn,yn,h) （9) 对于Euler方法,有2.

27、7 单步方法的收敛性单步方法的收敛性 y=(x,y) ,axb y(a)= 其中(x,y,h)称为增量函数增量函数. (x,y,h)=(x,y)对于改进的Euler方法,有 (x,y,h)=1/2(x,y)+(x+h,y+h(x,y) 设y(x)是初值问题（4)的解 ,yn是单步法 （9)产生的近似解.如果对任意固定的点xn,均有y(xn),则称单步法（9)是收敛的. 可见,若方法（9)是收敛的,则当h0时,整体截断误差en=y(xn)-yn将趋于零. 定理定理 设单步方法（9)是p1阶方法, 增量函数(x,y,h)在区域axb,-yn)的变化均不超过 ,则称此差分方法是绝对稳定绝对稳定的.

28、讨论数值方法的稳定性,通常仅限于典型的试验方程 y=y 其中是复数且Re()0. 在复平面上,当方法稳定时要求变量h的取值范围称为方法的绝对稳定域绝对稳定域,它与实轴的交集称为绝对稳定区间绝对稳定区间. . 将Euler方法应用于方程y=y, 得到 设在计算yn时产生误差n,计算值yn=yn+n,则n将对以后各节点值计算产生影响.记ym=ym+m ,mn,由上式可知误差m满足方程 m=(1+h)m-1=(1+h)m-nn , mn 对隐式单步方法也可类似讨论.如将梯形公式用于方程y=y,则有 yn+1=yn+h/2 (yn+yn+1) yn+1=(1+h)yn 可见,若要|m|n|,必须且只须

29、|1+h|1 ,因此Euler法的绝对稳定域为|1+h|1,绝对稳定区间是-2Re()h0.解出yn+1得 nnyhhy2121111类似前面分析,可知绝对稳定区域为1112121hh由于Re()0,所以此不等式对任意步长h恒成立,这是隐式公式的优点. 一些常用方法的绝对稳定区间为方 法方法的阶数稳 定 区 间Euler方法梯形方法改进Euler方法二阶R-K方法三阶R-K方法四阶R-K方法122234(-2 , 0)(- , 0)(-2 , 0)(-2 , 0)(-2.51 , 0)(-2.78 , 0) 解解 因y0=1,计算得y10=1024,而y(1)=9.35762310-14.例例

30、4 考虑初值问题 y=-30y , 0 x1 y(0)=1取步长h=0.1 ,利用Euler方法计算y10y(1). y(x)=e-30 x 这是因为h=-3不属于Euler方法的绝对稳定区间. 若取h=0.01,计算得y100=3.23447710-16. 若取h=0.001,计算得y1000=5.91199810-14. 若取h=0.0001,计算得y10000=8.94505710-14. 若取h=0.00001,计算得y100000=9.315610-14. 单步显式方法的稳定性与步长密切相关, 在一种步长下是稳定的差分公式,取大一点步长就可能是不稳定的. 收敛性是反映差分公式本身的截

31、断误差对数值解的影响;稳定性是反映计算过程中舍入误差对数值解的影响.只有即收敛又稳定的差分公式才有实用价值.2.9 线性多步方法线性多步方法 由于在计算yn+1时 ,已经知道yn ,yn-1 ,及(xn,yn), (xn-1,yn-1),利用这些值构造出精度高、计算量小的差分公式就是线性多步法.2.9.1 利用待定参数法构造线性多步方法利用待定参数法构造线性多步方法 r+1步线性多步方法的一般形式为当-10时,公式为隐式公式,反之为显式公式.参数i,i的选择原则是使方法的局部截断误差为 y(xn+1)-yn+1=O(h)r+2 选取参数,0,1,2,使三步方法 yn+1=yn+h(0n+1n-

32、1+2n-2) 这里,局部截断误差是指 ,在yn-i=y(xn-i),i=0,1,r的前提下,误差y(xn+1)-yn+1.为三阶方法. ririiniininfhyy011例例5 解解 设yn=y(xn),yn-1=y(xn-1),yn-2=y(xn-2),则有 n=(xn,y(xn)=y(xn) y(xn+1)=y(xn)+hy(xn)+1/2h2y(xn)+1/6h3y(xn) 于是有若使: y(xn+1)-yn+1=O(h4) ,只要,0,1,2满足: n-1=(xn-1,y(xn-1)=y(xn-1)=y(xn-h) =y(xn)-hy(xn)+1/2h2y(xn)-1/6h3y(4

33、)(xn)+O(h4) n-2=y(xn)-2hy(xn)+2h2y(xn)-4/3h3y(4)(xn)+O(h4) yn+1=y(xn)+h(0+1+2)y(xn)-h2(1+22)y(xn) +h3(1/21+22)y(xn)-h4/6(1+82)y(4)(xn)+O(h5) +1/24h4y(4)(xn)+O(h5) =1, 0+1+2=1, 1+22=-1/2 , 1+42=1/3于是有三步三阶显式差分公式设pr(x)是函数(x,y(x)的某个r次插值多项式,则有解之得: yn+1=yn+h/12(23n-16n-1+5n-2) 因为 125,34,1223, 12102.9.2 利用

34、数值积分构造线性多步方法利用数值积分构造线性多步方法 1)(,()()(1nnxxnndxxyxfxyxy1)()()(1nnxxnrnnRdxxpxyxy其中 选取不同的插值多项式pr(x),就可导出不同的差分公式.下面介绍常用的AdamsAdams公式公式. . 设已求得精确解y(x)在步长为h的等距节点xn-r,xn上的近似值yn-r ,yn , 记k=(xk,yk) ,利用r+1个数据(xn-r,n-r),(xn,n)构造r次Lagrange插值多项式由此,可建立差分公式 1.Adams 1.Adams显式公式显式公式 1)(1nnxxrnndxxpyy1)()(,(nnxxrndxx

35、pxyxfR其中 rjjnjnrfxlxp0)()(由此,可建立差分公式 由于 rjxxxxxlrjkkknjnknjn, 1 , 0)()()(0 hrj jnxxjnrjnnfdxxlyynn1)(01)()()()(110thxxdxxxxxdxxlnxxrjkkknjnknxxjnnnnn令100, 1 , 0,)()!( !) 1(rjkkjrjdtktjrjh则有 rjjnrjnnfhyy01称之为r+1r+1步步AdamsAdams显式公式显式公式. . 下面列出几个带有局部截断误差主项的Adams显式公式 r=0 yn+1= =yn+hn+(1/2)h2y(xn) 2.Adam

36、s 2.Adams隐式公式隐式公式 r=1 yn+1= =yn+(h/2)(3n-n-1)+(5/12)h3y(xn) r=2 yn+1= =yn+(h/12)(23n-16n-1+5n-2) +(3/8)h4y(4)(xn) r=3 yn+1= =yn+(h/24)(55n-59n-1+37n-2-9n-3) +(251/720)h5y(5)(xn) 如果利用r+1个数据(xn-r+1,n-r+1),(xn+1,n+1)构造r次Lagrange插值多项式pr(x),则可导出数值稳定性好的隐式公式,称为AdamsAdams隐式公式隐式公式,其一般形式为其中系数为 010*, 1 , 0,)()

37、!( !) 1(rjkkjrjrjdtktjrjrjjnrjnnfhyy01*1下面列出几个带有局部截断误差主项的Adams隐式公式 r=0 yn+1= =yn+hn+1-(1/2)h2y(xn) r=1 yn+1= =yn+(h/2)(n+n+1)-(1/12)h3y(xn) r=2 yn+1= =yn+(h/12)(5n+1+8n-n-1) -(1/24)h4y(4)(xn) r=3 yn+1= =yn+(h/24)(9n+1+19n-5n-1+n-2) -(19/720)h5y(5)(xn) 3.Adams 3.Adams预估预估- -校正公式校正公式 由显式公式提供一个预估值，再用隐式

38、公式校正一次，求得数值解，称为预估校正方法预估校正方法。 校正 yn+1= =yn+(h/24)(9n+1+19n-5n-1+n-2) 一般预估公式和校正公式都采用同阶公式。例如： 预估 yn+1= =yn+(h/24)(55n-59n-1+37n-2-9n-3) n+1=(xn+1,yn+1) , n=3,4,称为四阶Adams预估校正公式.实际计算时通常用四阶单步方法(如四阶R-K公式)为它提供起始值y1,y2,y3 . 例例6 用四阶Adams预估校正公式求解初值问题 y=y-2x/y , 0 x1 y(0)=1取步长h=0.1. 解 用四阶R-K公式提供起始值,计算结果如下xnR-k法

39、yn预估值yn校正值yn精确值y(xn)00.10.20.30.40.50.60.70.80.9111.0954461.1832171.2649121.3415511.4140451.4830171.5489171.6121141.6729141.7315661.3416411.4142131.4832391.5491921.6124501.6733181.73204811.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.7320513 RK方法的方法的MATLAB实现实现对于微分方程（组）的初值问

40、题 龙格库塔方法可用如下MATLAB命令实现其计算：t，x=ode23(f，ts，x0，options)t，x=ode45(f，ts，x0，options)其中ode23用的是3级2阶龙格库塔公式，ode45用的是以Runge-Kutta-Fehberg命名的 5级4阶公式。TnTnTnxxxxtxfffxxxxtftx),(,)(),(,),(),()(00100011命令的输入f是待解方程写成的函数m文件：function dx=f(t,x)dx=f1;f2;fn;若ts=t0,t1,t2, ,tf，则输出在指定时刻t0,t1,t2, ,tf的函数值；等分点时用ts=t0:k:tf，输出在

41、t0,tf内等分点处的函数值。x0为函数初值（n维向量）。options可用于设定误差限（options缺省时设定相对误差10-3，绝对误差10-6），命令为： options=odeset(reltol,rt,abstol,at)其中rt，at分别为设定的相对误差和绝对误差。命令的输出t为指定的ts，x为相应的函数值（n维向量）。注意，计算步长是根据误差限自动调整的，并不是输入中指定的ts的分点。下面用MATLAB软件解决1提出的两个问题例例1 海上缉私（续）海上缉私（续）模型的数值解模型的数值解 1. 设a=20 (海里/小时)，b=40 (海里/小时)，c=15 (海里)，由模型（1），

42、（2）求任意时刻缉私艇的位置及缉私艇航线。 对于给出的a, b, c用MATLAB求数值解时，记x(1)=x, x(2)=y, x=(x(1), x(2)T。编写如下m文件：function dx=jisi(t,x) % 建立名为jisi的函数m文件a=20; b=40; c=15;s=sqrt(c-x(1)2+(a*t-x(2)2); dx=b*(c-x(1)/s;b*(a*t-x(2)/s; % 以向量形式表示方程（1）然后运行以下程序：ts=0:0.05:0.5; % 设定t的起终点及中间的等分点，终点可先作试探，再按照x(t)c=15调整到0.5x0=0,0; % 输入x，y的初始值（

43、2）t,x=ode45(jisi,ts,x0); % 调用ode45计算t,x % 输出t, x(t), y(t)plot(t,x), grid, % 按照数值输出作x(t), y(t)的图形gtext(x(t), gtext(y(t), pauseplot(x(:,1),x(:,2), grid, % 作y( x) 的图形gtext(x), gtext(y)得到的数值结果x(t), y(t)为缉私艇的位置，列入表1。走私船的位置记作x1(t), y1(t)，显然x1(t)= c=15，y1(t)=at=20t，将y1(t)列入表1最后一列。可知当t=0.5（小时），x, y与x1, y1几乎

44、一致，认为缉私艇追上走私船。x(t), y(t)及y(x)的图形见图2，y(x)为缉私艇的航线。tx(t)y(t)y1(t)00000.051.99840.06981.00.103.98540.29242.00.155.94450.69063.00.207.85151.28994.00.259.67052.11785.00.3011.34963.20056.00.3512.81704.55527.00.4013.98066.17738.00.4514.74518.02739.00.5015.00469.997910.0a=20， b=40，c=15的数值解x(t), y(t)和y1(t) 00

45、.10.20.30.40.505101520 x(t)y(t)051015200246810 xya=20，b=40，c=15 时x(t), y(t) 和y( x)的图形 2. 设b，c不变，而a变大为30，35，接近40 (海里/小时)，观察解的变化 修改a的输入，并相应地延长t的终点。设a=35，t的终点经试探，调整为1.6合适。表2是计算结果，其中x(t), y(t)有两列数字，左边的是用“缺省”精度（即相对误差10-3，绝对误差10-6）计算的，中间的y1(t) =at=35t是走私船到达的位置。可知t=1.3, 1.4, 1.5时缉私艇的位置x15, 但y与y1(t)相差甚远，t=1

46、.6时x, y与x1, y1也有差距，这是累积误差造成的。可利用ode45的控制参数options提高精度（上面的“调用ode45计算”用以下程序代替），如设 opt=odeset(reltol,1e-6,abstol,1e-9);t,x=ode45(jisi,ts,x0,opt);得到的x(t), y(t)，与走私船到达的位置x1(t)，y1(t) 相对照，可知t=1.6时x, y与x1, y1几乎一致，认为缉私艇追上走私船。x(t), y(t)及y(x)的图形见图，y(x)为缉私艇的航线，当x接近15时航线几乎是正北方向，形成沿走私船逃向的追赶态势。 读者可尝试寻求缉私艇追上走私船的判断标

47、准，设计程序并计算。00.511.520102030405060 x(t)y(t)051015200204060 xy a=35，b=40，c=15 时x(t), y(t) 和y(x)的图形 模型的解析解模型的解析解 要想得到缉私艇拦截到走私船的精确时间和位置，必须对模型（1）做进一步的分析，设法得到某种形式的解析解。 将（1）的两式相除，消去dt得到 （27）（27）式对x求导得 （28）为消去式中的dt，利用曲线弧长的微分，缉私艇的速度 ，以及微分关系，得atydxdyxc )(dxdtadxydxc22)( （29）这是关于y(x)的二阶微分方程，若令 ，可化为p(x)的一阶微分方程 （

48、30）依题意其初始条件为 （31）（30）为可分离变量方程，容易求得在（31）下的解为 （32）222)(1)(dxdybadxydxcdxdyp bakpkdxdpxc/,1)(20)0(p)()(21kkcxccxcp设，即走私船速度a小于缉私艇速度b，对（32）积分并注意到得 （35）这就是缉私艇航线的解析表达式。 由（35）式知，当x=c时， 这也是走私船的y坐标。因为走私船速度是a，所以缉私艇拦截到走私船的时间为2111)(11)(112kkccxckcxckcykk2221ababckkcy221abbct例例2 弱肉强食（续）弱肉强食（续）模型的数值解模型的数值解 为了考察模型（

49、3）描述的食饵和捕食者密度随时间的演变过程，选取模型参数r=1，d=0.5， a=0.1，b=0.02，x0=25，y0=2，用MATLAB求（3）的数值解，记x(1)=x, x(2)=y, x=(x(1), x(2)T。编写如下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; % 以向量形式表示方程（3） 然后运行以下程序：ts=0:0.1:15; % t的终值经试验后定为15。x0=25,2;t,x=ode45(shier,ts,x0);t,xplot(t,x),grid,gtext(fontsize12x(t),gtext(fontsize12y(t), % 将标记x(t) , y(t)的字体放大pause,plot(x(:,1),x(:,2),grid,xlabel(x),ylabel(y)