1、1线性多步法线性多步法2011.10.25liukean微分方程数值解哈工大计算数学23线性多步法线性多步法 3-1线性多步法的一般形式线性多步法的一般形式 对于方程对于方程 y=f(x,y)x a,b y(a)=y0 (1)等分分划等分分划a,b,h=(b-a)/n,xi=a+ih,i=1,n 想法：由于想法：由于y=f(x,y)形式已知，故可以形式已知，故可以构造构造 yk+1=a0yk+a1yk-1+apyk-p +h(b-1yk+1+b0yk+b1yk-1+bpyk-p)(2)来逼近来逼近 y(xk+1)微分方程数值解哈工大计算数学3说明：1）ap,bp不同时为零；不同时为零；2）b-
2、1=0,（2）式称为显格式；）式称为显格式；b-10,存在存在hb-1yk+1=hb-1f(xk+1,yk+1)项，称项，称为隐格式；为隐格式；3）（）（2）式是已知）式是已知y0,yp,1+p个点的解值，个点的解值，外推点外推点yk+1的值，外推过程；的值，外推过程；4）实际上，（）实际上，（2）式是用）式是用p+1阶差分方程阶差分方程逼近微分方程（逼近微分方程（1）；差分方程的理论一点）；差分方程的理论一点也不比微分方程容易，但是，有了初始值，也不比微分方程容易，但是，有了初始值，可以递推计算。可以递推计算。微分方程数值解哈工大计算数学43-2误差与精度误差与精度 定义定义1 代数精度代数
3、精度(用试验多项式定义的用试验多项式定义的)定义定义2 设设yk=y(xk)是（是（1）的解，代入（）的解，代入（2），），称称 Tk+1=y(xk+1)-yk+1 =y(xk+1)a0yk+a1yk-1+apyk-p +h(b-1yk+1+b0yk+b1yk-1+bpyk-p)为（为（2）式从）式从xk到到xk+1一步的局部误差。一步的局部误差。k=p，p+1，n（3）微分方程数值解哈工大计算数学521()(1)2()()()()()()2!()()()()0()!(1)!k ikkkpppppkkihy xy xih y xy xihihyxyxhpp 21()(1)1()()()()()
4、()2!()()()()0()!(1)!k ikkkpppppkkihy xy xih y xy xihihyxyxhpp 欲使欲使Tk+1最小，将最小，将y(xk-i)，y(xk-i)，i=-1,0，p在在xk处作处作Taylor展开展开微分方程数值解哈工大计算数学6 代入（代入（3）式，）式，按按h 的幂次合并同类项、整理后得到的幂次合并同类项、整理后得到 101()111111(1)21(1)()1()()()!1()(1)!(1)()()0()rkikijprrjjjiikjiiprpiirpppikiTa y xhi aib yxjhiappi b yxh 微分方程数值解哈工大计算数
5、学7 代入（代入（3）式，）式，按按h 的幂次合并同类项、整理后得到的幂次合并同类项、整理后得到 Tk+1=c0y(xk)+c1h y(xk)+c2h2 y(xk)+cqhq y（q）(xk)+其中其中0011111111().11()()!piippiiiippqqqiiiicaci abciaibq （5）（4）(2)(2)式的精度取决于式的精度取决于ci是否为零！是否为零！微分方程数值解哈工大计算数学8 若若c0=c1=cr=0，cr+10 则则Tk+1=cr+1hr+1 y(r+1)(xk)+cr+2hr+2 y(r+2)(xk)+可见可见y(x)Pr x时，时，Tk+1 0，而当而当
6、y(x)=xr+1 时，时，Tk+1=cr+1hr+1(r+1)！0 即（即（2）式具有）式具有r 阶代数精度阶代数精度,具有具有r阶计算阶计算精度。精度。对于低于对于低于r 次的多项式而言，次的多项式而言，当然有当然有误差阶主部误差阶主部主部系数主部系数()1()0iiiki rc h yx 0微分方程数值解哈工大计算数学9 定义定义3 称称 c0=c1=0，的线性多步法（的线性多步法（2）式是相容的。）式是相容的。综上：1）线性多步法的构造，无外乎是在（6）成立的情况下，构造精度更高的算法，即确定ai，bi问题。2）后面的任务：在精度确定的条件下，去确定ai，bi，由（5）式知道问题变成了
7、线性方程组求根问题。算法给出后，要讨论其：误差估计、收敛性、稳定性。00111()piippiiiiai ab（6 6）微分方程数值解哈工大计算数学103-3 线性多步法的构造方法线性多步法的构造方法 3-3-1Taylor公式方法公式方法 令令c0=c1=cr=0，由（，由（5）得到关于）得到关于ai，bi的的2p+3个待定参数的线性方程组（个待定参数的线性方程组（r+1个个方程）。方程）。可见可见r=2p+2时，（时，（7）式解存在、唯一。）式解存在、唯一。得到得到p+1步方法。步方法。00111111010()1.0()()1piippiiiippqqriiiicaci abciaib（
8、7）微分方程数值解哈工大计算数学11 实践中，一般取实践中，一般取r a0=1 a0=1 c1=0=b-1+b0=1 b-1=1/2 c2=0=2b-1=1 b0=1/2 得到 yk+1=yk+(h/2)(f(xk+1,yk+1)+f(xk,yk)梯形隐格式，p+1单步法、隐格式，r=2阶。c3=-1/12 微分方程数值解哈工大计算数学13例题6-3 P=1 yk+1=a0yk+a1yk-1+h(b-1yk+1+b0yk+b1yk-1)确定确定a1为自由参数。令为自由参数。令c0=c1=c2=c3=0得到线得到线性代数方程组性代数方程组 1-(a0+a1)=0 a0=1-a1 1+a1-(b-
9、1+b0+b1)=0 b-1=(5-a1)/12 1-a1-2(b-1-b1)=0 b0=2(1+a1)/3 1+a1-3(b-1+b1)=0 b1=(5a1-1)/12=c4=(a1-1)/24 c5=-(1+a1)/180 取取a1=1,c5=-1/190,得到得到4阶阶Simpson方法方法 yk+1=yk-1+(yk+1+4yk+yk-1)h/3微分方程数值解哈工大计算数学14 类似地类似地,可以得到可以得到Milne格式格式 yk+1=yk-3+(2yk+yk-1+2yk-2)4h/3 Tk=h5 y(5)(xk)14/45+显式显式 Hamming方法方法 yk+1=(9yk-yk
10、-2)/8+(yk+1+2yk-yk-1)3h/8 Tk=-h5 y(5)(xk)/40+隐式隐式 都是都是4阶精度，但是误差常数不同，隐格式小于阶精度，但是误差常数不同，隐格式小于显式格式。显式格式。预估预估-校正法校正法:显显,隐式交替使用技术隐式交替使用技术 yk+1=yk-3+(2yk+yk-1+2yk-2)4h/3 yk+1=(9yk-yk-2)/8+(f(xk+1,yk+1)+2yk-yk-1)3h/8微分方程数值解哈工大计算数学153-4线性多步法的收敛性线性多步法的收敛性 定义定义4 线性多步法线性多步法(2)式的初值式的初值yk=yk(h),k=0,1,p 满足满足 其中其中
11、g0是问题是问题(1)的初值，任给的初值，任给x a,b有有 则称（则称（2）式是收敛的。）式是收敛的。00lim(),0,1,khy hg kp0,lim()kkkhkxakhyy x微分方程数值解哈工大计算数学16 问题：问题：1）如何判断（）如何判断（2）式的收敛性？）式的收敛性？2）方程不同，收敛性不同。）方程不同，收敛性不同。讨论的方法适用范围要大。讨论的方法适用范围要大。3）好思想）好思想 选一个试验方程作尺度。选一个试验方程作尺度。方程方程 y=y，y(a)=g0 xa,b 其解析解为其解析解为 y(x)=g0 e(x-a)对试验方程讨论好的方法对试验方程讨论好的方法,对一般方程
12、对一般方程也应该好也应该好!微分方程数值解哈工大计算数学17 用（用（2）式求实验方程的近似解）式求实验方程的近似解 yk+1=a0yk+a1yk-1+apyk-p +h(b-1yk+1+b0yk+b1yk-1+bpyk-p)(2)令令 yk=rk为之解，为之解，(目的：变成多项式根的性质目的：变成多项式根的性质问题问题)则有则有 rk+1=a0rk+a1 rk-1+aprk-p+h(b-1 rk+1+b0 rk+b1 rk-1+bprk-p)令令 k=p rp+1=a0rp+a1 rp-1+ap-1r+ap+h(b-1 rp+1+b0 rp+b1 rp-1+bp-1r+bp)（差分方程变成了
13、（差分方程变成了p+1次多项式）。次多项式）。微分方程数值解哈工大计算数学18 记记(r,h)=rp+1(a0rp+a1 rp-1+ap-1r+ap)h(b-1 rp+1+b0 rp+b1 rp-1+bp-1r+bp)=(r)h(r)(r)=)=rp+1(a0rp+a1 rp-1+ap-1r+ap)(r)=)=b-1 rp+1+b0 rp+b1 rp-1+bp-1r+bp 分别称为（分别称为（2）式的特征多项式、第一、）式的特征多项式、第一、第二特征多项式。第二特征多项式。决定了线性多步法的决定了线性多步法的收敛性、稳定性收敛性、稳定性微分方程数值解哈工大计算数学19 假设特征多项式假设特征多
14、项式(r,h)=rp+1(a0rp+a1 rp-1+ap-1r+ap)h(b-1 rp+1+b0 rp+b1 rp-1+bp-1r+bp)有有p+1个根个根 r0(h),r1(h),rp(h),则其通解可表成则其通解可表成 yk=d0rk0(h)+d1rk1(h)+dprkp(h)定理定理6.2 线性多步法是相容的线性多步法是相容的(c0=c1=0)，则，则由上式确定的根由上式确定的根yk有如下三条性质：有如下三条性质：微分方程数值解哈工大计算数学20 I(r,h)=0有一个根有一个根r0(h)r0(h)=1+h+0(h 2),(h0)II 当当h0时，时，ykg0，k=0,p 则当则当h0时
15、，有时，有 d0y0=g0 di0,i=1,2,p III当当h0时，时，ykg0，k=0,p 则当则当h0,k时，有时，有 d0r0(h)k g0 e(x-a)=y(x)微分方程数值解哈工大计算数学21定理定理2的内涵的内涵 I (r,h)=0的通解的通解 yk=d0rk0(h)+d1rk1(h)+dprkp(h)g0 e(x-a)=y(x)主部,h0,y(a+kh)=g0 ekh,(k),无界II r0(h)=1+h+0(h 2),(h0)蕴涵蕴涵 eh=1+h+0(h 2)=r0(h)Taylor展式展式微分方程数值解哈工大计算数学22III(r,h)=(r)h(r)g0 e(x-a)=
16、y(x)|r0(h)|1稳定稳定y(a+kh)=g0 ekh y*(a+kh)=g0 ek(h)*|g0 ekh-g0 ek(h)*|1+h+0(h 2)-(1+(h)*+0(h 2)|h-(h)*|主部(h0)微分方程数值解哈工大计算数学23 定义定义5 若若(r)=rp+1(a0rp+a1 rp-1+ap-1r+ap)的所有根的模均不大于的所有根的模均不大于1，且模为，且模为1的的根是单根，则称根是单根，则称(r)以及相应的线性以及相应的线性多步法满足根条件。多步法满足根条件。(r)=0的根的根|ri|1;1;|ri|=1,=1,ri是单是单根根.微分方程数值解哈工大计算数学24 定理定理
17、3 线性多步法收敛线性多步法收敛=满足根条件满足根条件 定理定理4 线性多步法相容线性多步法相容 =(1)=0，(1)=(1)定理定理5 线性多步法收敛线性多步法收敛=满足相容条件满足相容条件 定理定理6 线性多步法收敛线性多步法收敛=满足相容满足相容+根条根条件。件。即即c0=c1=0 (1)=0，(1)=(1)=(r)=0的根的根|ri|1;1;|ri|=1,=1,ri是单是单根根 可以计算，使收敛问题变成简单的计算可以计算，使收敛问题变成简单的计算微分方程数值解哈工大计算数学25例题6-4 用线性用线性2步法步法 yk+2-3yk+1+2yk=h(fk+1-2fk)y0=s0(h),y1
18、=s1(h)解初值问题解初值问题y=2x,y(0)=0.注意：注意：(r)=r2-3r+2，(r)=r-2(1)=0,(1)=(1)=-1,所以方法是相所以方法是相容的。容的。但是但是(r)=0的根，的根，r1=1,r2=21 故方法故方法不收敛不收敛。微分方程数值解哈工大计算数学26例题例题6-5 P=1 yk+1=a0yk+a1yk-1+h(b-1yk+1+b0yk+b1yk-1)确定确定a1为自由参数。令为自由参数。令c0=c1=c2=c3=0得到线性代数方得到线性代数方程组程组 1-(a0+a1)=0 a0=1-a1 1+a1-(b-1+b0+b1)=0 b-1=(5-a1)/12 1
19、-a1-2(b-1-b1)=0 b0=2(1+a1)/3 1+a1-3(b-1+b1)=0 b1=(5a1-1)/12=c4=(a1-1)/24 c5=-(1+a1)/180 yk+1=(1-a1)yk+a1yk-1+h(5-a1)/12 yk+1+2(1+a1)/3 yk+(5a1-1)/12 yk-1)问题：问题：a1=？方法方法收敛收敛。微分方程数值解哈工大计算数学27例题例题6-5 yk+1=(1-a1)yk+a1yk-1+h(5-a1)/12 yk+1+2(1+a1)/3 yk+(5a1-1)/12 yk-1)问题：问题：a1=？方法收敛。？方法收敛。注意到注意到(r)=rp+1(a
20、0rp+a1 rp-1+ap-1r+ap)(r)=r2(1-a1)r-a1=(r-1)(r+a1)=0 r 1=1,r 2=-a1 当当-1a11时格式收敛时格式收敛微分方程数值解哈工大计算数学283-5线性多步法的稳定性线性多步法的稳定性 线性多步法的收敛性取决于线性多步法的收敛性取决于(r)=0的根的根|ri|1;1;|ri|=1,=1,ri是单根。在是单根。在线性多步法线性多步法收敛的条件下，是否存在一个区间收敛的条件下，是否存在一个区间ha,b,线性多步法收敛线性多步法收敛?如果有，计算舍入误差只要不出此范围，如果有，计算舍入误差只要不出此范围，线性多步法就收敛。线性多步法就收敛。引入
21、计算稳定性的概念。引入计算稳定性的概念。线性多步法的稳定性取决于线性多步法的稳定性取决于(r,h)=0的根的性质。的根的性质。微分方程数值解哈工大计算数学29 定义定义 5 线性多步法收敛，线性多步法收敛，ri(h),i=0,1,p 是是(r,h)=0的根，的根，r0(h)=1+h+0(h 2),(h0)I 相对稳定区间相对稳定区间 若对于任意的若对于任意的ha,b,有有|ri(h)|r0(h)|,i=1,p 且当且当|ri(h)|r0(h)|,时，时，ri(h)是单根是单根 则称方法在则称方法在a,b上为相对稳定的。称上为相对稳定的。称a,b为相对为相对稳定区间稳定区间 II 绝对稳定区间绝对稳定区间 若对于任意的若对于任意的h(a,b),有有|ri(h)|0,yk,不存在稳定区间了。不存在稳定区间了。II(r,h)=rp+1(a0rp+a1 rp-1+ap-1r+ap)h(b-1 rp+1+b0 rp+b1 rp-1+bp-1r+bp)=(r)h(r)为什么这样定义稳定性？为什么这样定义稳定性？看书吧？看书吧？微分方程数值解哈工大计算数学31例题6-66-7稳定区间讨论看书，学会微分方程数值解哈工大计算数学32
