数值分析ch4
第四章第四章 常微分方程的数值解法常微分方程的数值解法 Numerical Solutions to Ordinary Differential Equations 概述概述 () () :P roblem 00 dy = fx, y dx yx= y 121212 L, y ,y , |f(x,s.y )- f(x,y )|L|y|t.- y Ø一阶常微分方程初值问题一阶常微分方程初值问题: Øf(x,y)在在a,b上连续上连续,且满足且满足Lipschitz条件条件: Ø实际工程技术、生产、科研上会出现大量的微分方程问题实际工程技术、生产、科研上会出现大量的微分方程问题 很难得到其解析解，有的甚至无法用解析表达式来表示，很难得到其解析解，有的甚至无法用解析表达式来表示， 因此只能依赖于数值方法去获得微分方程的数值解。因此只能依赖于数值方法去获得微分方程的数值解。 求函数求函数y=y(x)满足：满足： y(x) =f(x, y(x) 则初值问题则初值问题Problem I有唯一解有唯一解 Ø微分方程的微分方程的数值解法数值解法： p不求不求y=y(x)的精确表达式的精确表达式,而求离散点而求离散点x0,x1,xn处的函数值处的函数值 p设设Problem I的解的解y(x)的存在区间是的存在区间是a,b，初始点，初始点x0=a，取，取 a,b内的一系列节点内的一系列节点x0, x1,.,xn。a= x04)阶阶龙格库塔龙格库塔法。法。 m4时时：计算计算量量太太大，精确大，精确度度不一不一定提定提高高，有，有时时会会降降低低。 龙格库塔龙格库塔法法 对于经典的四阶对于经典的四阶Runge-Kutta法给出如下算法：法给出如下算法： 求解求解： dy/dx=f(x,y) axb y (a)=y0 ØStep1: 输入输入a,b,y0及及N ØStep 2: (b-a)/N=h,a=x,y0=y ØStep3: 输输出出 (x,y) ØStep4: For I=1 T0 N phf(x,y)=k1 phf(x+h/2,y+ k1/2)= k2 phf(x+h/2,y+k2/2)=k3 phf(x+h,yk3)=k4 py+(k1+2k2+2k3+k4)/6=y px+h=x p输输出出(x,y) ØEND 龙格库塔龙格库塔法法 例例：用四阶经典：用四阶经典RungeKutta方法解初值问题：方法解初值问题： ( ) = = 10 2 y y x y dx dy 2 . 0=h (1)求，求， 1 y2 . 0, 1, 0 00 =hyx () 0 1000 0 2 ,()0.2 x Khfxyh y y = 2001 0 01 01 1 , 22 2 12 ()0.18363636 1 2 2 h KhfxyK h x h yK yK =+ + =+= + 3002 0 02 02 1 , 22 2 12 ()0.1817275 1 2 2 h KhfxyK h x h yK yK =+ + =+= + () () 4003 0 03 03 , 2 0.16864798 Khfxh yK xh h yK yK =+ + =+= + () 101234 1 22 6 1.1832293 yyKKKK=+ = ( )1832160. 14 . 112 11 =+=xxy 5 111 ()1.3 10ey xy =× 龙格库塔龙格库塔法法 (2)求求， 2 y 1 0.2,0.2xh=1832293. 1 1 =y () 1 1111 1 2 ,()0.16903428 x Khf x yh y y = () 1 211122 ,0.15893312 h Khf xyK=+= () 1 311222 ,0.1574989 h Khf xyK=+=() 4113 ,0.1488075Khf xh yK=+= () 211234 1 221.3416803 6 yyKKKK=+= ()3416408. 112 22 =+=xxy 5 2 100 . 4 ×=e 变变步长步长龙格库塔龙格库塔法法 Ø用户提出问题用户提出问题I ： 问题问题：如如何何判判断断|y(xn)-yn| 精确值精确值y(xn)未知未知。 y=f(x,y) y(x0)=y0 要要求求误差误差10-8求数值解求数值解 xn h/2 xn+1 hØ如如 xn-àxn+1 令令 yn=y(xn) yn+1(h) 则则 y(xn+1)-yn+1(h)chp+1 步长折步长折半半xnàxn+h/2àxn+1分分两两步计算步计算y(xn+1)的的近似近似值值yn+1(h/2)。 则则y(xn+1)-yn+1(h/2)2c(h/2)p+1 h 2 n+1 ( ) n+1 (h)p n+1n+1 y(x)- y 1 y(x)- y2 h 2 h 2 nn+1+1 ( ) (h) ) n+1 p n+1 1 2 y - yyy 1 (x)- ：如如何何取取h=? 解解：如如用用p阶阶龙格库塔龙格库塔法法计算计算，局部截断误差为局部截断误差为O(hp+1) 变变步长步长龙格库塔龙格库塔法法 定理定理：对于问题：对于问题I若用若用P阶龙格库塔法计算阶龙格库塔法计算y(xn+1)在步长折半在步长折半 前后的近似值分别为前后的近似值分别为yn+1(h), yn+1(h/2)则有误差公式则有误差公式 p h 2 n+1 h 2 n+1 1 2 ( ) n ( ) (h) n+ 1 11+ y(x|)- yyy |-=| 注注：10误差误差的的事事后后估估计计法法 20停机准停机准则：则：(可可保保证证|y(xn+1)-yn+1(h/2)|）将步长折将步长折半反复半反复计算计算，直直至至为为止止， 取取h为为最最后后一一次次的的步长步长， yn+1为为最最后后一一次次计算计算的的结果结果。 Else if (为为止止， 最最后后一一次运次运算算的的前前一一次次计算计算结果结果即为所需即为所需。 HW: p.116 #3（用用经典经典四阶四阶Runge-Kuttta法法） #8（选选作）作）
