数值计算方法第六章

1、第六章第六章 常微分方程初值问题的常微分方程初值问题的数值解法数值解法6.1 欧拉方法欧拉方法 6.2 计算公式的误差分析计算公式的误差分析 6.3 龙格库塔方法龙格库塔方法 6.4 向一阶方程组与高阶方程的推广向一阶方程组与高阶方程的推广问题的提出问题的提出的的数数值值解解求求初初值值问问题题)16(,)(),(000 bxxyxyyxfy数值求解方法数值求解方法的的近近似似值值数数值值逐逐个个求求解解出出节节点点上上的的函函生生成成节节点点作作等等距距分分割割对对区区间间)(,),(),()2(), 1 , 0(,)1(21000NixyxyxyNxbhniihxxbx 6.1 欧拉方法欧

2、拉方法6.1.1 欧拉公式与改进欧拉公式欧拉公式与改进欧拉公式算法：算法：Nxbhniihxxbxi000), 1 ,0(, 生生成成节节点点作作等等距距分分割割对对区区间间处处的的近近似似值值在在处处的的精精确确值值在在iiiixxyyxxyxy)()()( 111)(,()()()(,()(),(,11iiiiiixxiixxxxiidxxyxfxyxydxxyxfdxxyyxfyxx即即：积积分分上上对对在在)46()(,()()(11 iixxiidxxyxfxyxy化化为为积积分分方方程程问问题题从从而而把把微微分分方方程程问问题题转转选择不同的数值积分公式来求选择不同的数值积分公式

3、来求近似值就得到初值问题的各种数值解法近似值就得到初值问题的各种数值解法 1)(,(iixxdxxyxf1.欧拉公式欧拉公式)(,()()(1iiiixyxfhxyxy 得得，时时由由当当)46()(,()(,(1 iixxxyxfhdxxyxfiioxy),(yxfY 0 x1x1 ixbxn ixh)26()1, 1 , 0()(),(001 Nixyyyxfhyyiiii由由此此可可建建立立计计算算公公式式这称为这称为欧拉公式欧拉公式2.后退欧拉公式后退欧拉公式oxy),(yxfY 0 x1x1 ixbxn ixh得得，时时由由当当)46()(,()(,(111 iixxxyxfhdxx

4、yxfii)36()1, 1 , 0()(),(00111 Nixyyyxfhyyiiii由由此此可可建建立立另另一一公公式式)(,()()(111 iiiixyxfhxyxy这称为这称为后退欧拉公式后退欧拉公式后退欧拉公式是一个隐式公式，通常采用迭代法求解。后退欧拉公式是一个隐式公式，通常采用迭代法求解。比比较较并并与与精精确确解解的的数数值值解解xxexxyxyyxey )2(21)(,1 , 01)0(2例例6.1 以以 h=0.1为步长，用欧拉法求常微分方程初值问题为步长，用欧拉法求常微分方程初值问题)9 ,2 , 1 , 0(1)(01 iyyexhyyixiiii解解：由由欧欧拉拉

5、公公式式得得xiyiy(xi)y(xi)-yi01100.10.9000000.9093629.362x10-3.立表立表 (见表见表6-1（书（书121页）页）)6.1.2 梯形公式与改进欧拉公式梯形公式与改进欧拉公式3.梯形公式梯形公式)56()(),(),(2)(,()(,(2)()()(,()(,(2)(,(00111111111 xyyyxfyxfhyyxyxfxyxfhxyxyxyxfxyxfhdxxyxfiiiiiiiiiiiiiiiixxii从而导出梯形公式从而导出梯形公式得得由由-梯形公式梯形公式也是隐式单步法公式也是隐式单步法公式oxy),(yxfY 0 x1x1 ixbx

6、n ixh图图1 梯形公式梯形公式 用梯形公式计算时，通常取欧拉公式的解作为迭代初值进行迭用梯形公式计算时，通常取欧拉公式的解作为迭代初值进行迭代计算，即采用下式代计算，即采用下式),1 ,0(),(),(2),()(11)1(1)0(1 kyxfyxfhyyyxfhyykiiiiikiiiii )1(1)(1)2()3(1)2()2(1)2()1(1)2()0(1)1(0kkyyyyyyy(1)(2) )1(2)(2)2()3(2)2()2(2)2()1(2)2()0(2)1()0(1kkyyyyyyy11)(1)1(1)(1)(,|, iikikikiyxyyyyk的的近近似似值值值值作作

7、为为函函数数取取时时使使得得迭迭代代计计算算到到某某个个 ) 1, 1 , 0()2(),(),(2) 1 (),(,11111)1(1 Niyxfyxfhyyyxfhyyyyiiiiiiiiiiii则则得得到到显显式式计计算算公公式式代代替替取取直直接接简简单单地地4.改进欧拉公式改进欧拉公式这称为这称为改进欧拉公式改进欧拉公式NNyyyyyyyyy)2(3)2(3)1(2)2(2)1(1)2(1)1(0例例6.2 仍取步长仍取步长h = 0.1，采用改进欧拉法重新计算例，采用改进欧拉法重新计算例 6.1 的的常微分方程初值问题。常微分方程初值问题。)9 ,1 , 0(1)()(2)(011

8、111 iyyexyexhyyyexhyyixiixiiiixiiiiii(见表见表6-2（书（书125页）页）)这时改进欧拉公式为这时改进欧拉公式为解解xiyiy(xi)y(xi)-yi01100.10.9095240.9093621.363x10-4.立表立表6.2 计算公式的误差分析计算公式的误差分析 定义定义6.1 若若 yi+1 是是 yi=y(xi) 从计算得到的近似解，则称从计算得到的近似解，则称y(xi+1) yi+1为所用公式的为所用公式的局部截断误差局部截断误差，简称为，简称为截断误差截断误差。 定理定理6.1 若单步法 yi+1 = yi+h (xi , yi , h)

9、的局部截断误差为 O (h p+1) ，且增量函数 (x , y , h) 关于 y 满足李普希兹条件，即存在常数 L0，使对 成立不等式yy,| ),(),(|yyLhyxhyx则其整体截断误差 y(xi) yi=O(hp) 截断误差的估计截断误差的估计(基本假设：基本假设： yi = y( xi ) )设设 y(x) C 3 x0 , b , 则则 )76()()(2)()()()(321 hOxyhxyhxyhxyxyiiiii（1）对欧拉公式，有）对欧拉公式，有)86()()()(2)()()()(,()(232111 hOhOxyhyxyxyhxyxyxfhxyyiiiiiiiii故

10、因此，欧拉公式的局部截断误差为 O (h2)（2）对后退欧拉公式，有）对后退欧拉公式，有)96()()()(2)()(,()(23211111 hOhOxyhyxyxyxfhxyyiiiiiii故因此，后退欧拉公式的局部截断误差为 O (h2)（3）对梯形公式，注意到其公式可改写为）对梯形公式，注意到其公式可改写为 ),(),(21111 iiiiiiiyxfhyyxhfyy故由式（故由式（6-8）和（）和（6-9）得）得)106()()()(2)()(221),()(),()(21)(33232111111 hOhOxyhhOxyhyxfhyxyyxfhyxyyxyiiiiiiiiiiii因

11、此，梯形公式的局部截断误差为因此，梯形公式的局部截断误差为 O ( h3 )（4）对改进欧拉公式，有）对改进欧拉公式，有 ),(),(2)()()(1111iiiiiiiiiyxfyxfhxyyxyhxyy而由 ，故有),(),(),(yxfyyxfyyxfyyx 得)()(2)()()()()()()(,()()(,(),(32122)(,(11hOxyhxyhxyyhOxyhxyhOfyhfhxyxfxyhxyhxfyxfiiiiiixyxyxiiiiiiiii 所以所以与式（与式（6-7）比较得）比较得 y(xi+1) yi+1 = O ( h3 ) 因此，改进欧拉公式的局部截断误差为因

12、此，改进欧拉公式的局部截断误差为 O ( h3 ) 定义定义6.2 若一种求解常微分方程初值问题的数值计算方法若一种求解常微分方程初值问题的数值计算方法的局部截断误差为的局部截断误差为 O ( hp+1 ) ，则称该方法为，则称该方法为 p阶精度阶精度，或称该，或称该方法为方法为 p阶方法阶方法。 由此定义知，欧拉方法与后退欧拉方法为一阶精度，梯由此定义知，欧拉方法与后退欧拉方法为一阶精度，梯形法与改进欧拉方法为二阶精度。形法与改进欧拉方法为二阶精度。6.3 龙格龙格-库塔方法库塔方法由中值定理，有由中值定理，有),(,)(,()()()()(111 iiiiiixxyfhyxxxyxy 因此

13、，以上介绍的各种单步法本质上都是对平均斜因此，以上介绍的各种单步法本质上都是对平均斜率率 f( , y( ) 进行近似，龙格进行近似，龙格-库塔据之提出了适当选取若库塔据之提出了适当选取若干点上的斜率值作近似以构造高精度计算公式的方法，其干点上的斜率值作近似以构造高精度计算公式的方法，其基本思想是基于泰勒展式的待定系数法。基本思想是基于泰勒展式的待定系数法。6.3.1 二阶二阶Rung-KuttaRung-Kutta公式公式问题：问题：建立二阶精度的计算格式形为建立二阶精度的计算格式形为待待定定系系数数为为这这里里babhKyahxfKyxfKKKhyyiiiiii,)126(),(),()(

14、2112122111 在在 y(xi) = yi 的假设下，有的假设下，有)()()()(,()(,()()(,(2)(,(2)(,(1121hOfybfahxyhOfbhKfahxyxfbhKxyahxfKxyxyxfKiiiixyxyxixyxyxiiiiiii )()()()(3)(,(222211hOfybfahxyhxyyiixyxyxiii 故故解解)()()()(3)(,(222211hOfybfahxyhxyyiixyxyxiii 变形欧拉公式变形欧拉公式根据格式为二阶精度，即根据格式为二阶精度，即 y(xi+1) yi+1 = O(h3) 比较两式系数得比较两式系数得)136

15、(2/12/112221 ba而而)(2121)()()()()2/()()()(3)(,(2321hOfyfhxyhxyhOxyhxyhxyxyiixyxyxiiiiii 系数满足系数满足(6-13)的形为的形为(6-12)计算格式统称为二阶计算格式统称为二阶R-K公式。公式。当令当令1=1/2时，解得时，解得 2=1/2 ，a=b=1，即为改进欧拉公式。，即为改进欧拉公式。若令若令 1=0，解得，解得 2=1，a=b=1/2，则得另一计算公式，则得另一计算公式)146()2/, 2/(),(12121 hKyhxfKyxfKhKyyiiiiii6.3.2 四阶四阶 R-K R-K 公式公式

16、每一步需计算的 f 值的个数1234567n8精度阶1234456n-2 1965年，年，Butcher研究发现显式研究发现显式R-K公式的精度与需要组公式的精度与需要组合的斜率值的个数具有如下关系合的斜率值的个数具有如下关系 可见，超过四阶精度的可见，超过四阶精度的R-K公式效率并不高，实际计算通公式效率并不高，实际计算通常选用如下四阶格式常选用如下四阶格式)156(),()2/, 2/()2/, 2/(),()22(6342312143211hKyhxfKhKyhxfKhKyhxfKyxfKKKKKhyyiiiiiiiiii经典经典R-KR-K公式公式这时经典这时经典R-K公式为公式为 例

17、例6.3 取步长取步长h = 0.2，采用经典，采用经典R-K法计算例法计算例 6.1 的常微的常微分方程初值问题。分方程初值问题。 )()()2()2()2()2()22(63)(42)2/(31)2/(2143211hKyehxKKhyehxKKhyehxKyexKKKKKhyyihxiihxiihxiixiiiiiii 取取 h=0.2 计算得到表计算得到表6-4(书书133页）。页）。 与例与例6.1和例和例6.2比较可见，用经典比较可见，用经典R-K法计算得到的解比法计算得到的解比用欧拉法和用欧拉法和改进欧拉法改进欧拉法所得到的解精确得多。所得到的解精确得多。解解6.3.3 步长的自动选择步长的自动选择 对于对于 p 阶精度的计算格式，当取步长为阶精度的计算格式，当取步长为 h 时，记时，记 为为从从 y(xi) 计算得到的计算得到的 y (xi+1) (xi+1= xi+h) 的近似解，则有的近似解，则有1hiy )(121)(21)()(,)(22)()()(12/12/11112/1112/11111hihiphiiphiihiiphiiphiiyyyxyyxyyxyCChChCyxyCChyxy 所以所以故有故有较小时较小时为常数为常数为常数为常数 为便于进行事后误差估计，实际计算时通常采用