第五讲第4章线性多步法(续)

第4章线性多步法4.1线性多步法的一般公式前面给出了求解初值问题(1.2.1)的单步法，其特点是计算时只用到的值，此时的值均已算出.如果在计算时除用的值外，还用到的值，这就是多步法.若记,h为步长，，则线性多步法可表示为(4.1.1)其中为常数，若(即不同时为零)，称(4.1.1)为线性k步法.计算时用到前面已算出的k个值.当时，(4.1.1)为显式多步方法，当则称(4.1.1)为隐式多步法.隐式方法与梯形方法一样，计算时要用迭代法求.多步法(4.1.1)的局部截断误差定义也与单步法类似.举例来说，对于初值问题'1,yyx(0)1y，步数k=2时，线性多步法表示为101111011(),1,2,nnnnnnyyyhfffn当时，格式为显示的：10110111[(1)(1],1,2,nnnnnnnyyyhyxyxn，而时，格式为隐式的：10111110111[(1)(1)(1],1,2,nnnnnnnnnyyyhyxyxyxn。定义4.1设y(x)是初值问题(1.2.1)的精确解，线性多步法(4.1.1)在处的局部截断误差定义为(4.1.2)若，则称线性多步法(4.1.1)是p阶的.如果我们希望得到的多步法是p阶的，则可利用Taylor公式展开，将在处展开到阶，它可表示为(4.1.3)注意，(4.1.2)式按Taylor展开可得经整理比较系数可得(4.1.4)若线性多步法(4.1.1)为p阶，则可令于是得局部截断误差(4.1.5)右端第一项称为局部截断误差主项.称为误差常数.要使多步法(4.1.1)逼近初值问题(1.2.1),方法的阶p≥1,当p=1时，则，由(4.1.4)得(4.1.6)称为相容性条件.公式(4.1.1)当k=1时即为单步法，若,由(4.1.6)则得式(4.1.1)就是，即为Euler法.此时，方法为p=1阶.若，由得,为确定及，必须令,由(4.1.4)得及此时(4.1.1)就是即为梯形法.由故p=2,方法是二阶的，与3.1节中给出的结果相同.实际上，当k给定后，则可利用(4.1.4)求出公式(4.1.1)中的系数及，并求得的表达式(4.1.5).4.2Adams显式与隐式方法形如(4.2.1)的k步法称为Adams方法，当时为Adams显式方法，当时，称为Adams隐式方法.对初值问题(1.2.1)的方程两端从到积分得显然只要对右端的积分用插值求积公式，求积节点取为即可推出形如(4.2.1)的多步法，但这里我们仍采用Taylor展开的方法直接确定(4.2.1)的系数.对比(4.1.1)可知，此时，只要确定即可.现在若k=4且，即为4步的Adams显式方法其中为待定参数，若直接用(4.1.4)，可知此时自然成立，再令可得解此方程组得05524，15924，23724，3924.由此得到于是得到四阶Adams显式方法及其余项为(4.2.2)(4.2.3)若，则可得到p=4的Adams隐式公式，则k=3并令，由(4.1.4)可得解得，而，于是得到四阶Adams隐式方法及余项为(4.2.4)(4.2.5)一般情形，k步Adams显式方法是k阶的，k=1即为Euler法，k=2为k=3时，.k步隐式方法是(k+1)阶公式，k=1为梯形法，k=2为三阶隐式Adams公式k步的Adams方法计算时必须先用其他方法求出前面k个初值才能按给定公式算出后面各点的值，它每步只需计算一个新的f值，计算量少，但改变步长时前面的也要跟着重算，不如单步法简便.例4.1用四阶显式Adams方法及四阶隐式Adams方法解初值问题,步长h=0.1用到的初始值由精确解计算得到.解本题直接由公式(4.2.2)及(4.2.4)计算得到.对于显式方法，将直接代入式(4.2.2)得到其中.对于隐式方法，由式(4.2.4)可得到直接求出，而不用迭代，得到计算结果如表所示.表4-1Adams方法和Adams隐式方法的数值解与精确解比较4.3Adams预测-校正方法上述给出的Adams显式方法计算简单，但精度比隐式方法差，而隐式方法由于每步要做迭代，计算不方便.为了避免迭代，通常可将同阶的显式Adams方法与隐式Adams方法结合，组成预测-校正方法.以四阶方法为例，可用显式方法(4.2.2)计算初始近似，这个步骤称为预测(Predictor)，以P表示，接着计算f值(Evaluation),，这个步骤用E表示，然后用隐式公式(4.2.4)计算，称为校正(Corrector)，以C表示，最后再计算，为下一步计算做准备.整个算法如下：(4.3.1)公式(4.3.1)称为四阶Adams预测-校正方法(PECE).其matlab程序如下functiony=DEYCJZ_adms(f,h,a,b,y0,varvec,type)formatlong;N=(b-a)/h;y=zeros(N+1,1);x=a:h:b;y(1)=y0;y(2)=y0+h*Funval(f,varvec,[x(1)y(1)]);y(3)=y(2)+h*Funval(f,varvec,[x(2)y(2)]);y(4)=y(3)+h*Funval(f,varvec,[x(3)y(3)]);fori=5:N+1v1=Funval(f,varvec,[x(i-4)y(i-4)]);v2=Funval(f,varvec,[x(i-3)y(i-3)]);v3=Funval(f,varvec,[x(i-2)y(i-2)]);v4=Funval(f,varvec,[x(i-1)y(i-1)]);t=y(i-1)+h*(55*v4-59*v3+37*v2-9*v1)/24;ft=Funval(f,varvec,[x(i)t]);y(i)=y(i-1)+h*(9*ft+19*v4-5*v3+v2)/24;endformatshort;问题1（4.3.1）中的第四步在程序中那一行实现了？利用(4.2.2)和(4.2.4)的局部截断误差(4.2.3)和(4.2.5)可对预测-校正方法(4.3.1)进行修改，在(4.3.1)中的步骤P有5(5)11251()()720pnnnyxyhyx720问题2这个约等式是怎么得来的？什么方法？对于步骤C有5(5)1119()720pnnnyyhyx两式相减可得5(5)11720()()270pnnnhyxyy于是有若用代替上式,并令显然比更好，但注意到的表达式中是未知的，因此改为下面给出修正的预测-校正格式(PMECME).(4.3.2)经过修正后的PMECME格式比原来PECE格式提高一阶.问题3试着编出该程序！！4.4Milne方法与Hamming方法与Adams显式方法不同的另一类四阶显式方法的计算公式形如(4.4.1)这里为待定常数，此公式也是k=4步方法，即计算时要用到4个值.为了确定，当然可以利用公式(4.2.1)直接算出，但下面我们直接利用Taylor展开式确定,,使它的阶尽量高.方法(4.4.1)的局部截断误差为将它在点展成Taylor级数，得要使公式的阶尽量高，要令前3项系数为0.即解得,代入公式，的系数为0，故(4.4.2)于是得四阶方法(4.4.3)称为Milne公式，它的局部截断误差为(4.4.2).与(4.4.3)配对的隐式方法为k=3的多步法，它的一般形式可表示为要求公式的阶p=4，可直接用(4.2.1)，并令,可得(4.4.4)若令，可解出，于是得到下列四阶方法(4.4.5)称为Simpson公式，它的局部截断误差为(4.4.6)用Simpson公式与Milne公式(4.4.3)相匹配，用(4.4.3)做预测，(4.4.5)做校正，由于(4.4.5)的稳定性较差，因此通常较少使用.为了改善稳定性，可重新选择四阶的隐式公式，Hamming通过试验，发现在(4.4.4)中若令,得到的公式稳定性较好，此时(4.4.4)的解为,于是得四阶多步法(4.4.7)称为Hamming公式，它的局部截断误差为(4.4.8)用Milne公式(4.4.3)与Hamming公式(4.4.7)相匹配，并利用截断误差公式(4.4.2)与(4.4.8)改进计算结果.131212114(22)31[93(2)]8nnnnnnnnnnnyyfffyyyhfff(4.4.7’)该算法称为Hamming预测-校正法。类似Adams预测-校正格式(4.3.2)，可得以下的修正的milne-Hamming预测-校正格式：(4.4.9)附hamming程序functiony=DEYCJZ_hm(f,h,a,b,y0,varvec,type)formatlong;N=(b-a)/h;y=zeros(N+1,1);x=a:h:b;y(1)=y0;y(2)=y0+h*Funval(f,varvec,[x(1)y(1)]);y(3)=y(2)+h*Funval(f,varvec,[x(2)y(2)]);y(4)=y(3)+h*Funval(f,varvec,[x(3)y(3)]);iftype==1%hamming预测校正法fori=5:N+1v1=Funval(f,varvec,[x(i-3)y(i-3)]);v2=Funval(f,varvec,[x(i-2)y(i-2)]);v3=Funval(f,varvec,[x(i-1)y(i-1)]);t=y(i-4)+4*h*(2*v3-v2+2*v1)/3;ft=Funval(f,varvec,[x(i)t]);y(i)=(9*y(i-1)-y(i-3)+3*h*(2*v3+ft-v2))/8;endelse%修正的hamming预测校正法p0=0;c=0;fori=5:N+1v1=Funval(f,varvec,[x(i-3)y(i-3)]);v2=Funval(f,varvec,[x(i-2)y(i-2)]);v3=Funval(f,varvec,[x(i-1)y(i-1)]);p=y(i-4)+4*h*(2*v3-v2+2*v1)/3;M=p-112*(p0-c)/121;F=Funval(f,varvec,[x(i),M]);c=(9*y(i-1)-y(i-3)+3*h*(2*v3+F-v2))/8;y(i)=c+9*(p-c)/121;p0=p;endendformatshort;例4.2用四步四阶显式Milne公式及三步四阶隐式Hamming公式解初值问题,步长h=0.1初值仍由精确解给出，要求计算到为止，给出计算结果及误差，并与例4.1结果比较.解直接用公式(4.4.3)及(4.4.7)计算.用Milne法计算公式为其中误差用Hamming方法(4.4.7)计算公式为可解得，n=2,3,4误差从所得结果可见Milne方法误差比显式Adams方法误差略小，而Hamming方法与隐式Adams方法误差相当.例4.3将例4.2的初值问题用修正的Milne-Hamming预测-校正公式计算及，初值，仍用已算出的精确解，即，给出计算结果及误差.解根据修正的Milne-Hamming预测-校正公式(4.4.9)得从结果看，此方法误差比四阶Adams隐式法和四阶Hamming方法小，这与理论分析一致.讲解：线性多步法（4.1.1）的局部截断误差定义为与单步法相似，可表示为（4.1.2），即只要直接将右端各项在处展成Taylor公式，根据公式阶数为阶，即按的幂整理，令各项系数为0，则可求得相应的线性多步法及其局部截断误差，这里只用到一元函数的Taylor展开.因此不必记系数满足的公式（7.5.4），只要直接展开即可，它不但可以求出Adams显式与隐式公式以及Milne公式，Hamming公式等，还可以求出任何需要的多步法公式，下面再给出两个例题，说明如何直接用Taylor展开的方法.例4.4解初值问题用显式二步