chap7常微分方程的数值解法ppt课件

1、第第7章章 常微分方程组的数常微分方程组的数值解法值解法刘东毅刘东毅天津大学理学院数学系天津大学理学院数学系第第7章章 常微分方程组的数值解法常微分方程组的数值解法主要目标：主要目标： 掌握常微分方程初掌握常微分方程初值问题数值解法值问题数值解法的基本理论的基本理论掌握计算机上的常掌握计算机上的常用算法用算法主要内容：主要内容：初值问题计算格式的建初值问题计算格式的建立立Runge-Kutta 方法方法一阶常微分方程组与高一阶常微分方程组与高阶方程的数值解法阶方程的数值解法第第7章章 常微分方程组的数值解常微分方程组的数值解法法 在科学研究和工程实践中会遇到很多微分方程，虽然从理论上可以证明其

2、解的存在性，但其解的解析表达式往往是很难求解的，或者即使可以写出来，但也难于计算，此时，只能借助数值解来解决问题. 常微分方程组定解问题是自然科学和工程技术领域中常见的数学模型. 本章介绍求解此类问题的基本理论和数值解法。 yR定义定义7.0.1 若存在常数若存在常数 L 0, 使得对一切的使得对一切的 xa , b 及及 y , , 均有均有 |( , )( , )|,f x yf x yL yy则称则称 f (x, y) 在在 D 上关于上关于 y 满足满足 Lipschitz 条件条件, 其中其中 L 称为称为 Lipschitz 常数常数. 我们首先考虑一阶常微分方程初值问题0( ,

3、) ,( ).yf x yaxby ay 其中其中 f (x , y) 是区域是区域( , )|,Dx yaxbyR上的实值函数上的实值函数.(7.0.1) 我们首先给出常微分方程初值问题解的存在惟一性定理。定理定理 假设假设 f (x, y)C(D),且关于且关于 y 满足满足 Lipschitz 条件条件, 则一阶常微分方程初值问题则一阶常微分方程初值问题(7.0.1) 存在唯一解存在唯一解. 下面在此前提下下面在此前提下, 我们讨论上述初值问我们讨论上述初值问题题 (7.0.1) 的数值解法。的数值解法。0( , ) ,( ).yf x yaxby ay 012,Naxxxxb然后在节点

4、上建立逼近于原初值问题的计算格然后在节点上建立逼近于原初值问题的计算格式式 (或差分格式或差分格式), 由此计算出原问题的解由此计算出原问题的解 y( x ) 在节点在节点 x1 , x2 , . . . , xN 处的近似值：处的近似值： y1 , y2 , . . ., yN， 称它们为常微分方程初值问题的称它们为常微分方程初值问题的数值解数值解. 相邻两个节点的距离相邻两个节点的距离 hn= xn+1 - xn 称称为步长为步长, 通常取定步长通常取定步长h 0, 即节点即节点 xn = x0 + nh , n = 0,1, , N. 其基本思想是在区间 a , b 上引入一系列节点7.

5、1 初值问题计算格式的建立初值问题计算格式的建立1. 数值微分方法数值微分方法在等距节点下讨论问题在等距节点下讨论问题. 利用两点数值微分公式利用两点数值微分公式7.1.1 计算格式的建立计算格式的建立将上式代入初值问题将上式代入初值问题 (7.0.1),)()(0yayx,yfy1()()(, ()() ,2nnnnny xy xhf xy xyh有有11()()()(),()2nnnnnnny xy xhy xyxxh(7.1.1) 略去余项略去余项, 并以数值解并以数值解 yn , yn+1 替代替代 y (xn) 及及 y (xn+1), 则得差分方程则得差分方程上式称为上式称为 Eu

6、ler 公式公式. 利用此式可由初值利用此式可由初值 y0 出发按出发按“步进式步进式” 方法方法, 逐步求得数值逐步求得数值解解y1 , y2 , . . . , yN .)2 . 1 . 7(),(1nnnnyxhfyy由于计算由于计算yn+1时时, 只只用到它前一步的结用到它前一步的结果果yn , 这类公式称为这类公式称为单步法单步法. 又因为其关又因为其关于于yn+1是显式形式是显式形式, 故称该故称该Euler公式为公式为显格式显格式.如果利用下列数值微分公式如果利用下列数值微分公式111()()()(),()2nnnnnnny xy xhy xyxxh由由 类似的可导出类似的可导出

7、)()(xx,yfxy),(111nnnnyxhfyy上述公式称为后退的上述公式称为后退的 Euler 公式公式, 此公式为此公式为单步法公式单步法公式. 又因为它关于又因为它关于 yn+1 成隐式形成隐式形式式, 所以该公式为隐式公式，简称隐格式所以该公式为隐式公式，简称隐格式.2111()()()(),()26nnnnnnny xy xhy xyxxh类似地，可导出类似地，可导出上述公式称为上述公式称为Euler两步法公式两步法公式. 这因为，当这因为，当计算计算 yn+1 时时, 要用到要用到 yn -1 与与 yn . 显然它也是显然它也是显格式显格式.如果利用下列三点数值微分公式如果

8、利用下列三点数值微分公式112).(,nnnnyyhf xy(7.1.3) 设设 y (x)C2a , b, 由由 Taylor 公式有公式有由于由于 故上式即为故上式即为 ( )( , ( ) ,y xf x y x略去余项略去余项, 并以并以 yn , yn+1 替代替代 y (xn) 及及 y (xn+1), 得到的差分方程正是得到的差分方程正是Euler 公式公式.211()()()()().2nnnnnnnhy xy xhy xyxx(7.1.4) )5 . 1 . 7()(! 2)(,()()(21nnnnnyhxyxhfxyxy 2. Taylor 展开法展开法3. 数值积分方法

9、数值积分方法对对 ，在区间，在区间 xn , xn+1 上积分，得上积分，得( )( , )y xf x y11( )( , ( ),nnnnxxxxy x dxf x y x dx则有则有11()()( , ( ).nnxnnxy xy xf x y x dx对上式中的积分采用不同的数值积分公式可对上式中的积分采用不同的数值积分公式可得到不同的差分方程得到不同的差分方程. 例如例如, 对上式的积分采对上式的积分采用左矩形公式用左矩形公式, 可得到可得到 Euler 公式公式. 11131()() (, ()(, ()2(, ()().12nnnnnnnnnnnhy xy xf xy xf x

10、y xhfyxx若对此式的积分采用梯形公式若对此式的积分采用梯形公式, 11()()( , ( ).nnxnnxy xy xf x y x dx则有则有若略去余项若略去余项, 以以 yn , yn+1 替代替代 y (xn) 及及 y (xn+1), 得到的差分方程得到的差分方程111 (,)(,2.)nnnnnnhyyf xyf xy上式称为梯形公式上式称为梯形公式. 由于它关于由于它关于 yn+1 成隐式成隐式形式形式, 故其为隐格式故其为隐格式. 隐格式求解比较困难隐格式求解比较困难, 当当 yn 已知时已知时, 要求要求 yn+1 ，需解关于，需解关于 yn+1 的非线性方程的非线性方

11、程. 在实际在实际应用时应用时, 上式常与上式常与 Euler 公式联合使用公式联合使用, 构成构成如下计算格式如下计算格式:111 (,)(,2.)nnnnnnhyyf xyf xy(7.1.6)(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk隐式梯形公式的迭代格式隐式梯形公式的迭代格式(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk(7.1.7)由上式可以得到一个序列由上式可以得到一个序列: , k = 0,1, , 关于此序列的收敛性

12、关于此序列的收敛性, 有如下的定理有如下的定理.( )1kny定理定理 7.1.1 设设 f (x , y) 在区域在区域 D 上关于上关于 y 满满足足 Lipschitz 条件条件, 即即|( , )( , )|()|.f x yf x yLyy其中其中 L 为为 Lipschitz 常数常数, 当步长当步长 时时, 对任意的初值对任意的初值 按格式按格式7.1.7）2hL(0)1ny(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk生成的序列生成的序列 收敛于梯形公式收敛于梯形公式7.1.6）( )1kny111

13、 (,)(,)2nnnnnnhyyf xyf xy1ny的解的解 .为了减少计算量, 可采用预测-校正格式. 方法是先用Euler公式求得一个初始近似值 称为预测值, 再把 带入梯形公式右端计算一次求得 yn+1 称之为校正值, 即1ny1ny预测预测:校正校正:1(,) ,nnnnyyhf xy111 (,)(,).2nnnnnnhyyf xyf xy上式称为预测上式称为预测 - 校正公式或改进的校正公式或改进的 Euler 公式公式. 上式也可写成如下形式上式也可写成如下形式:11 (,)(,(,)2.nnnnnnnnhyyf xyf xyhf xy例例7.1.1：利用：利用Euler公式

14、与改进的公式与改进的Euler公公式求解初值问题步长式求解初值问题步长h = 0.1）. 1)0(, 10,2yxyxyy解：由步长解：由步长h=0.1，知节点，知节点 设数值解为设数值解为 利用利用Euler公公式得式得,1 . 00nnhxxn.10, 2 , 1 , 0n,nynnnnnnnnnnnyxyyxyyyxhfyy2 . 01 . 121 . 0),(1计算结果见下表计算结果见下表(见书见书P227表表7.1) 此初值问题的解析解为 , 从上表可以看出, 数值解 yn与解析解 y(xn) 比较, yn精度较差. xy21解此问题的改进的解此问题的改进的Euler公式为公式为 )

15、,(1nnnnyxhfyy),(),(2111nnnnnnyxfyxfhyy.1 . 011. 0105. 12221 . 0,2 . 01 . 121 . 01111111nnnnnnnnnnnnnnnnnnnnnyxyxyyxyyxyyyyxyyxyyy同同Euler公式比较公式比较, 改进的改进的Euler法显然精度提高了法显然精度提高了.由于误差大小是评价计算格式优劣的重要依据由于误差大小是评价计算格式优劣的重要依据, 故需要给出有关误差的概念故需要给出有关误差的概念. 计算结计算结果见下果见下表表(见书见书P228表表7.2)7.1.2 截断误差与方法的精度截断误差与方法的精度定义定

16、义 7.1.1 称误差称误差 en+1 = y ( xn+1 ) - yn+1为数值方法在点为数值方法在点 xn+1 的截断误差的截断误差, 又称整体截又称整体截断误差断误差. 设 yk= y ( xk ) (k = 0,1,. . . , n)，那么 为数值方法在点 xn+1 的局部截断误差.111()nnny xy整体截断误差整体截断误差 en+1 是在没有引进舍入误差的是在没有引进舍入误差的情况下情况下, 纯粹因为不准确的计算格式造成的纯粹因为不准确的计算格式造成的, 故又称为方法误差故又称为方法误差.它不仅与它不仅与 x = xn+1 这一这一步的计算有关步的计算有关, 而且和而且和

17、xn , xn-1 ,. . . , x1 这几这几步的计算都有关系步的计算都有关系. 局部截断误差是假设局部截断误差是假设 xn 之前各数值解没有误之前各数值解没有误差差, 仅由仅由 xn 到到 xn+1 这一步计算由计算格式这一步计算由计算格式引起的误差引起的误差.如如Euler公式公式1(,)nnnnyyhf xy在点在点 xn+1 的整体截断误差的整体截断误差 en+1 = y (xn+1)- yn+122()(, ()()(,)2(, ()(,)()2nnnnnnnnnnnnnhy xhf xy xyyhf xyheh f xy xf xyy局部截断误差局部截断误差2111()()2

18、nnnnhy xyy定义定义7.1.2 若某数值方法的局部截断误差为若某数值方法的局部截断误差为 则称该方法具有则称该方法具有 P 阶精度阶精度, 或称其为或称其为 P 阶方法阶方法.11(),pnO h可以证明可以证明:Euler 方法的局部截断误差方法的局部截断误差 其具有其具有一阶精度一阶精度. 梯形方法的局部截断误差梯形方法的局部截断误差 其具有二阶精度其具有二阶精度.改进的改进的 Euler 方法的局部截断误方法的局部截断误差差 具有二阶精度具有二阶精度.21(),nO h31(),nO h31(),nO h7.2 Runge-Kutta 方法方法 继续讨论前面的 Taylor 展开

19、法。设设 y (x)C2a , b, 由由 Taylor 公式有公式有211()()()()().2nnnnnnnhy xy xhy xyxx 由 故上式即为 ( )( , ( ),y xf x y x)(! 2)(,()()(21nnnnnyhxyxhfxyxy 略去余项略去余项, 并以并以 yn , yn+1 替代替代 y (xn) 及及 y (xn+1), 得到得到Euler 公式公式.进一步假设设进一步假设设 y (x)Cp+1a , b,由由 Taylor 公式有公式有) 1 . 2 . 7(,)(!)(! 2)()()()(21nnppnnnnRxyphxyhxyhxyxy )2

20、. 2 . 7.(),()()!1(11)1(1nnnpnppnxxhOyphR其中其中 由 故式(7.2.1)即为 ( )( , ( ),y xf x y x7.2 Runge-Kutta 方法略去余项略去余项, 并以数值解并以数值解 yn , yn+1 替代替代 (7.2.3) 中中的的 解析解解析解y (xn) 及及 y (xn+1), 可得到一个差分方可得到一个差分方程，即程，即nnppnnnnRxyphxyhxyhxyxy )(!)(! 2)()()()(21) 3 . 2 . 7(,)(,(!)(,(! 2)(,()()()1(21nnnppnnnnnnRxyxfphxyxfhxy

21、xhfxyxy., )()(,()!1(11)(1nnnpnnppnxxhOyfphR其中余项可写成其中余项可写成)(,(),()()(xyxfdxdyxfkkk注：这里注：这里在在(7.2.3)中略去余项中略去余项,用用yn , yn+1 替代替代 y (xn) 及及 y (xn+1)4 . 2 . 7(.! 2! 2)1(1)1(21pnpnnnpnpnnnnfphfhfhyfphfhfhyy, ),(nnnyxff . ) 1, 2 , 1(),()()(piyxffnniin其中其中称称 (7.2.4) 式为求解常微分方程初值问题数值解式为求解常微分方程初值问题数值解Taylor的格式

22、的格式 .21(1)()()(,()(,()2!(,(),(7.2.3)!nnnnnnppnnnhy xy xhf xy xfxy xhfxy xRp由于局部截断误差由于局部截断误差 ),(O)(1111pnnnhyxy可知它是一个可知它是一个 p 阶方法。当阶方法。当p=1时，上式正是时，上式正是Euler 公式。但当公式。但当 p 2 时，需要计算时，需要计算f (x, y(x) ) 的高阶导数，特别是对于复杂函数的高阶导数，特别是对于复杂函数 f (x, y(x) 的的求导，这无疑是大大增加计算量，这是它最大求导，这无疑是大大增加计算量，这是它最大的缺点。因此高阶的的缺点。因此高阶的Ta

23、ylor方法是不实用的。方法是不实用的。 德国数学家德国数学家C.Runge 及及M.W.Kutta提出了提出了一种改进策略，得到了至今还被作为高精度的一种改进策略，得到了至今还被作为高精度的单步法广泛使用龙格单步法广泛使用龙格-库塔法库塔法 ( Runge- Kutta method )。 7.2.1 Runge-Kutta方法的基本思想方法的基本思想.! 2)1(11pnpnnnnfphfhfhyyRunge-Kutta方法是利用 f 在某些点处函数值的线性组合替代7.2.4步长 h 后面括号中的因子来构造差分方程, 从而避免了高阶导数的计算, 这就是 Runge-Kutta 方法的基本思

24、想. )1(1! 2pnpnnfphfhf用用f 在某些点处函数值的线性组合替代这一部分在某些点处函数值的线性组合替代这一部分其一般形式为其一般形式为:11111,(,),(,)(2,3, )rnniiinniininijjjyyhkkf xykf xh yhkir其中其中 r 是上式中调用是上式中调用 f 的个数的个数, r 称为级数，称为级数， 为待定参数为待定参数, 适当确定这些参数适当确定这些参数, 可使可使之具有尽可能高的精度之具有尽可能高的精度. 如局部截断误差满足如局部截断误差满足,iii j . )(O)(1111rnnnhyxy 7.2.2 二阶二阶 Runge-Kutta

25、方法方法考虑考虑 r = 2 的情况的情况, 此时有此时有11 122121(),(,),(,).nnnnnnyyhkkkf xykf xh yh k利用二元函数的利用二元函数的 一阶一阶Taylor 公式，即全微分公式公式，即全微分公式 希望适当选择参数希望适当选择参数 使上式的局部截断误差为使上式的局部截断误差为12, , 3111()()nnny xyO h即为二阶方法即为二阶方法.,)()()(,()(,(),(),(22nnnnnynnnxnnyyxxOyyyxfxxyxfyxfyxf下面将下面将yn+1与与y(xn+1)作比作比较较2222( ,)(,)(,)()(,)()()()

26、() ()() ()()()nnxnnnynnnnndefnxnnynnnnf x yf xyfxyxxfxyyyOxxyyffxxfyyOxxyy. )O()()()(212hkfhfhfknynxn从而有, )(22111kkhyynn将上式代入),(12khyhxfknn再由得到. )()()()()O()()()()(322212121122111hOfffhfhyhkfhfhfhkhykkhyynynnxnnnynxnnnn在下面要将yn+1与y(xn+1)作比较，使它们的局部截断误差满足. )O()()()(212hkfhfhfknynxn, )(O)(3111hyxynnn为此考

27、虑y(xn+1)。再根据 y(xn+1) 在点 xn的一元 3 阶 Taylor 展开式. )()()(! 2)(),(! 2)()(! 2)()()(323),(2321hOfffhhfyhOyxfhhfyhOxyhxyhxyxynynnxnnyxnnnnnnnn )()()()(322211hOfffhfhyynynnxnnn由刚才已得到的, )(O)(3111hyxynnn让它们满足即由122211212左式含有四个未知元三个左式含有四个未知元三个方程方程, 因此解不唯一因此解不唯一. 参数参数满足左式的一族公式统称满足左式的一族公式统称二阶二阶 Runge-Kutta 公式公式. 可得

28、参数应满足下列方程组:. )()()(! 2)(321hOfffhhfyxynynnxnnn)()()()(322211hOfffhfhyynynnxnnn 取上式称为中点公式上式称为中点公式 .1210,1,212,nnyyhk1(,),nnkf xy21(,).22nnhhkf xyk 取上式称为上式称为 Heun 公式公式 .112(3),4nnhyykk1(,),nnkf xy2122(,).33nnkf xh yhk12132,443 可见, 二阶 Runge-Kutta 公式, 每计算一步需要 两次调用 f 的函数值. 取121,1,2112121(),2(,),(,).nnnnn

29、nhyykkkf xykf xh yhk得这正是改进的这正是改进的 Euler 公式公式. 7.2.3 四阶Runge-Kutta方法当当 r = 4 时时, 类似地可导出四阶类似地可导出四阶 Runge-Kutta 公公式式, 这种公式也有一族这种公式也有一族, 其中常用地有其中常用地有:规范规范 (经典经典) 的的 Runge-Kutta 方法方法11234(22),6nnhyykkkk1(,),nnkf xy21(,),22nnhhkf xyk32(,),22nnhhkf xyk43(,).nnkf xh yk11234(22)(22),6nnhyykkkk1(,),nnkf xy21(

30、,),22nnhhkf xyk3122122(,),222nnhkf xyhkhk423222(,).22nnkf xh yhkhk Gill 公式公式 Gill 公式是标准的 Runge-Kutta 公式的改进形式, 这种算法可节省存储单元, 并能控制舍入误差的增长.四阶四阶 Runge-Kutta 公式公式, 每一步计算需四次调用每一步计算需四次调用 f 的函数值的函数值, 计算量较大计算量较大, 但其局部截断误差可达但其局部截断误差可达 O(h5), 精度较高精度较高.例例7.2.1 用标准的用标准的 四阶四阶Rung-Kutta 法解初值问题法解初值问题,取取步长步长h=0.2.2,

31、(0)1,(01).xyyyxy 解解: 解此问题的计算公式为解此问题的计算公式为112340.2(22),6nnyykkkk12,nnnxkyy21122,22()nnnhxhkykhyk32222,22()nnnhxhkykhyk4332,2()nnnxhhkykyhkxn0.20.40.60.81.0ynyn- y (xn)1.183 21.341 71.483 3 1.612 5 1.732 10.000 00.000 00.000 0 0.000 10.000 1计算结果如下计算结果如下:显然在计算量大致相同的情显然在计算量大致相同的情况下况下, 标准的标准的 Runge-Kutta

32、方方法比改进的法比改进的 Euler 方法精确方法精确度更高度更高. （参见（参见p227和和p228的表）的表）7.5一阶常微分方程组与高阶方程一阶常微分方程组与高阶方程初值问题的数值解法初值问题的数值解法 7.5.1 一阶常微分方程组初值问题一阶常微分方程组初值问题 )y,y,y(x,fdxdy)y,y,y(x,fdxdy)y,y,y(x,fdxdym21mmm2122m2111( )( )( )1122mmy a = s ,ya = s , ya = s ,xa,b写成向量形式：写成向量形式： 0( ,)( ) (7.5.2)dYF x YdxY a =Y,)()()(Yn21n21xy

33、xyxyyyy.s ,s ,sYTn210,)()()()(21212211mmmm,y,yx,yf,y,yx,yf,y,yx,yfx,YF其中其中注意注意: 在形式上在形式上 (7.5.2) 与与 (7.0.1) 一样一样, 所以可所以可以把求解常微分方程初值问题的各种数值方以把求解常微分方程初值问题的各种数值方法推广到方程组上来法推广到方程组上来. 利用向量值函数的微积分理论利用向量值函数的微积分理论, 很容易推导很容易推导出一阶常微分方程组初值问题的数值解法出一阶常微分方程组初值问题的数值解法. )3 . 5 . 7(),(1nnnnYxhFYY如如Euler公式公式,),(T,2, 1

34、nmnnnyyyY其中其中 .),(),(),(),(, 2, 1, 2, 12, 2, 11nmnnnmnmnnnnmnnnnnyyyxfyyyxfyyyxfYxF(7.5.3)的分量形式为的分量形式为 )4 . 5 . 7(. ), 2, 1(),(, 2, 1,1,miyyyxhfyynmnnninini),(),(),(, 2, 1, 2, 12, 2, 11, 2, 11,1, 21, 1nmnnnmnmnnnnmnnnnmnnnmnnyyyxfyyyxfyyyxfhyyyyyy或或 四阶标准的四阶标准的Runge-Kutta公式公式 )5 . 5 . 7()22(643211KKK

35、KhYYnn),(1nnYxFK )2,2(12KhYhxFKnn)2,2(23KhYhxFKnn),(34hKYhxFKnn设设 那么那么 (7.5.5)的分量形式为的分量形式为 ,)4, 3, 2, 1(),(T, 2, 1jkkkKjmjjj四阶标准的四阶标准的Runge-Kutta公式的分量形式公式的分量形式 , )2,2,2,2(),2,2,2,2(),2,2,2,2(),(),22(63 ,3 , 23 , 14,2,2, 22, 13 ,1 ,1 , 21 , 12, 2, 11 ,4,3 ,2,1 ,1,inmininniiinmininniiinmininniinmnnnii

36、iiiininikhykhykhyhxfkkhykhykhyhxfkkhykhykhyhxfkyyyxfkkkkkhyy., 2, 1mi其中其中例例7.5.1：试写出用中点公式解下列初值问：试写出用中点公式解下列初值问题的计算公式题的计算公式. 3)0(, 1)0(,35,643zyzyzzyxy解解:令令 ,zyY那么那么 .nnnzyY.35643),(zyzyxYxF,)2 , 1(, 2, 1jkkKjjj再取再取由向量形式的中点公式由向量形式的中点公式中点公式的向量形式中点公式的向量形式 )2,2(),(12121KhYhxFKYxFKhKYYnnnnnn上述中点公式的分量计算形式

37、为上述中点公式的分量计算形式为2, 22, 111kkhzyzynnnnnnnnnzyzyxkk356431 , 21 , 1)2( 3)2( 5)2(6)2(4)2( 31 , 21 , 11 , 21 , 12, 22, 1khzkhykhzkhyhxkknnnnn分量计算形式为分量计算形式为整理得分量计算格式整理得分量计算格式. 3, 1, )2( 3)2(5, )2(6)2(4)2( 3,35,643,001 , 21 , 12, 21 , 21 , 12, 11 , 21 , 12, 212, 11zyhkzhkykhkzhkyhxkzykzyxkkhzzkhyynnnnnnnnnnnnnn7.5.2 高阶常微分方程初值问题的数值解高阶常微分方程初值问题的数值解法法高阶常微分方程初值问题的一般形式为高阶常微分方程初值问题的一般