数值计算9章.doc

第9章常微分方法的数值解法科学技术与工程问题常常需要建立微分方程形式的数学模型,下面是这类问题的例子.设N(t)为某物种的数量,为该物种的出生率与死亡率之差,为生物的食物供给及它们所占空间的限制,描述该物种增长率的数学模型是设Q是电容器上的带电量,C为电容,R为电阻,E为电源的电动势,描述该电容器充电过程的数学模型是以上两个例子是常微分方程的初值问题,下面是一个两点边值问题的例子.设一根长为L的矩形截面的梁,两端固定.E是弹性模量,S是端点作用力,I(x)是惯性矩,q是均匀载荷强度,梁的桡度y(x)满足如下方程针对实际问题建立的数学模型,要找出模型解的解析表达式往往是困难的,甚至是不可能的.因此,需要研究和掌握微分方程的数值解法,即计算域内离散点上解的近似值的方法.本章讨论常微分方程数值解的基本方法和理论.9.1Euler方法9.1.1Euler方法及其有关的方法考虑一阶常微分方程的初值问题(9.1)设f(x,y)是连续函数,对y满足Lipschitz条件,即存在正数L,使得对于任意两点(x,y)与(x,),有|f(x,y)-f(x,)|L|y-|,这样初值问题的解是存在唯一的,而且连续依赖于初始条件.为了求得离散点上的函数值,将微分方程的连续问题(9.1)式进行离散化.一般是引入点列xn,这里xn=xn-1+hn,n=1,2,.称hn为步长,经常考虑定长的情形,即hn=h,xn=x0+nh, n=0,1,.记y(xn)为初始问题(9.1)的准确解y(x)在xn处的值,用均差近似代替(9.1)式中的导数得令yn为y(xn)的近似值,将上面两个近似式写成等式,整理后得yn+1=yn+hf(xn,yn),n=0,1,(9.2)yn+1=yn+hf(xn+1,yn+1),n=0,1,.(9.3)从x0处的初值y0开始,按(9.2)式可逐步计算以后各点上的近似值.称(9.2)式为显式Euler公式.由于(9.3)式的右端隐含有待求函数值yn+1,不能逐步显式计算,称(9.3)式为隐式Euler公式或后退Euler公式.如果将(9.2)式和(9.3)式作算术平均,就得梯形公式yn+1=yn+f(xn,yn)+f(xn+1,yn+1),n=0,1,.(9.4)梯形公式也是隐式公式.以上公式都是由yn去计算yn+1,故称它们为单步法.例9.1取h=0.1,用Euler方法隐式Euler方法和梯形方法解解本题有f(x,y)=x-y+1,y0=1.如果用Euler方法,由(9.2)式并代入h=0.1得yn+1=0.1xn+0.9yn+0.1.同理,用隐式Euler方法有yn+1=(0.1xn+1+yn+0.1).用梯形公式有yn+1=(0.1xn+0.95yn+0.105).3种方法及准确解y(x)=x+e-x的数值结果如表9-1所示.从表中可看到,在xn=0.5处,Euler方法和隐式Euler方法的误差|y(xn)-yn|分别是1.610-2和1.410-2,而梯形方法的误差却是2.510-4.表9-1xnEuler方法隐式Euler方法梯形法准确解011110.11.0000001.0090911.0047621.0048370.21.0100001.0264461.0185941.0187310.31.0290001.0513151.0406331.0408180.41.0561001.0830131.0700961.0703200.51.0904901.1209211.1062781.106531在例9.1中,由于f(x,y)对y是线性的,所以对隐式公式也可方便地计算yn+1.但是,当f(x,y)是y的非线性函数时,如y=5x+,其隐式Euler公式为yn+1=yn+h(5xn+1+).显然,它是yn+1的非线性方程,可以选择非线性方程求根的迭代法求解yn+1.以梯形公式为例,可用显式Euler公式提供迭代初值,用迭代公式反复迭代,直到,其中,步长h成为迭代参数,它需要满足一定的条件迭代公式才能收敛.若将(9.4)式减去该迭代公式,得yn+1-=(f(xn+1,yn+1)-f(xn+1,).假设f(x,y)关于y满足Lipschitz条件,则有yn+1-yn+1-,这里,L是Lipschitz常数.由上式可见,当1即h时,迭代序列收敛于yn+1.对于隐式公式,通常采用预估校正技术,即先用显式公式计算,得到预估值,然后以预估值为隐式公式的迭代初值,用隐式公式迭代一次得到校正值,称为预估校正技术.例如,用显式Euler公式作预估,用梯形公式作校正,即(9.5)称该公式为改进的Euler公式.它等价的显式公式为yn+1=yn+(f(xn,yn)+f(xn+1,yn+hf(xn,yn),(9.6)也可表示为下列平均化的形式例9.2取h=0.1,用改进的Euler方法解解按(9.5)式,改进的Euler公式为.由于y0=1,h=0.1得计算结果如表9-2.该初值问题的准确解为y(x)=.表9-2xn0.10.20.30.40.50.60.70.8yn1.09591.18411.26621.34341.41641.48601.55251.6153y(xn)1.09541.18321.26491.34161.41421.48321.54921.61659.1.2局部误差和方法的阶初值问题(9.1)式的单步法可以写成如下统一形式yn+1=yn+h(xn,xn+1,yn,yn+1,h),(9.7)其中与f有关.若中不含yn+1,则方法是显式的,否则是隐式的,所以一般显式单步法表示为yn+1=yn+h(xn,yn,h).(9.8)例如,Euler方法中,有(x,y,h)=f(x,y).对于不同的方法,计算值yn与准确解y(xn)的误差各不相同.所以有必要讨论方法的截断误差 我们称en=y(xn)-yn为某一方法在xn点的整体截断误差.显然,en不单与xn这步的计算有关,它与以前各步的计算都有关,所以误差称为整体的.分析和估计整体截断误差en是复杂的.为此,我们假设xn处的yn没有误差,即yn=y(xn),考虑从xn到xn+1这一步的误差,这就是如下的局部误差的概念.定义9.1设y(x)是初值问题(9.1)式的准确解,则称Tn+1=y(xn+1)-y(xn)-h(xn,xn+1,y(xn),y(xn+1),h)为单步法(9.7)式的局部截断误差.定义9.2如果给定方法的局部截断误差Tn+1=O(hp+1),其中p1为整数,则称该方法是p阶的,或者具有p阶精度.若一个p阶单步法的局部截断误差为Tn+1=g(xn,y(xn)hp+1+O(hp+2),则称其第一个非零项g(xn,y(xn)hp+1为该方法的局部截断误差的主项.对于Euler方法,由Taylor展开式有Tn+1=y(xn+1)-y(xn)-hf(xn,y(xn)=y(xn+1)-y(xn)-hy(xn)=y(xn)+y(xn)+O(h4)=O(h2).所以,Euler方法是一种一阶方法,其局部截断误差的主项为y(xn).对于隐式Euler方法,其局部截断误差为Tn+1=y(xn+1)-y(xn)-hf(xn+1,y(xn+1)=y(xn+1)-y(xn)-hy(xn+1)=-y(xn)+O(h3)=O(h2).所以,隐式Euler方法也是一种一阶方法,该方法的局部截断误差的主项为-y(xn),仅与显式Euler方法的局部截断误差的主项相差一个符号.梯形方法也是一种隐式单步法,类似可得其局部截断误差Tn+1=y(xn+1)-y(xn)-(f(xn,y(xn)+f(xn+1,y(xn+1)=-y(xn)+O(h4)=O(h3).可见,梯形方法是2阶精度的.例9.3证明yn+1=yn+(4f(xn,yn)+2f(xn+1,yn+1)+hf(xn,yn)定义的隐式单步法公式是3阶的.证设yn=y(xn),则由f(xn,yn)=y(xn),f(xn,yn)=y(xn),f(xn+1,yn+1)=y(xn+1)=y(xn)+hy(xn)+y(xn)+y(4)(xn)+O(h4),可得yn+1=y(xn)+(4f(xn,yn)+2f(xn+1,yn+1)+hf(xn,yn)=y(xn)+hy(xn)+y(xn)+y(xn)+y(4)(xn)+O(h5).将此式与y(xn+1)在xn 处的Taylor展开式比较得y(xn+1)-yn+1=-172h4y(4)(xn)+O(h5)=O(h4).由此可见,所给公式是3阶的.9.2Runge-Kutta方法9.2.1Runge-Kutta方法的基本思想显式Euler方法是最简单的单步法,它是一阶的,它可以看作Taylor展开后取前两项.因此,得到高阶方法的一个直接想法是用Taylor展开,如果能计算y(x)的高阶导数,则可写出p阶方法的计算公式yn+1=yn+hyn+yn+,其中是y(j)(xn)的近似值,j=0,1,2,p.若将f(x,y),分别记成f,fx,fy,则对于2阶和3阶导数可表示为y=fx+fyf,y=fxx+2fxyf+fxfy+fyyf2+fy2f.这个方法并不实用,因为一般情况下,求f(x,y)的导数相当麻烦.从计算高阶导数的公式知道,方法的截断误差提高一阶,需要增加的计算量很大.例9.4取步长h=0.25,用2阶Taylor展开法求初值问题y=x2+y2,y(1)=1的解在x=1.5时的近似值.解2阶Taylor展开公式为y(xn+1)=y(xn)+hy(xn)+y(xn)+O(h3).用y=x2+y2,y=2x+2yy=2x+2y(x2+y2)代入上式并略去高阶项O(h3),则得求解公式yn+1=yn+h(xn2+yn2)+(2xn+2yn(xn2+yn2).由y(1)=y0=1,计算得y(1.25)y1=1.6875,y(1.50)y2=3.333298.尽管如此,Taylor展开法启发我们用区间上若干个点的导数f,而不是高阶导数,将它们作线性组合得到平均斜率,将其与解的Taylor展开式相比较,使前面若干项吻合,从而得到具有一定阶的方法.这就是Runge-Kutta方法的基本思想,其一般形式为(9.9)其中K1=f(xn,yn),Ki=f(xn+cih,yn+cih),i=2,3,L,ci1,=1,=1.它的局部截断误差是Tn+1=y(xn+1)-y(xn)-,(9.10)其中,K*i与Ki的区别在于:用微分方程准确解y(xn)代替Ki中的yn就得到K*i.参数i,ci和aij待定,确定它们的原则和方法是:将(9.10)式中的y(xn+1)在xn处作Taylor展开,将K*i在(xn,y(xn)处作二元Taylor展开,将展开式按h的幂次整理后,令Tn+1中h的低次幂的系数为零,使Tn+1首项中h的幂次尽量高,比如使Tn+1=O(hp+1),则称(9.9)式为L级p阶Runge-Kutta方法(简称R-K法).类似于显式R-K公式(9.9),稍加改变,就得到隐式R-K公式yn+1=yn+hLi=1iKi,其中Ki=f(xn+cih,yn+cih),i=1,2,L.它与显式R-K公式的区别在于:显式公式中,对系数aij求和的上限是i-1,从而aij构成的矩阵是一个严格下三角阵.而在隐式公式中,对系数aij求和的上限是L,从而aij构成的矩阵是方阵,需要用迭代法求出近似斜率Ki(i=1,2,L).推导隐式公式的思路和方法与显式公式类似.9.2.2几类显式Runge-Kutta方法对于L=2,则K1=f(xn,yn),K2=f(xn+c2h,yn+c2hK1),yn+1=yn+h(1K1+2K2),其局部截断误差是Tn+1=y(xn+1)-y(xn)-h(1K*1+2K*2).(9.11)将Tn+1中的各项作Taylor展开,并利用y(xn)=f(xn,y(xn),y=fx+fyf,则有y(xn+1)=y(xn)+hy(xn)+y(xn)+y(xn)+O(h4),K*1=f(xn,y(xn)=y(xn),K*2=f(xn+c2h,y(xn)+c2hy(xn)=y(xn)+c2hy(xn)+(fxx+2fxyf+fyyf2)+O(h3),将它们代入(9.11)式,整理后得Tn+1=(1-1-2)hy(xn)+12-2c2h2y(xn) +O(h4).选取1,2和c2,使方法的阶尽可能高,就是使h和h2的系数为零,因为h3的系数一般不为零.于是得到方程组1+2=1,2c2=.显然,该方程组有无穷多组解,从而得到一族二级二阶R-K方法.若以c2为自由参数,取c2=得中点公式yn+1=yn+hf(xn+,yn+f(xn,yn),(9.12)取c2=得Heun公式yn+1=yn+(f(xn,yn)+3f(xn+h,yn+hf(xn,yn),(9.13)取c2=1得改进的Euler公式(9.6).对于L=3的情形,要计算3个斜率的近似值:K1=f(xn,yn),K2=f(xn+c2h,yn+c2hK1),K3=f(xn+c3h,yn+c3h(a31K1+a32K2).类似于二阶方法的推导,可以得三阶的方法,所得系数应满足的方程组是1+2+3=1,a21=1,2c2+3c3=,2c22+3c23=,3c2c3a32=,a31+a32=1.该方程组的解也是不唯一的.常见的一种三级三阶方法是K1=f(xn,yn),K2=f(xn+h2,yn+h2K1),K3=f(xn+h,yn-hK1+2hK2),yn+1=yn+(K1+4K2+K3).对于L=4的情况,可进行类似推导.最常用的四级四阶方法是如下的经典R-K方法K1=f(xn,yn),K2=f(xn+,yn+K1),K3=f(xn+,yn+K2),(9.14)K4=f(xn+h,yn+hK3),yn+1=yn+(K1+2K2+2K3+K4).为了分析经典R-K公式的计算量和计算精度,将四阶经典R-K公式(9.14)与一阶显式Euler公式(9.2)及2阶改进的Euler公式(9.6)相比较.一般来说,公式的级数越大,计算右端项f的次数越多,计算量越大.在同样步长的情况下,Euler方法每步只计算一个函数值,而经典方法要计算4个函数值.4阶R-K法的计算量差不多是改进的Euler公式的2倍,是显式Euler公式的4倍.下面的例子中Euler方法用步长h1,2阶改进的Euler法用步长2h1,而四阶的经典公式用步长4h1.这样,从xn到xn+4h1,3种方法都计算了4个函数值,计算量大体相当.例9.5考虑初值问题y=-y+1,y(0)=0.其解析解为y(x)=1-e-x.分别用h=0.025的显式Euler方法,h=0.05的改进Euler法和h=0.1的经典R-K方法计算到x=0.5.3种方法在x方向每前进0.1都要计算4个右端函数值,计算量相当.计算结果列于表9-3.从计算结果看,在工作量大致相同的情况下,还是经典R-K方法比其他两种方法的结果好得多.在x=0.5处,3种方法的误差分别是3.810-3,1.310-4,2.810-7.经典R-K方法对多数好条件问题(fy,我们反复将步长折半进行计算,直到为止,这时取最终得到的作为结果;如果为止,这时再将前一次步长折半的结果作为所要的结果.这种通过加倍或折半处理步长的方法称作变步长方法.虽然为了选择步长,每一步的计算量有所增加,但总体考虑还是值得的.9.3单步法的收敛性和稳定性9.3.1单步法的收敛性数值解法的基本思想就是要通过某种离散化方法,将微分方程转化为某种差分方程(例如(9.8)式)来求解.这种转化是否合理,还要看差分方程的解yn是否收敛到微分方程的准确解y(xn).定义9.3对于任意固定的xn=x0+nh,若对于初值问题(9.1)的显式单步法(9.8)产生的近似解yn,均有yny(xn)(h0,同时n),则称该方法是收敛的.在定义中,xn是固定的点,h0时有n,n不是固定的.显然,若方法是收敛的,则在固定点xn处的整体截断误差en=y(xn)-yn趋于零.下面给出方法收敛的条件.定理9.1设初值问题(9.1)的单步法(9.8)是p阶的(p1),且函数满足对y的Lipschitz条件,即存在常数L0,使(x,y1,h)-(x,y2,h)L|y1-y2|,对一切y1,y2成立,则方法(9.8)收敛,且y(xn)-yn=O(hp).证仍记en=y(xn)-yn,根据局部截断误差的定义y(xn+1)=y(xn)+h(xn,y(xn),h)+Tn+1.将此式与(9.8)式相减得en+1=en+h(xn,y(xn),h)-(xn,yn,h)+Tn+1.因为(9.8)式是p阶的,所以存在h0,当0hh0时有|Tn+1|chp+1.再用的Lipschitz条件有|en+1|en|+hL|en|+chp+1.为了方便,记=1+hL,=chp+1,即有|en+1|en|+.由此可推得|en|en-1|+2|en-2|+n|e0|+(n-1+n-2+1).利用关系式eLh=1+Lh+1+Lh,n=(1+Lh)nenLh=eL(xn-x0),可以得到|en|e0|eL(xn-x0)+(eL(xn-x0)-1)chpL-1.现在取y0=y(x0),有e0=0,于是en=O(hp).定理得证.容易证明,如果(9.1)式的f满足对y的Lipschitz条件,且初值是准确的,则显式Euler法改进的Euler法和R-K方法是收敛的.由定理9.1说明,f关于y满足Lipschitz条件是使单步收敛的充分条件,而且,还说明一个方法的整体截断误差比局部截断误差低一阶.所以,常常通过求出局部截断误差去了解整体截断误差的大小.单步法的显式形式(9.8)可写成(xn,yn,h)=.(9.15)称(xn,yn,h)为增量函数.对于收敛的方法,固定x=xn,有yny(xn)(h0),从而y(xn)(h0).对于(9.15)式,我们自然要考虑(xn,yn,h)f(xn,y(xn)(h0)是否成立.这就是相容性问题.定义9.4若方法(9.8)的增量函数满足(x,y,0)=f(x,y),则称方法(9.8)与初值问题(9.1)是相容的.相容性说明数值计算的差分方程(9.15)趋于(9.1)式中微分方程.我们本章讨论的数值方法都是与原初值问题相容的.9.3.2单步法的稳定性对于一种收敛的相容的差分方程,由于计算过程中舍入误差总会存在,我们需要讨论其数值稳定性.一个不稳定的差分方程会使计算解失真或计算失败.为了讨论方便起见,将(9.1)式中的f(x,y)在解域内某一点(a,b)作Taylor展开并局部线性化,即y=f(x,y)f(a,b)+(x-a)fx(a,b)+(y-b)fy(a,b)=fy(a,b)y+c1x+c2.令=fy(a,b)u=y+,利用线性化的关系,可得uu.因此,我们通过如下的试验方程y=y(9.16)讨论数值方法的稳定性.当某一步yn有舍入误差时,若以后的计算中不会逐步扩大,则称这种稳定性为绝对稳定性.现讨论显式Euler法的稳定性.将显式Euler法用于试验方程(9.16),有yn+1=(1+h)yn.当yn有舍入误差时,其近似值为,从而有=(1+h).令n=yn-,得到误差传播方程n+1=(1+h)n.令E(h)=1+h,只要|E(h)|1,则显式Euler方法的解和误差都不会恶性发展,即-2h0时,显式Euler方法是稳定的,即是条件稳定的.对于梯形方法,应用于试验方程后,有,同理,有误差方程n+1=E(h)n,其中E(h)=.因此当0时,梯形方法是稳定的.一般地,在试验方程(9.16)中,我们只考虑0的情形,我们认为微分方程是不稳定的.比如,将显式Euler方法用于(9.1)式中的方程,有n+1=n+h(f(xn,yn)-f(xn,)=(1+hfy(xn,)n.当fy(xn,)0时,有1+hfy(xn,)1.对于每一种单步法应用于试验方程(9.16),可得yn+1=E()yn,(9.17)然而,对于不同的单步法,E()有不同的表达式.定义9.5若(9.17)式中的|E()|1,则称对应的单步法是绝对稳定的.在复平面上,满足|E()|1的区域,称为方法的绝对稳定区域,它与实轴的交称为绝对稳定区间.一些单步法的E()表达式和它们的绝对稳定区间列于表9-4.从表中可见,隐式方法比显式方法的绝对稳定性好.表9-4方法E()绝对稳定区间Euler法1+h-2h0改进的Euler法-2h0三阶R-K法-2.51h0四阶R-K法-2.785h0隐式Euler法-h-h0梯形式-0例9.6分别取h=1,2,4,用经典R-K方法求解y=-y+x-e-1,y(1)=0,其准确解为y(x)=e-x+x-1-e-1.解本题=-1,h分别为-1,-2,-4.由表94可知,当h2.785时,该方法才稳定.计算结果列于表9-5.表9-5xnh=1的解h=2的解h=4的解准确解53.63943.67305.47153.638997.63237.636716.82917.63221311.632111.632657.617111.6321由表95可见,h=1和h=2时,计算结果确实稳定,h=4,结果发散.此外,h为1的计算精度比h为2的计算精度高.因为h越小,方法的截断误差越小.但若h过分小的话,计算步数非常多,其累积误差会增加.所以,实际计算时,应选取合适的步长,常常采用自动变步长的R-K方法.9.4线性多步法常微分方程初值问题(9.1)的数值解法中,除了像Runge-Kutta型公式等单步法之外,还有另一种类型的解法,即某一步解的公式不仅与前一步解的值有关,而且与前若干步解的值有关,利用前面多步的信息预测下一步的值,这就是多步法的基本思想,可能期望获得较高的精度.构造多步法有多种途径,下面先讨论基于数值积分的方法.9.4.1基于数值积分的方法将(9.1)式中的方程在区间xn,xn+1上积分,可以得到y(xn+1)=y(xn)+.(9.18)如推导Newton-Cotes求积公式一样,用等距节点的插值多项式来替代被积函数,再对插值多项式积分,这样就得到一系列求积公式.例如,用梯形方法计算积分项(f(xn,y(xn)+f(xn+1,y(xn+1),代入(9.18)式有y(xn+1)y(xn)+(f(xn,y(xn)+f(xn+1,y(xn+1).据此即可导出(9.4)式.一般地,设由r+1个数据点(xn,fn),(xn-1,fn-1),(xn-r,fn-r)构造插值多项式Pr(x),这里fk=f(xk,yk),xk=x0+kh.运用插值公式有将(9.18)式离散化即得下列计算公式yn+1=yn+h,(9.19)其中rj=,j=0,1,r.由此可得(9.19)式中的系数,其具体数值见表96.(9.19)式是一个r+1步的显式公式,称为Adams,显式公式.r=0时,即为Euler公式.表9-6j012340j121j3-1122j23-165243j55-5937-97204j1901-27742616-1274251在上述Adams,显式公式的推导中,选用了xn,xn-1,xn-r作为插值节点.这样的插值多项式Pr(x)在求积区间xn,xn+1上逼近f(x,y(x)是一个外推结果.为了改善逼近效果,我们变外推为内插,即改用xn+1,xn,xn-r+1为插值节点,用数据点(xn+1,fn+1),(xn,fn),(xn-r+1,fn-r+1)构造插值多项式Pr(x),则有于是,我们有如下的计算公式yn+1=yn+h, (9.20)其中rj=,j=0,1,r,其具体数值见表97.(9.20)式是隐式公式,称为Adams,隐式公式.r=0,1时,分别为隐式Euler公式和梯形公式.表9-7j012340j121j11122j58-1243j919-517204j251646-264106-19对于隐式公式(9.20),需要用迭代求解.确定yn+1的迭代公式为迭代收敛条件为h|r0|L1,其中L为f关于y的Lipschitz常数.利用插值多项式的余项,可以求出Adams,方法的局部截断误差.当然也可以从得到的显式和隐式Adams,公式,由局部截断误差的定义来求出方法的局部截断误差.表98中列出了它们的局部截断误差的主项,由表98可以看出,Adams,隐式方法的局部截断误差要小.表9-8r0123Adams,显式公式Adams,隐式公式9.4.2基于Taylor展开的方法基于数值积分可以构造出一系列求解常微分方程的数值计算公式,计算公式由插值多项式唯一确定.下面介绍基于Taylor展开的待定系数法,它可灵活地构造出线性多步法.对固定的步数,可以选取待定系数使线性多步法的阶尽可能高.还可以根据需要,确定显式还是隐式.设构造如下具有p阶精度的线性多步公式yn+1=0yn+1yn-1+ryn-r+h(-1fn+1+0fn+rfn-r).(9.21)当-1=0时,则(9.21)式为显式多步式.当-10时,(9.21)式为隐式多步式.它们的局部截断误差为利用原微分方程,有(9.22)现利用Taylor展开定理,确定线性多步公式(9.21)中的待定参数k,k,使它达到p阶精度,即Tn+1=O(hp+1).对(9.22)式的右端各项在xn点处作Taylor展开有将它们代入(9.22)式整理后得使y(xn),h,h2,hp的系数为零,得到关于k和k的线性方程组 (9.23)而且得到线性多步法的局部截断误差当参数k和k满足(9.23)式时,线性多步法公式(9.21)达到p阶精度,满足(9.23)式的k和k可能有多组解.下面我们构造几个著名的四阶线性多步公式,考虑下列形式的公式(9.24)(9.25)由于r=3,p=4,由(9.23)式得到5个方程,而(9.24)式中有9个未知量,因此,(9.24)式中有4个自由度.若取-1=0,1=2=3=0,由(9.23)式得到其他5个待定参数的方程组,解之得0=1,0=,1=,2=,3=.代入(9.24)式和(9.25)式,得到常用的四步四阶显式Adams,公式和它的余项(9.26)(9.27)若取-1=0,0=1=2=0,由(9.23)式得到其他5个待定参数的方程组,解之得3=1,0=,1=,2=,3=0.由此构造成著名的四步四阶显式Milne公式和它的余项yn+1=yn-3+h(2fn-fn-1+2fn-2),(9.28)Tn+1=h5y(5)(xn)+O(h6).(9.29)若取1=2=3=0,3=0,由(9.23)式得到其他5个待定参数的方程组,从而得三步四阶隐式Adams,公式及其余项yn+1=yn+(9fn+1+19fn-5fn-1+fn-2),(9.30)Tn+1=h5y(5)(xn)+O(h6)(9.31)若取1=3=0,2=3=0,求解(9.23)式得著名的三步四阶隐式Hamming公式及其余项yn+1=(9yn-yn-2)+(fn+1+2fn-fn-1),(9.32)Tn+1=h5y(5)(xn)+O(h6).(9.33)若取0=2=3=0,3=0,求解(9.23)式得到隐式Simpson公式及其余项yn+1=yn-1+(fn+1+4fn+fn-1),Tn+1=h5y(5)(xn)+O(h6).例9.7分别取h=0.2,2,用4阶显式Milne公式和4阶隐式Hamming公式求解例9.6所给的初值问题.解我们用单步法提供多步法的开始值.由4阶经典R-K公式为Milne公式提供开始值y1,y2,y3,为Hamming公式提供开始值y1,y2.h=0.2和h=2时的计算结果及与准确解之间的误差分别列于表9-9和表9-10.表9-9(h=0.2)xnMilne方法误差Hamming方法误差2.20.94294268-1.810-50.942919554.210-62.41.122833495.010-61.122833864.610-62.61.30643214-3.810-51.306389304.810-62.81.492916251.410-51.492925824.810-63.01.68195450-4.710-51.681902994.610-6表9-10(h=2) xnMilne方法误差Hamming方法误差75.645745-1.310-25.645745-1.310-297.3823252.510-17.637126-4.910-31110.905316-1.39.635636-3.510-3134.1438317.511.632261-1.410-41558.310717-4.510113.632240-1.110-317-249.6626722.710215.6316904.310-4从表99看出,两种多步法的计算精度都很高,Hamming公式比Milne公式更精确.这是因为Hamming公式的截断误差主项的系数比Milne公式小.从表910看到,当计算步长变大后,显式多步法Milne公式的计算结果误差增大,不稳定,而隐式多步法Hamming公式的计算结果仍然是稳定的,这说明隐式公式的稳定性比同阶的显式公式好.经典R-K法和上述四阶线性多步公式都是4阶精度,但每前进一步,前者要计算4次微分方程的右端函数值,而后者只要计算一次新的右端函数值,计算量减小了.9.4.3预估校正算法显式多步法容易计算,但其精度和稳定性没有相应的隐式方法好.然而,隐式多步法需解方程,如果初值选得不当,则计算量较大.因此,设法选取好的迭代初值是必要的.初值的自然选取是采用同阶显式多步法计算得到的解作为隐式方法迭代的初值.这样,迭代次数不会多.若只迭代一次,则这样的算法就是预估校正算法.对于线性多步法,常用的预估校正方法有四阶Adams,显隐式预估校正公式和Milne-Hamming方法.(1) Adams,预估校正公式.由(9.26)式作为预估公式,由(9.30)式作为校正公式,构成Adams,预估校正公式ypn+1=yn+(55fn-59fn-1+37fn-2-9fn-3),yn+1=yn+(9f(xn+1, ypn+1)+19fn-5fn-1+fn-2).若需作进一步的修正,则记上式所得的yn+1= ycn+1由(9.27)式和(9.31)式有Tpn+1=y(xn+1)- ypn+1h5y(5)(xn), Tcn+1=y(xn+1)- ycn+1h5y(5)(xn).于是得到y(xn+1)- ypn+1-( ypn+1- ycn+1),y(xn+1)- ycn+1(ypn+1- ycn+1).由此可见,若记= ypn+1+(ycn+1- ypn+1),= ycn+1-(ycn+1- ypn+1),则,分别比ypn+1, ycn+1更好.但注意到,的表达式中, ycn+1是未知的,因此改为= ypn+1+(ycn-ypn).这样,得到下面修正的Adams,预估校正公式预估: ypn+1=yn+(55fn-59fn-1+37fn-2-9fn-3),修正:= ypn+1+(ycn-ypn),校正: ycn+1=yn+(9f(xn+1,)+19fn-5fn-1+fn-2),修正: yn+1= ycn+1-(ycn+1- ypn+1).在计算时,可调节计算步长h,使|-(ycn+1- ypn+1)|,其中是要求达到的计算精度.开始值y1,y2,y3由同阶单步法提供,当计算y4时,可取.(2) 修正Hamming公式.将Milne公式(9.28)和Hamming公式(9.32)结合,构成Milne-Hamming预估校正公式ypn+1=yn-3+h(2fn-fn-1+2fn-2),yn+1=(9yn-yn-2)+(f(xn+1,ypn+1)+2fn-fn-1).若需作进一步的修正,则记上式所得的yn+1=ycn+1,由(9.29)式和(9.33)式有Tpn+1=y(xn+1)- ypn+1h5y(5)(xn),Tcn+1=y(xn+1)- ycn+1-h5y(5)(xn).于是得到y(xn+1)- ypn+1(ycn+1- ypn+1),y(xn+1)- ycn+1-(ycn+1- ypn+1).由此分别得Milne和Hamming公式的修正公式= ypn+1+(ycn+1- ypn+1),= ycn+1-(ycn+1- ypn+1).从而构成如下的修正Hamming公式预估: ypn+1=yn-3+h(2fn-fn-1+2fn-2),修正:= ypn+1+(ycn-ypn),校正: ycn+1=yn+(9yn-yn-2)+(f(xn+1,)+2fn-fn-1),修正: yn+1= ycn+1-(ycn+1- ypn+1).在计算时,可调节计算步长h,使|-(ycn+1- ypn+1)|.开始值y1,y2由同阶单步法提供,当计算y3时,可取.例9.8取h=0.2,用Milne-Hamming预估校正公式和修正Hamming公式求解例9.6所给的初值问题.解用经典R-K法提供开始值,计算结果列于表911.将表99与表911所示的计算结果进行比较,它们的计算精度排列次序是:修正Hamming公式的精度最好,其次是隐式Hamming公式,再次是Milne-Hamming预估校正公式,最后是Milne公式.表9-11 xnMilne-Hamming误差修正Hamming误差2.20.942916257.510-60.94292449-7.810-72.41.122828729.810-61.12283955-1.010-62.61.306382711.110-51.30639537-1.210-62.81.492918161.210-51.49293184-1.210-63.01.681894671.310-51.68190879-1.210-69.5一阶方程组的数值解法9.5.1一阶方程组和高阶方程考虑一阶常微分方程组的初值问题(9.34)若将其中的未知函数方程的右端项都表示成向量形式y=(y1,y2,yN)T,f=(f1,f2,fN)T,初始条件表示成y(x0)=y0=(y01,y02,y0N)T,那么,(9.34)式可以写成 (9.35)可见,(9.35)式在形式上与一个方程的初值问题一样.关于一个方程的初值问题的数值方法均适用于方程组.相应的理论问题也可类似地讨论.下面仅写出两种数值方法作说明.梯形方法yn+1=yn+(f(xn,yn)+f(xn+1,yn+1),或表达为yn+1,i=yni+(fi(xn,yn)+fi(xn+1,yn+1),i=1,2,N,其中yni是第i个因变量yi(x)在节点xn处的近似值,相应地,fi(xn,yn)=fi(xn,yn1,yn2,ynN).经典R-K方法yn+1=yn+(K1+2K2+2K3+K4),其中K1=f(xn,yn),K2=f,K3=f,K4=f(xn+h,yn+hK3),或表达为yn+1,i=yni+(K1i+2K2i+2K3i+K4i),i=1,2,N,其中K1i=fi(xn,yn1,yn2,ynN),K2i=fi(xn+,yn1+K11,y n2+K12,ynN+K1N),K3i=fi(xn+,yn1+K21,yn2+K22,ynN+K2N),K4i=fi(xn+h,yn1+hK31,yn2+hK32,ynN+hK3N).对于高阶方程,可把它转化为一阶方程组.例如,考察下列m阶微分方程 (9.36)只要引进新的变量y1=y,y2=y,ym=y(m-1),则可将m阶方程(9.36)化为如下的一阶方程组(9.37)因此,可用求解方程组形式的方法来求解(9.37).例如,对2阶微分方程初值问题y=f(x,y,y),y(x0)=y0,y(x0)=y0,令z=y,则可将该初值问题化为一阶微分方程组的初值问题y=z,z=f(x,y,z),y(x0)=y0,z(x0)=y0. 对此应用Euler方法有,即得数值计算公式yn+1=yn+hyn, yn+1=yn+nf(xn,yn,yn),n=0,1,2,.9.5.2刚性方程组先考虑两个简单的初值问题.问题1.问题2.这两个问题有同样的解.采用4阶经典R-K方法来计算上面的两个问题,以相同的误差要求来自动选取步长,计算从x=0到x=10.第一个问题可用相当大的步长,而第二个问题能使用的步长小到难以接受.如果改用某种低阶隐式公式,那么这两个问题均可用较大的步长,计算出大致符合要求的解来.上述显示出来的现象称为刚性.问题2是刚性的,问题1是非刚性的.由于这两个问题的解是相同的,因此这种现象不是问题的解的作用,而是方程组的一种特性所引起的.基于这个事实,较为正确的应称之为刚性方程组而不是刚性问题.考虑方程组的通解.对于问题1,方程组的系数矩阵的特征值为1=-1和2=-3,其通解为,(9.38)其中1,2