ch6-3Runge-Kutta方法.ppt

1、Runge-Kutta方法方法微分方程数值解哈工大计算数学4Runge-Kutta方法方法(单步法单步法)(高精度、适合于变步长、多步法起步高精度、适合于变步长、多步法起步)4-1 RK方法方法 思想：对于思想：对于 y=f(x,y)y(a)=g0 xa,b (1)由由Lagrange中值公式有中值公式有 y(xk+1)-y(xk)=hy(k)y(xk+1)=y(xk)+hf(,y()(4-1)令令K=f(,y()，为了获得高精度的方法，为了获得高精度的方法，必须选择准确的必须选择准确的K。微分方程数值解哈工大计算数学 做法：用做法：用xk,xk+1 上若干个点的斜率上若干个点的斜率ki的线性

2、组合的线性组合逼近逼近K，于是有，于是有 k1=f(xk,yk)k2=f(xk+c2h,yk+a2,1hk1)k3=f(xk+c3h,yk+a3,1hk1+a3,2hk2)ks=f(xk+csh,yk+as,1hk1+as,s-1hks-1)把它们作线性组合把它们作线性组合 K=b1k1+b2k2+bsks (4-2)作为作为yk到到yk+1的的“增量增量”形成一类显式单步法，称为形成一类显式单步法，称为RK方法。方法。11skkkiiiyyhKybk(4-3)微分方程数值解哈工大计算数学4-2 参数选取原则参数选取原则(4-3)式中参数式中参数ci,ai,j,bi 的选取原则的选取原则 将将

3、(4-3)式右端项式右端项ki在在(xk,yk)处作二元函数的处作二元函数的Taylor展展开，按开，按h的幂次重新整理，使得的幂次重新整理，使得 与微分方程与微分方程(1)的解的解y=y(x)在在x=xk处的处的Taylor展开式展开式 比较比较h的的系数系数,希望希望有尽可能多的重合项，即要求有尽可能多的重合项，即要求 r1=fk ,r2=fk,r3=fk,r4=fk,这样得到这样得到(4-3)式叫做式叫做 S 级递推公式。级递推公式。如果如果 rj=fk(j-1),j=1,p而而rp+1fk(p)，则递推公式，则递推公式(4-3)具有具有 p 阶精度。阶精度。234112341112!3

4、!4!kkyyr hr hr hr h 2341111()()2!3!4!kkkkkky xy xf hf hf hf h(4-4)(4-5)微分方程数值解哈工大计算数学4-3 RK方法推演方法推演 通常使用较多的是通常使用较多的是 S=2,3,4 情况。情况。当当s=2时，二级时，二级RK算法算法 k1=f(xk,yk)k2=f(xk+c2h,yk+a2,1hk1)yk+1=yk+b1hk1+b2hk2 将将k2在在(xk,yk)作二元函数作二元函数Taylor展开展开222223222,11222,112,1121(,)(,)(,)0()2!Tkkkkffc hc hxx ykf x yf

5、x yc h a hkha hka hkffy xy 222222322,11222,1 12,111(2)0()2!xyxxxyyyfc hfa hk fc h fc h a k fa h k fh将将k1，k2 代入（代入（3-6），注意到），注意到k1=f，得，得112222,1322242222,12,1()()1(2)0()2!kkxyxxxyyyyyh bbfb h c faf fb h c hfc hafaffh右乘上右乘上hb2(3-6)微分方程数值解哈工大计算数学欲使欲使Tk+1的阶数尽可能高，应该选取的阶数尽可能高，应该选取b1,b2 xk,c2,a2,1使使h1,h2 的

6、的系数为零，即系数为零，即再将再将y(xk+1)在在xk作一元函数作一元函数Taylor展开展开2124()()2!1(2()0()3!kkxyxxxyyyxxyhy xy xhfff fffffffff fh2111122 22,1 2322242 22 2,1 22 2,111()(1)()()221 111 11()()()()0()6 236 26kkkxyxxxyyyyxyTy xyhb b fhbc fa b f fhbc fca b f fbaf ff ff fh 12222,1 210102102bbb ca bc2自由取值，自由取值，得到得到222,121212112bcac

7、bc不会等不会等于零于零微分方程数值解哈工大计算数学 由局部误差看，无论怎样选取由局部误差看，无论怎样选取c2都不会使都不会使h3系数为零，可见系数为零，可见S=2的的RK方法的精度阶数方法的精度阶数不会超过不会超过2阶。阶。常用的是：常用的是：取取c2=1/2,得到中点法得到中点法 取取c2=1/3,得到得到Heun法法 取取c2=1,得到改进的得到改进的Newton法法1(,(,)22kkkkkkhhyyh f xyf xy122(,)3(,(,)433kkkkkkkkhhhyyf xyf xyf xy1(,)(,(,)2kkkkkkkkhyyf xyf xh yhf xy微分方程数值解哈

8、工大计算数学4-4 四阶四阶RK方法方法 S=4时，可达到时，可达到4阶精度，有经典阶精度，有经典4阶阶KR方方法法121(,)1(,)22kkkkkhf x yhkhf xyk33431(,)22(,)kkkkhkhf xykkhf xh yk112341(22)6kkyykkkk书上例题书上例题6,76,7改进牛顿方法改进牛顿方法.4.4阶龙格库塔方法，并且给出阶龙格库塔方法，并且给出两种方法的比较。两种方法的比较。一句话，解决好一个实际课题不容易，实践性非常强，教一句话，解决好一个实际课题不容易，实践性非常强，教科书不会全部都写出来的！科书不会全部都写出来的！微分方程数值解哈工大计算数学

9、例题例题6-8,6-9 用改进的用改进的Euler法和经典四阶法和经典四阶R-K方法求解方法求解 PP298-2992,01(0)1,0.25xyyxyyh微分方程数值解哈工大计算数学4-5 局部误差局部误差 Tk+1=arhrr+1+o(hr+2)实际上Tk+1、ar都依赖于f(x,y)。例如。例如二级二级2阶方法阶方法211122 22,1 2322242 22 2,1 22 2,111()(1)()()221 111 11()()()()0()6 236 26kkkxyxxxyyyyxyTy xyhbb fhbc fa b f fhbc fca b f fbaf ffff fh 222,

10、121212112bcacbc222111()(2)()646xxxyyyyxyacff ff ffff f微分方程数值解哈工大计算数学4-64-6稳定性稳定性 单步法稳定多项式是一次的，只有一个根。单步法稳定多项式是一次的，只有一个根。只讨论绝对稳定区间只讨论绝对稳定区间。该方法用于检验方程该方法用于检验方程y=y，y(a)=g0 xa,b时，时，当当0时时 y(x)=g0 e(x-a),(x）当当yk=g0 ekh,(k+）时，发散；）时，发散；当当0时时 y(x)=g0 e(x-a)0,(x）可见当可见当yk 0,(k+）时，方法对于确定）时，方法对于确定的的h是绝对收敛的。例如是绝对收

11、敛的。例如6-9微分方程数值解哈工大计算数学二级二级2阶法用于检验方程阶法用于检验方程 y=y k1=f(xk,yk)=y k k2=f(xk+c2h,yk+a2,1hk1)=(y k+a2,1hy k)yk+1=yk+b1hk1+b2hk2=y k+b1hy k+b2h(y k+a2,1hy k)=(1+b1h+b2h+b2ha2,1h)y k =(1+h+hh/2)y k=(1+h+hh/2)k+1 y 0 可见可见|1+h+hh/2|1,稳定。稳定。222,121212112bcacbcS=2时时微分方程数值解哈工大计算数学 可见可见|1+h+hh/2|1,稳定。稳定。-11+h+hh/

12、21-2 h(1+h/2)0 h0 h-2 h0 1+h/20 h-2 不可能不可能 故故 绝对稳定区间是绝对稳定区间是 （2，0）微分方程数值解哈工大计算数学对于一般的 r 阶KR方法_1()()()1000()()kkkxaxh axahhkky xy ey ey e ee y x _23_(1)()2!3!khhhy x yk=y(x)时，时，y(xk+1)-yk+1=O(hr+1)所以所以S级级r价价KR方法应用于实验方程有方法应用于实验方程有 yk+1=Rs,r(h)yk其中其中 _2_,1()12!rsis rii rhhRhhv hr 称为RK方法的传递函数传递函数。微分方程数值

13、解哈工大计算数学其中系数vk依赖于特殊方法中的系数。若s=r则RK方法的绝对稳定区间的条件绝对稳定区间的条件为|Rs,r(h)|1 于是对于所有s=r的r阶方法，他们的绝对稳定区间是相同的，例如,r=s=3时 绝对稳定区间为（-2.51,0）例如,r=s=4时 绝对稳定区间为（-2.78,0）_23_3,3()12!3!hhRhh_234_3,3()12!3!4!hhhRhh4,4微分方程数值解哈工大计算数学6预估-校正法 隐格式线性多步法隐格式线性多步法 yk+1=a0yk+a1yk-1+apyk-p +h(b-1yk+1+b0yk+b1yk-1+bpyk-p)(2)I 起步时，yi,i=0

14、,1,p 要用单步法计要用单步法计算算启动启动；II 之后，每一步都要解一个非线性方程之后，每一步都要解一个非线性方程的根，初始值怎样给定才会有的根，初始值怎样给定才会有(1)()11!10(,)(,)pjjkkkik iik ik iiyb hf xyaybhf xy微分方程数值解哈工大计算数学 一般取一般取j=0，只迭代一次，例如，只迭代一次，例如(0)1(0)1!1(,)1(,)(,)2kkkkkkkkkkyyhf x yyyf x yf xyEuler格式格式梯形公式梯形公式问题：问题：I 预估预估效正法的特点：显格式与隐格式效正法的特点：显格式与隐格式交替使用交替使用II显格式与隐格式精度应该相同配对使用；显格式与隐格式精度应该相同配对使用；III启动阶段用精度不小于效正法的精度。启动阶段用精度不小于效正法的精度。微分方程数值解哈工大计算数学数学实验数学实验 求下列方程组的数值解求下列方程组的数值解 y1=0.04(1-y1)-(1-y2)y1+0.0001(1-y2)2 y2=-104y1+3000(1-y2)2 y1(0)=0,y2(0)=1,0 x100100