计算方法4常微分方程数值解法

1、常微分方程的数值解法Numerical Solutions to Ordinary Differential Equations 对象对象一阶常微分方程初值问题一阶常微分方程初值问题: 00y)x(y)y,x( fdxdy一阶常微分方程组初值问题一阶常微分方程组初值问题: 000n0n202101n21nnn2122n2111y)x(y,y)x(y,y)x(y)y,y,y,x(fy)y,y,y,x(fy)y,y,y,x(fy高阶常微分方程初值问题高阶常微分方程初值问题: )n(00)n(0000)n()n(y)x(y,y)x(y,y)x(y)y,y,y,x( fy(4.1)一阶常微分方程初值问

2、题一阶常微分方程初值问题:实际工程技术、生产、科研上会出现大量的微分方程问实际工程技术、生产、科研上会出现大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表题很难得到其解析解，有的甚至无法用解析表达式来表示，因此只能依赖于数值方法去获得微分方程的数值解。示，因此只能依赖于数值方法去获得微分方程的数值解。 00y)x(y)y,x( fdxdy常数常数称为称为有唯一解有唯一解则初值问题则初值问题使得使得条件：条件：上连续，且满足上连续，且满足在在若若定理定理LipschitzL,)1.4(|yy|L|)y,x( f)y,x( f|,y,y,LLipschitzb,a)y,x( f21

3、2121 ab x0 x1 x2 . xn-1 xn用数值方法，求得用数值方法，求得y(x)在每个节点在每个节点xk 的值的值y(xk ) 的近似的近似值，用值，用yk 表示，即表示，即yk y(xk)，这样，这样y0, y1,.,yn称为微分称为微分方程的数值解方程的数值解求求y(x)求求y0, y1,.,yn 微分方程的微分方程的数值解法数值解法：不求不求y=y(x)的精确表达式的精确表达式,而求离散点而求离散点x0,x1,xn处的处的函数值函数值设设(4.1) 的解的解y(x)的存在区间是的存在区间是a,b，初始点，初始点x0=a，取取a,b内的一系列节点内的一系列节点x0, x1,.,

4、xn。a= x0 x1 xn =b，一般采用等距步长，一般采用等距步长思路思路计算过程计算过程：方法方法：采用步进式和递推法：采用步进式和递推法将将a,bn等分等分, a= x0 x14)阶龙格库塔法。阶龙格库塔法。 m4时：计算量太大，精确度不一定提高，有时会降低。时：计算量太大，精确度不一定提高，有时会降低。 Gill公式公式 )hK)221(hK22y,hx( fK)hK)221(hK212y,2hx( fK)K2hy,2hx( fK)y,x( fK)KK)22(K)22(K(6hyy32nn421nn31nn2nn14321n1n节省存储单元节省存储单元控制舍入误差控制舍入误差 对于经

5、典的四阶对于经典的四阶Runge-Kutta法给出如下算法：法给出如下算法：算法算法4.2求解：求解： dy/dx=f(x,y) axb y (a)=y0 Step 1: 输入输入a,b,y0 及及NStep 2: (b-a)/N=h,a=x,y0=yStep 3: 输出输出 (x,y)Step 4: For I=1 T0 Nhf(x,y)=k1hf(x+h/2,y+ k1/2)= k2hf(x+h/2,y+k2/2)=k3hf(x+h,yk3)=k4y+(k1+2k2+2k3+k4)/6=yx+h=x输出输出(x,y)Step 5 : 停止停止例例4.3用四阶经典用四阶经典RungeKutt

6、a方法解初值问题：方法解初值问题： 102yyxydxdy2 . 0h(1)求求 ， 1y2 . 0, 1, 000hyx12,000001yxyyxfK9181818. 022222,2100101002KhyhxKhyKhyhxfK9086375. 022222,2200202003KhyhxKhyKhyhxfK8432399. 02,300303004hKyhxhKyhKyhxfK1832293. 126432101KKKKhyy 1832160. 14 . 11211xxy51103 . 1e(2)求求 ， 2y10.2,0.2xh1832293. 11y8451714. 02,111

7、111yxyyxfK211122,0.7946656hhKf xyK311222,0.7874945hhKf xyK4113,0.7440375Kf xh yhK21123421.34168036hyyKKKK3416408. 11222xxy52100 . 4e自适应龙格库塔自适应龙格库塔法法用户提出问题I ：问题：问题：如何判断：如何判断|y(xn)-yn| 精度值精度值y(xn)未知。未知。 ：如何取如何取h=?解解：如用：如用p阶龙格库塔法计算，局部截断误差为阶龙格库塔法计算，局部截断误差为O(hp+1) y=f(x,y) y(x0)=y0 要求误差要求误差10-8求数值解求数值解 x

8、n h/2 xn+1h如如 xn-xn+1 令令 yn=y(xn) yn+1(h) 则则 y(xn+1)-yn+1(h)chp+1步长折半步长折半xnxn+h/2xn+1分两步计算分两步计算y(xn+1)的近似值的近似值yn+1(h/2)。则则 y(xn+1)-yn+1(h/2)2c(h/2)p+1 h2n+1( )n+1(h)pn+1n+1y(x)-y1y(x)-y2h2h2nn+1+1( )(h)n+1pn+112y -yyy1(x)-定理：对于问题I若用P阶龙格库塔法计算y(xn+1)在步长折半前后的近似值分别为yn+1(h), yn+1(h/2)则有误差公式 ph2n+1h2n+112

9、( )n( )(h)n+111+y(x|)-yyy |-=|注：注：10 误差的事后估计法误差的事后估计法 20 停机准则：停机准则： (可保证可保证|y(xn+1)-yn+1(h/2)|） 将步长折半反复计算，直至将步长折半反复计算，直至为止，为止， 取取h为最后一次的步长为最后一次的步长， yn+1为最后一次计算的结果为最后一次计算的结果。 Else if (为止，最后为止，最后一次运算的前一次计算结果即为所需。一次运算的前一次计算结果即为所需。 4.5 收敛与稳定性收敛与稳定性 对于常微分方程初值问题对于常微分方程初值问题(4.1)求求y=y(x)-求求y(xn)- 求求yn xn=x0

10、+nh离散化离散化某种公式某种公式例：例：Euler法，改进的法，改进的Euler法，龙格库塔法都是单步法。法，龙格库塔法都是单步法。数值解法：数值解法： 00y = fx, yyx= y单步法：计算单步法：计算yn+1时只用到前一步的结果时只用到前一步的结果yn。显式单步法：显式单步法： yn+1=yn+h(xn,yn,h) (x,y,h)为增量函数，它依赖于f,仅是xn,yn, h的函数Def:若某数值方法对于任意固定的若某数值方法对于任意固定的xn=x0+nh,当当 h-0(n-)时，时， yn-y(xn)，则称该法是则称该法是收敛收敛的。的。0(0)(limnhnny xy即即 (xn

11、=x0+nh为固定值为固定值)改进改进Euler法的收敛性法的收敛性 判定改进判定改进Euler法的收敛性：法的收敛性： yn+1=yn+hf(xn,yn)+f(xn+h,yn+hf(xn,yn)/2 则则(x,y,h)=f(x,y)+f(x+h,y+hf(x,y)/2若：若：yo=y(x0) _|f(x,y)-f(x,y)| L|y-y|f(x,y)关于关于y满足李满足李-条件条件:12_ |(x,y,h)- (x,y,h)|f(x,y)-f(x,y)|+|f(x+h,)-f(x+y+hf(x,y)y+hh,f(x,y)|则：则： 12_L|y- y|+L|y+hf(x,y)- y-hf(x

12、,y)|_21L|y- y|+L|y- y|+L hhL(1+L)2|y- y|=|y- y|2限定限定hh0，则则_0_| (x,y,h)- (x,y,hhL(1+L)2)|y-y|即即 (x,y,h) 满足李普希兹条件满足李普希兹条件,由由Th改进的改进的Euler法收敛法收敛同样可验证，若同样可验证，若f(x,y)关于关于y满足李普希兹条件满足李普希兹条件,且且 y0=y(x0)时，时， Euler法，标准四阶龙格库法也收敛。法，标准四阶龙格库法也收敛。4.5.2 稳定性稳定性在讨论收敛性时，一般认可：数值方法本身计算在讨论收敛性时，一般认可：数值方法本身计算过程是准确的。实际上，并非如

13、此：过程是准确的。实际上，并非如此： 初始值初始值y0有误差有误差0y0-y(x0). 计算的每一步有舍入误差。计算的每一步有舍入误差。这些初始误差在计算过程的传播中，是逐步衰减这些初始误差在计算过程的传播中，是逐步衰减的，还是恶性增长，这就是稳定性问题。的，还是恶性增长，这就是稳定性问题。Def4.1若一种数值方法在节点若一种数值方法在节点xn处的数值解处的数值解yn的扰动的扰动n0，而在以后的各节点值而在以后的各节点值ym(mn)上有扰动上有扰动m.当当|m|n|，(m=n+1,n+2,)，则称该数值算法则称该数值算法是是稳定稳定的。的。分析某算法的数值稳定性，通常考察模型方程分析某算法的

14、数值稳定性，通常考察模型方程 y=y ， (0)Def4.1 ： 设在节点设在节点xn处用数值法得到的理想数值解为处用数值法得到的理想数值解为yn，而实际计算得到的近似值为，而实际计算得到的近似值为 ，称，称 为第为第n步数值解的扰动步数值解的扰动nynnnyy Euler法的稳定性法的稳定性Euler法：yn+1=yn+hf(xn,yn)考察模型方程考察模型方程 y=y，(0) 即即yn+1=(1+h)yn假设在节点值假设在节点值yn上有扰动上有扰动n，在，在yn+1上有扰动上有扰动n+1，且且n+1仅由仅由n引起引起(计算过程不再引进新的误差计算过程不再引进新的误差)Euler法稳定法稳定

15、Euler法的稳定的条件为法的稳定的条件为 0h-2/nn1n1n1ny)h1(y)h1(yy nnn1n)h1()yy)(h1( 1h1 n1n 隐式隐式Euler法稳定性法稳定性隐式隐式Euler法法：yn+1=yn+hf(xn+1,yn+1)对于模型方程对于模型方程 y=y( yn+1=yn/(1-h)yn+1的扰动的扰动0 上式均成立上式均成立，隐式隐式Euler法无条件稳定法无条件稳定h1h1yh1yyynnn1n1n1n 1h11Euler 法法稳稳定定隐隐式式梯形公式稳定性梯形公式稳定性梯形公式梯形公式 yn+1=yn+hf(xn,yn)+f(xn+1,yn+1)/2，设模型方程

16、为设模型方程为y=y(0)则则 yn+1=yn+hyn+yn+1/2 n+1nh1 +2y=yh1 -20时恒成立时恒成立 梯形公式恒稳定。梯形公式恒稳定。 yn+1的扰动的扰动n1n2h12h1 12h12h1 梯梯形形公公式式稳稳定定 00y = fx, yyx= yTaylor公式公式公式公式单步否单步否显式否显式否精度精度 n+1nnny=y +hf(x ,y )单步单步显式显式一阶一阶n+1nn+1n+1y=y +hf(x,y)单步单步隐式隐式一阶一阶n+1n-1nny=y+2hf(x ,y )二步二步显式显式二阶二阶数值积分法数值积分法hn+1nnnn+1n+12y=y + f(x

17、 ,y )+f(x,y)单步单步隐式隐式二阶二阶11n+1n12221nn2nn1y= y +k +kk = hf(x ,y )k = hf(x +h,y +k )单步单步显式显式二阶二阶RungeKutta方法方法局部截断误差局部截断误差整体截断误差整体截断误差收敛与稳定性收敛与稳定性4.6 多步法多步法某一步解和前若干步解的值有关某一步解和前若干步解的值有关1N,m,1mn)y,x( fbhyaymj1nmj1n1m0jm0jjmj1nj1n m步线性多步法的一般形式步线性多步法的一般形式其中其中a0,am-1, b0,bm为和为和h无关的常数无关的常数初值为初值为y0, y1, ,ym-

18、1若若bm0，方法是隐式的方法是隐式的(implicit)否则是显式的否则是显式的(explicit)(4.7)AdamsBashforth四阶方法四阶方法1N,4,3n)y,x( f9)y,x( f37)y,x( f59)y,x( f55(24hyy,y,y,y,y3n3n2n2n1n1nnnn1n3322110 AdamsMoulton四阶方法四阶方法1N,3,2n)y,x( f)y,x( f5)y,x( f19)y,x( f9(24hyy,y,y,y2n2n1n1nnn1n1nn1n22110 多步法的截断误差多步法的截断误差1N,m,1mn)x(y,x( fb)x(ya)x(y(h1)

19、h(7.4m0jmj1nmj1nj1m0jmj1nj1n1n ）式式的的截截断断误误差差为为定定义义（4.7 常微分方程组和高阶微分方程的数值解法常微分方程组和高阶微分方程的数值解法形如形如 0000z)x(z),z,y,x(gzy)x(y),z,y,x( fy欧拉方法的计算公式欧拉方法的计算公式 00nnnn1n00nnnn1nz)x(z),z,y,x(hgzzy)x(y),z,y,x(hfyy隐式欧拉方法的计算公式隐式欧拉方法的计算公式 001n1n1nn1n001n1n1nn1nz)x(z),z,y,x(hgzzy)x(y),z,y,x(hfyy梯形公式梯形公式 )z,y,x(g)z,y

20、,x(g(2hzz)z,y,x( f)z,y,x( f (2hyy1n1n1nnnnn1n1n1n1nnnnn1nAdams公式公式R-K方法方法对于高阶微分方程，总可以化成方程组的形式，对于高阶微分方程，总可以化成方程组的形式，例如例如 0000y)x(y,y)x(y)y,y,x( fy可以化成可以化成 0000z)x(z,y)x(y)z,y,x( fzzy例例6.0)0(y,4.0)0(y,1x0,xsiney2y2yx2 令令y(x)=z，则有则有 z2y2xsinezzyx2用用4阶阶RK方法，有方法，有6.0z,4.0y,0 x000 04.0)z2y2xsine(h)z,y,x(h

21、fk06.0hz)z,y,x(hfk000 x200022,1000011,10 70324764475.0)2kz,2ky,2hx(hfk062.0)2kz(h)2kz,2ky,2hx(hfk2,101,10022,22,102,101,10011,2 70315240923.0)2kz(2)2ky(2)05.0 xsin(e(hk80616283223.0)2kz(hk2,201,200)05.0 x(22,32,201,30 80217863729.0)kz(2)ky(2)1.0 xsin(e(hk40631524092.0)kz(hk2,301,300)1.0 x(22,42,301,

22、40 6/)kkkk(zz6/)kkkk(yy2,42,32,22,1011,41,31,21,101 4.8 常微分方程边值问题的数值解法常微分方程边值问题的数值解法形如形如bxa),y,y,x( fy 第一类边界条件第一类边界条件 )b(y,)a(y(4.8)第二类边界条件第二类边界条件 )b(y,)a(y第三类边界条件第三类边界条件 )b(y,)a(y定理定理 假定问题（假定问题（4.8）式中的函数）式中的函数f在集合在集合 y,y,bxa)y,y,x(D上连续，且偏导数上连续，且偏导数fy和和fy也在也在D上连续。若有上连续。若有M)y,y,x(fM) ii(0)y,y,x(f ,D)

23、y,y,x() i (yy 使使得得存存在在常常数数对对于于所所有有则边值问题有唯一解。则边值问题有唯一解。例例 边值问题边值问题0)2(y)1(y,2x1,0ysineyxy 有有ysine)y,y,x( fxy 因此边值问题有唯一解。因此边值问题有唯一解。1ycos)y,y,x(f,0 xe)y,y,x(fyxyy 定理定理 若线性边值问题若线性边值问题 )b(y,)a(y,bxa),x( ry)x(qy)x(py满足满足0)x(qb,a) ii()x( r)x(q),x(pb,a) i ( 上上在在连连续续和和上上在在则边值问题有唯一解。则边值问题有唯一解。打靶法打靶法（Shooting Method）)x(y)b(y)b(y)x(y)x(y2211 其中，其中，y1(x)是初值问题是初值问题的解的解0)a(y,)a(y,bxa),x( ry)x(qy)x(py 的解，的解，y2(x)是初值问题是初值问题1)a(y,0)a(y,bxa,y)x(qy)x(py 有限差分法有限差分法(Finite-Differe