第9章常微分方程数值解

1、（Numerical Methods for Ordinary Differential Equations ）问题驱动：蝴蝶效应问题驱动：蝴蝶效应 洛伦兹吸引子洛伦兹吸引子(Lorenz attractor)是由是由MIT大学的气象学家大学的气象学家Edward Lorenz在在1963年给出的，他给出第一个混沌现象年给出的，他给出第一个混沌现象蝴蝴蝶效应。蝶效应。 图图9.1.1蝴蝶效应示意图蝴蝶效应示意图 洛伦兹方程是大气流体动力学模型的一个简化的常微分方程组：洛伦兹方程是大气流体动力学模型的一个简化的常微分方程组： dxxydtdyrxyxzdtdzbzxydt 该方程组来源于模拟大气

2、对流，该模型除了在天气预报中有显该方程组来源于模拟大气对流，该模型除了在天气预报中有显著的应用之外，还可以用于研究空气污染和全球侯变化。洛伦著的应用之外，还可以用于研究空气污染和全球侯变化。洛伦兹借助于这个模型，将大气流体运动的强度兹借助于这个模型，将大气流体运动的强度x与水平和垂直方与水平和垂直方向的温度变化向的温度变化y和和z联系了起来。参数联系了起来。参数 称为普兰特数，称为普兰特数， r是规范是规范 化的瑞利数，化的瑞利数， b和几何形状相关。洛伦兹方程是非线性方程组，和几何形状相关。洛伦兹方程是非线性方程组，无法求出解析解，必须使用数值方法求解上述微分方程组。洛无法求出解析解，必须使

3、用数值方法求解上述微分方程组。洛伦兹用数值解绘制结果图伦兹用数值解绘制结果图9.1.1，并发现了混沌现象。，并发现了混沌现象。 微分方程数值解一般可分为：常微分方程数值解和偏微分微分方程数值解一般可分为：常微分方程数值解和偏微分方程数值解。自然界与工程技术中的许多现象，其数学表达式方程数值解。自然界与工程技术中的许多现象，其数学表达式可归结为常微分方程（组）的定解问题。一些偏微分方程问题可归结为常微分方程（组）的定解问题。一些偏微分方程问题也可以转化为常微分方程问题来（近似）求解。也可以转化为常微分方程问题来（近似）求解。Newton最早采最早采用数学方法研究二体问题，其中需要求解的运动方程就

4、是常微用数学方法研究二体问题，其中需要求解的运动方程就是常微分方程。许多著名的数学家，如分方程。许多著名的数学家，如 Bernoulli（家族），（家族），Euler、Gauss、Lagrange和和Laplace等，都遵循历史传统，研究重要等，都遵循历史传统，研究重要的力学问题的数学模型，在这些问题中，许多是常微分方程的的力学问题的数学模型，在这些问题中，许多是常微分方程的求解。作为科学史上的一段佳话，海王星的发现就是通过对常求解。作为科学史上的一段佳话，海王星的发现就是通过对常微分方程的近似计算得到的。本章主要介绍常微分方程数值解微分方程的近似计算得到的。本章主要介绍常微分方程数值解的若干

5、方法。的若干方法。一一、初值问题的数值解法初值问题的数值解法1、一阶常微分方程初值问题的一般形式、一阶常微分方程初值问题的一般形式0( )( , ( ),(1)( )y xf x y xaxby ay常微分方程的数值解法分为常微分方程的数值解法分为（1 1）初值问题的数值解法）初值问题的数值解法 （2 2）边值问题的数值解法）边值问题的数值解法(2) 一般构造方法：一般构造方法：2. 迭代格式的构造迭代格式的构造(1) 构造思想：构造思想：将连续的微分方程及初值条件离散为线性方程将连续的微分方程及初值条件离散为线性方程组加以求解。由于离散化的出发点不同，产生出各种不同的数组加以求解。由于离散化

6、的出发点不同，产生出各种不同的数值方法。基本方法有：有限差分法（数值微分）、有限体积法值方法。基本方法有：有限差分法（数值微分）、有限体积法（数值积分）、有限元法（函数插值）等等。（数值积分）、有限元法（函数插值）等等。 (3) 如何保证迭代公式的稳定性与收敛性如何保证迭代公式的稳定性与收敛性?3. 微分方程的数值解法需要解决的主要问题微分方程的数值解法需要解决的主要问题(1) 如何将微分方程离散化，并建立求其如何将微分方程离散化，并建立求其数值解的迭代公式？数值解的迭代公式？(2) 如何估计迭代公式的局部截断误差与整体误差？如何估计迭代公式的局部截断误差与整体误差？称称 在区域在区域D上对上

7、对 满足满足Lipschitz条件条件是指是指:1212120. .( ,)( ,), , , ( ),( )Ls tf x yf x yL yyxa byyy xy x ( , )f x yy( , ), ( )( )Dx y axb y xyy x记记4、相关定义、相关定义二、初值问题解的存在唯一性二、初值问题解的存在唯一性 考虑一阶常微分方程的初值问题考虑一阶常微分方程的初值问题 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy| ),(),(|2121yyLyxfyxf 则上述则上述IVP存在唯一解。存在唯一解。只要只要 在在 上连续上

8、连续, 且关于且关于 y 满足满足 Lipschitz 条件，条件，( , )f x y1, a bR即存在与即存在与 无关的常数无关的常数 L 使使, x y对任意定义在对任意定义在 上的上的 都成立，都成立，, a b 12,yxyx 求函数求函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值 的方法称为微分方程的数值解法。的方法称为微分方程的数值解法。() (1,., )iiyy xin称节点间距称节点间距 为步长，为步长，通常采用通常采用等距节点等距节点，即取，即取 hi = h (常数常数)。) 1,., 0(1 nixxhiii1,nyy

9、称为微分方程的称为微分方程的数值解数值解。三、初值问题的离散化方法三、初值问题的离散化方法 离散化方法的基本特点是依照某一递推公式，离散化方法的基本特点是依照某一递推公式，值值 ,取取 。按节点从左至右的顺序依次求出按节点从左至右的顺序依次求出 的近似的近似( )iy x(1,., )iyin0y 如果计算如果计算 ，只用到前一步的值，只用到前一步的值 ，则称这则称这类方法为类方法为单步方法单步方法。1iyiy如果计算如果计算 需用到前需用到前r步的值步的值 , ,则称这类方法为则称这类方法为r步方法步方法。1iy11,ii ryy iy 欧拉公式欧拉公式(单步显示公式）：单步显示公式）：向前

10、差商近似导数向前差商近似导数hxyxyxy)()()(010 ),()()()(000001yxfhyxyhxyxy 1y记为记为x0 x1) 1,., 0(),(1 niyxfhyyiiii亦称为亦称为欧拉折线法欧拉折线法 /* Eulers polygonal arc method*/ 在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提步计算是精确的前提下，考虑的截断误差下，考虑的截断误差 Ri = y(xi+1) yi+1 称为称为局部截断局部截断误差误差 /* local truncation error */。定义定义2.2 若某算法的局部截断误差为若某算法的局部

11、截断误差为O(hp+1)，则称该，则称该 算法有算法有p 阶精度。阶精度。定义定义2.1 欧拉法的局部截断误差：欧拉法的局部截断误差：11()iiiRy xy23()()2ihyxO hRi 的的主项主项/* leading term */欧拉法具有欧拉法具有 1 阶精阶精度。度。232 ( )( )( )() ( ,)hiiiiiiy xhy xy xO hyhf x y( )iiyy x( )( , ( )iiiy xf x y x2()O h例例1:1: 用欧拉公式求解初值问题用欧拉公式求解初值问题 2201.201yxyxy （）取步长取步长 。 0.1h 解解: : 应用应用Eule

12、rEuler公式于题给初值问题的具体形式为：公式于题给初值问题的具体形式为： 2120,1,.,1101iiiiyyhx yiy 其中其中 。0.1ixi计算结果列于下表：计算结果列于下表： iixiy iy xiiy xy 1234567891011120.10.20.30.40.50.60.70.80.91.01.11.21.0000000.9800000.9415840.8883890.8252500.7571470.6883540.6220180.5601130.5036420.4529110.4077830.9900990.9615380.9174310.8630690.800000

13、0.7352940.6711410.6097560.5524860.5000000.4524890.4098360.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.002053可用来检验近似解的准确程度。可用来检验近似解的准确程度。 进行计算，数值解已达到了一定的精度。进行计算，数值解已达到了一定的精度。这个初值问题的准确解为这个初值问题的准确解为 ， 21 1y xx从上表最后一列，我们看到取步长从上表最后一列，我们看到取步长0.1h 隐式欧拉法隐式欧拉法(向后向后

14、Euler法）法） /* implicit Euler method */向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy 111,0,1iiiiyhfxiyyn注：注：由于未知数由于未知数 yi+1 同时出现在等式的两边，不能同时出现在等式的两边，不能直接得到，故称为直接得到，故称为隐式隐式 /* implicit */ 欧拉公式，而欧拉公式，而前者称为前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧

15、拉法的局部截断误差：11()iiiRy xy23( )()2ihy xO h即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。 梯形公式梯形公式 / /* *trapezoid formula trapezoid formula * */ / 显、隐式两种算法的显、隐式两种算法的平均平均) 1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：梯形公式的局部截断误差梯形公式的局部截断误差 ，311iiiRy xyO h即梯形公式即梯形公式具有具有2 阶精度阶精度，比欧拉方法有了进步。，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式公式隐式公式，计算时不

16、得不用到，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。欧拉公式的改进欧拉公式的改进 中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假设假设 , 则可以导出则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中

17、点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多一个初值多一个初值, 可能影响精度可能影响精度 Cant you give me a formula with all the advantages yet without any of the disadvantages? Do you think it possible? Well, call me greedy OK, lets make it possible. 改进欧拉法改进欧拉法 /* modified Eulers method */Step

18、 1: 先用显式欧拉公式作预测，算出先用显式欧拉公式作预测，算出Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 iy ) 1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii1,iiiiyyhf x y111,2iiiiiiyyyyhf xf x注注: :此法亦称为此法亦称为预测预测- -校正法校正法 / /* * predictor-corrector method predictor-corrector method * */ /可以证明该算法可以证明该算法具有具有 2 阶精度阶精度，同时可以看到它，同时可以看

19、到它是个是个单步单步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将看到，它的。后面将看到，它的稳定性高于稳定性高于显式欧拉法。显式欧拉法。在实际计算时，可将欧拉法与梯形法则相结合，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为计算公式为 ( )1( )( )111(1)( )( )111( ,)(,)10,1, 2,.2kiiiickiiiikkciiiyyhf x yyyhf xyyyyk应用改进欧拉法应用改进欧拉法, ,如果序列如果序列 收敛收敛, ,(0)(1)11,iiyy它的极限便满足方程它的极限便满足方程111( ,)(,)2iiiiii

20、hyyf x yf xy改进欧拉法的截断误差改进欧拉法的截断误差311()()iiy xyO h因此，改进欧拉法公式具有因此，改进欧拉法公式具有 2 2 阶精度阶精度二元函数的二元函数的n阶阶Taylor展式：展式：*1*1(,)1(,)(,)( , )!knnik iik ikf xh yf xyf xyhOhkx y ( )( , ( )( )( )( )( )( ) ( )( )xyxxxyyxyyyy xf x y xyxff y xyxff y xff y xy xf yx例例2:2: 用改进用改进Euler公式求解例公式求解例1中的初值问题，中的初值问题， 取步长取步长 。0.1h

21、 解：解：对此初值问题采用改进对此初值问题采用改进EulerEuler公式，公式， 其具体形式为其具体形式为 0212111111122()0,1,.,1112piiiicpiiiipciiiyyyhx yyyhxyiyyy 计算结果列于下表：计算结果列于下表：iixiy 1piy 1ciyiiy xy改进的改进的Euler法法iiy xyEuler法法01234567891011120.00.10.20.30.40.50.60.70.80.91.01.11.21.0000000.9900000.9613660.9172460.8619540.8000340.7355270.6175870.6

22、103990.5532890.5009190.4534790.4108591.0000000.9703890.9243970.8667650.8025170.7360290.6706070.6084430.5507850.4981860.4507350.4082370.9800000.9523330.9100950.8571430.7975510.7350250.6725670.6123550.5557930.5036510.4562230.4134810.0000000.0000990.0001730.0001850.0001150.0000340.0002330.0004460.00064

23、30.0008030.0009190.0009900.0010230.0000000.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.002053通过计算结果的比较可以看出，改进的通过计算结果的比较可以看出，改进的Euler方法方法的计算精度比的计算精度比Euler方法要高。方法要高。建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，点出发，欧拉法及其各种变形所能达到的最高精度为欧拉法

24、及其各种变形所能达到的最高精度为2 2阶。阶。以以某一斜率某一斜率沿直线达到沿直线达到 点。点。11,iixy 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii l ll l首先希望能确定系数首先希望能确定系数 l l1、l l2、p，使得到的算法格，使得到的算法格式有式有2阶精度，即在阶

25、精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hOyxyRiii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiii

26、iiiii l ll ll ll ll lStep 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较)()(2)()()(321hOxyhxyhxyxyiiii )()()()(322211hOxyphxyhyyiiii l ll ll l要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii21,1221 pl ll ll l这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶阶龙格龙格 - 库塔格式库塔格式。注意到，注意到， 就是就是改进的欧拉

27、法改进的欧拉法。 21, 121 l ll lp二阶中点方法二阶中点方法取取 12121,(,),(,)22iiiiiiyyhkkfxyhhkfxyk1210,1,2pll二阶二阶HeunHeun方法方法取取 11212113(),44(,),22(,),330,1,.,1.iiiiiiyyhkkkfxyhhkfxykiN12132,443pllv记记2,2xyxxxyyyFff f Gff ff f22321()2kfphFp h GO h则则2234112221()()2iiyyhfph Fp h G O hllll23411()24iyhfh Fph GO h二级二级Runge-Kutt

28、aRunge-Kutta格式的精度格式的精度2341234()( )( )( )( )()2!3!()().26iiiiiiyhhy xy xhy xy xyxO hhhyhfFGf FO h因此局部截断误差只能达到因此局部截断误差只能达到3().O h二级二级Runge-KuttaRunge-Kutta方法的精度方法的精度不超过二阶！不超过二阶！ 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxf

29、KyxfKKKKhyyb bb bb ba ab bb ba ab ba al ll ll l其中其中l li ( i = 1, , m )，a ai ( i = 2, , m ) 和和 b bij ( i = 2, , m; j = 1, , i 1 ) 均为待定系数，确定这些均为待定系数，确定这些系数的步骤与前面相似。系数的步骤与前面相似。 11111,( ,),(,),2,3,.,niijjjiijjijijmmmyyhkkf x ykf xc h yha kjnln n级显式级显式Runge-KuttaRunge-Kutta公式公式11jjmjmac三级三级Runge-KuttaRun

30、ge-Kutta方法方法取取n=3n=3，得，得11 12233122213331132233132(),( ,),(,),(,),.iiiiiiiiyyhkkkkf x ykf xc h yc hkkf xc h ya hka hkcaalll记记2,2xyxxxyyyFff f Gff ff f222122322(,)1()2iikfxc h yc hkfc hFc h GO h31 132222332322321()2a ka kc fc a hFc a h GO h三级三级Runge-KuttaRunge-Kutta方法的精度方法的精度则则33311322331132222333113

31、222331132222332323(,)()12()2()()1()()2iixyxxxyyyykf xc h ya hka hkfh c fa ka kfhc fca ka kfa ka kfO hfc hFhc a F fc GO hv又由于又由于211232 23 332243 2322 23 3()()12() ().2iiyyyhfcc h Fhc a FfccGO hllllllll2341()()( ),26iiyhhy xyhfFGf FO h因此要使局部截断误差为因此要使局部截断误差为 ，必须，必须 123223322223332321,1/ 2,1/3,1/ 6.cccc

32、c allllllll4O h三阶三阶KuttaKutta方法方法取取1232332311/6,2/3,1/6,1/2,1,2,1ccaalll 1123121312121(),636(,),11(,),22(,2).iiiiiiiiyyhkkkkf xykf xh yhkkf xh yhkhk得得三阶三阶HeunHeun方法方法v取取1232332311/4,0,3/4,1/3,2/3,2/3,0ccaalll1131213213(),44(,),11(,),3322(,).33iiiiiiiiyyhkkkfxykfxh yhkkfxh yhk得得三级三级Runge-KuttaRunge-K

33、utta方法的精度不超过三阶方法的精度不超过三阶v完全类似于二级完全类似于二级Runge-KuttaRunge-Kutta方法的分析方法的分析三级三级Runge-KuttaRunge-Kutta方法的局部截断误差方法的局部截断误差只能达到只能达到将将 和和 都展开到都展开到 项项, ,易证易证1iy1()iy x4h4()O h 四级四级R-KR-K方法方法取取n=4,n=4,得得11 1223344122213331132244411422433331324414243(),(,),(,),(,),(,),.lllliiiiiiiiiiyyhkkkkkf x ykf xc h yc hkkf

34、 xc h ya hka hkkf xc h ya hka hka hkcaacaaa2233132341424341234223 344222333223 344223 3443 2324242343222332442433 2 3324424234342,1,1/2,1/3,1/4,()1/6,()1/12,()1/8,llllllllllllllllllllaacaaaccccccccccc ac ac ac ac ac ac c ac c ac ac a32431/24.a其中，其中，112341213243(22),6(,),(,),22(,),22(,).iiiiiiiiiihyy

35、kkkkkfxyhhkfxykhhkfxykkfxh yhk局部截断误局部截断误差为差为O(hO(h5 5) )四阶四阶经典龙格经典龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ：例例3 用四阶用四阶R-K方法求初值问题方法求初值问题 1( , ), (0)12f x yyy 取步长取步长 1h ，求，求 (1)y。解：取解：取 00100,1,1xyxxh，计算四个组成部分，计算四个组成部分 1213243(,) (,) 22(,)22(,) nnnnnnnnkf xyhhkf xykhhkf xykkf xh yhk将其代入将其代入 11234(

36、22)6nnhyykkkk得得 111( 0.3333 2 0.35292 0.35420.3780)0.645744444444446y 即使用四阶即使用四阶R-K方法在方法在 1x 处的估计值为处的估计值为 0.6457444444444。注：注：该问题有一个明显的解析解该问题有一个明显的解析解( )922y xx(1)720.64575131106459y易知易知11(1)0.0000068666Ryy截断误差为截断误差为附注：附注：二阶二阶Runge-KuttaRunge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到3();O h 五阶五阶Runge-KuttaRu

37、nge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 四阶四阶Runge-KuttaRunge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 三阶三阶Runge-KuttaRunge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 4();O h5();O h6().O h注：注： 龙格龙格-库塔法的主要运算在于计算库塔法的主要运算在于计算 的值，即计的值，即计算算 的值。的值。Butcher 于于1965年给出了计算量与可达年给出了计算量与可达到的最高精度阶数的关系：到的最高精度阶数的关系：iKf753可达到的最高精度可达到的最高精度

38、642每步须算每步须算Ki 的个数的个数3()O h4()O h5()O h7()O h8()O h6()O h)(2nhO8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精库塔法的导出基于泰勒展开，故精太好的解，最好采用太好的解，最好采用低阶算法低阶算法而将步长而将步长h 取小取小。度主要受解函数的光滑性影响。对于光滑性不度主要受解函数的光滑性影响。对于光滑性不 收敛性收敛性 /* Convergency */定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xi = x0 + i h，当，当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是，则称该

39、算法是收敛收敛的。的。 例例4：就初值问题就初值问题 考察欧拉显式格式的考察欧拉显式格式的 收敛性。收敛性。0(0)yyyyl 解：解：该问题的精确解为该问题的精确解为 xeyxyl l0)( 欧拉公式为欧拉公式为iiiiyhyhyy)1 (1l ll l 对任意固定的对任意固定的 x = xi = i h ，有，有iixhhxihyhyyl ll ll ll l)1()1(/10/0 ehhhl ll l/10)1(lim)(0ixxyeyi l l 稳定性稳定性 /* Stability */例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解上的解.分别用欧拉显、隐式格式和改

40、进的欧拉格分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.478

41、8 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong ?! An Engineer complains: Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!定义定义若某算法在计算过程中任一步产生的误差在若某算法在计算过程中任一步产生的误差在以后的计算中都以后的计算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定绝对稳定的的 /*absolutely stable

42、*/。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程 /* test equation */yyl l 常数，可以常数，可以是复数是复数我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定的绝对稳定区域比区域比 B 的的大大。当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在时，将某算法应用于上式，并假设只在初值产生误差初值产生误差 ，则若此误差以后逐步衰减，则若此误差以后逐步衰减,000yy 就称该算法相对于就称该算法相对于 绝对稳定，绝对稳定， 的全体构成的全体构成hl l h h绝对稳定区域绝对稳定区域。例例5：考察显式欧拉

43、法考察显式欧拉法110(1)iiiiyyhyhyll000yy 011)1 (yhyii 01111)1 ( iiiihyy0-1-2ReImg由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：hhl l 1|1| h例例6：考察隐式欧拉法考察隐式欧拉法11 iiiyhyyl liiyhy 11101111 iih可见绝对稳定区域为：可见绝对稳定区域为：1|1 | h210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。的显式法的好。例例7：隐式龙格隐式龙格-库塔法库塔法 ),.,

44、1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmiib bb ba al ll l其中其中2阶方法阶方法 的绝对稳定的绝对稳定区域为区域为 )2,2(1111KhyhxfKhKyyiiii0ReImg无条件稳定无条件稳定而而显式显式 1 4 阶方法的绝对稳定区域为阶方法的绝对稳定区域为k=1k=2k=3k=4-1-2-3-123ReImg5* 线性多步法（阅读材料）线性多步法（阅读材料） 通常可以采用线性多步法，它的构造思想是利用前面通常可以采用线性多步法，它的构造思想是利用前面多步信息进行计算，以获得较高精度的数值公式。多步信息进行计算，以获得较高精度的数值公式。

45、多步法中最常用的是线性多步法。多步法中最常用的是线性多步法。 一、一般迭代格式一、一般迭代格式设线性多步法的一般迭代格式为：设线性多步法的一般迭代格式为：101rrnin iin iiiyyhfab其中，其中， ,iia b均为常数，均为常数， (,) (1, ,)kkkff xyknnnr。10b为隐式。为隐式。 构造线性多步法的主要途径有基于数值积分法和基于构造线性多步法的主要途径有基于数值积分法和基于Taylor展开方法。这里介绍基于展开方法。这里介绍基于Taylor展开式的构造方法，即将线展开式的构造方法，即将线性多步公式在性多步公式在 nx处处Taylor展开，与展开，与 ()ny

46、x在在 nx处处Taylor展开式相展开式相比较，要求它们前面有尽可能多的项重合，由此确定参数比较，要求它们前面有尽可能多的项重合，由此确定参数 ,iia b。下面介绍几种常见的线性多步法。下面介绍几种常见的线性多步法。 10b 为显式，为显式， 若若01011()()1(1,2,3 4)riirrkkiiiiikikaab，可解得可解得 001235559379=1,24242424abbbb ，相应的线性多步公式为，相应的线性多步公式为1123(5559379)24nnnnnnhyyffff因因 1=0b，此式称为，此式称为Adams显式公式，是四阶公式。局部截断误显式公式，是四阶公式。局

47、部截断误差为差为53354(5)65(5)61112511()5()()()5!720niinniihRiiyO hh yO hab如果令如果令 12330aaab，由方程组可解得，由方程组可解得 019 =1,24ab，0121951,242424bbb 二、二、 Adams 格式格式取取 3r ，并令，并令 12310aaab，由方程组，由方程组 相应的线性多步公式为相应的线性多步公式为 1112(9195)24nnnnnnhyyffff称为称为四阶的四阶的Adams隐式公式隐式公式，其，其局部截断误差局部截断误差为为5(5)6119()720nnRh yO h 三、三、 Milne公式公

48、式取取 3r ，并令，并令 01210aaab，由方程组解得，由方程组解得 31a083b，143b 283b30b，相应的线性多步公式为，相应的线性多步公式为 13124 (22)3nnnnnyyhfff称为称为Milne公式公式，其，其局部截断误差局部截断误差为为 5(5)6114()45nnRh yO h,Milne公式是四阶四步显式公式。公式是四阶四步显式公式。 四、四、 Hamming公式公式取取 2r ，并令，并令 120ab，可得到，可得到Hamming公式公式 121113 (9)(2)88nnnnnnyyyh fff其局部截断误差为其局部截断误差为 5(5)611()40nn

49、Rh yO h ，Hamming公式是四阶三公式是四阶三步隐式公式。注意若只有一个初值，不能用多步法，需先用单步隐式公式。注意若只有一个初值，不能用多步法，需先用单步法算出步法算出 个值，再用多步法。个值，再用多步法。 1r 1ny 用显式公式计算预测值，然后用隐式公式进行校正，得到用显式公式计算预测值，然后用隐式公式进行校正，得到近似值近似值 ，这样的一组计算公式称为预测，这样的一组计算公式称为预测校正格式。一般校正格式。一般采用同阶的隐式公式与显式公式，常用的预测采用同阶的隐式公式与显式公式，常用的预测校正系统有两校正系统有两种：种：123111121(5559379) 249 (,) 1

50、95 24nnnnnnnnnnnnnhyyffffhyyf xyfff1Adams预测预测-校正格式校正格式2Milne- Hamming预测预测-校正格式校正格式3121121114(22)313(9) (,)2)88nnnnnnnnnnnnyyhfffyyyh f xyff 以上两种预测以上两种预测-校正格式均为四阶公式，其起步值通常用校正格式均为四阶公式，其起步值通常用四阶四阶R-K方法计算。有时为了提高精度，校正格式可以多次方法计算。有时为了提高精度，校正格式可以多次迭代，但通常迭代次数不超过迭代，但通常迭代次数不超过3次次.例例8 用用Adams预测预测-校正格式求解初值问题校正格式

51、求解初值问题2, (01)(0)1 xyyxyy ，先用四阶，先用四阶Runge-Kutta方法计算初始值方法计算初始值 123,y yy，然后使用然后使用Adams预测预测-校正格式，计算结果如表校正格式，计算结果如表9.5.1所示所示 表表9.5.1解：取步长解：取步长 0.1h xi0.00.10.20.30.40.50.60.70.80.91.01.34151.41411.48321.54911.61241.67331.7320 yi1.00001.05941.18321.26491.34161.41421.48321.54921.61241.67331.7320y(xi)1.0000

52、1.05941.18321.26491.34161.41421.48321.54921.61251.67331.7321iy其中，其中， iy和和 iy分别为预测值和校正值，分别为预测值和校正值， ( )iy x是准确是准确值，值，比较显示，计算结果相当精确。比较显示，计算结果相当精确。6 方程组与高阶方程的解法方程组与高阶方程的解法 一阶微分方程初值问题的各种数值解法，同样适用于一阶一阶微分方程初值问题的各种数值解法，同样适用于一阶微分方程组与高阶方程。下面介绍含有两个未知函数的方程组微分方程组与高阶方程。下面介绍含有两个未知函数的方程组和二阶方程的解法，可以将其推广到更高阶以及更大规模的方

53、和二阶方程的解法，可以将其推广到更高阶以及更大规模的方程组。程组。一、一、 一阶方程组一阶方程组设有一阶方程组初值问题设有一阶方程组初值问题0000( , , ), ()( , , ), ()dyf x y zy xydxdzg x y zz xzdx1欧拉法的计算格式欧拉法的计算格式110000(,)(,) 0,1,2,(), ()nnnnnnnnnnyyhf xyzzzhg xyznyy xzz x2改进欧拉法的预测改进欧拉法的预测-校正格式校正格式1111(0)(0)(1)( )( )111(1)(,) (,) (,)(,)2 (,)(2nnnnnnnnnnnnkkknnnnnnnknn

54、nnyyhf xyzzzhg xyzhyyf xyzf xyzhzzg xyzg( )( )111 0,1,2,)0,1,2, kknnnnxyzk3四阶标准四阶标准R-K公式公式11234112341111211222322341(22)61(22)6(,), (,)(,), 222(,),222(,),222(,),222(nnnnnnnnnnnnnnnnnnnnnnnyykkkkzzmmmmkhf xyzmhg xyzkmhkhf xyzkmhmhg xyzkmhkhf xyzkmhmhg xyzkhf x33433 0,1,2,),(,).nnnnnnh yk zmmhg xh yk

55、zm4四阶四阶Adams显式格式显式格式11231123(5559379)24 4,5,(5559379)24nnnnnnnnnnnnhyyffffnhzzgggg二、二、 高阶方程高阶方程1迭代格式的构造迭代格式的构造(1)基本思想：化高阶微分方程为一阶微分方程组，然后采用某基本思想：化高阶微分方程为一阶微分方程组，然后采用某种数值解法求出解函数及解函数各阶导数的近似值。种数值解法求出解函数及解函数各阶导数的近似值。(2) 构造方法：以二阶导数为例说明求解高阶微分方程的方法，构造方法：以二阶导数为例说明求解高阶微分方程的方法，并采用古典形式的龙格并采用古典形式的龙格-库塔公式求解。库塔公式求

56、解。0000( , ,) (), ()yf x y yy xyy xy设二阶微分方程的初值问题为设二阶微分方程的初值问题为引进新的变量引进新的变量 zy，即可化为下列一阶方程组的初值问题：，即可化为下列一阶方程组的初值问题： 0000, ( , , ) (), () yz zf x y zy xyy xy利用四阶经典形式的龙格利用四阶经典形式的龙格-库塔公式求解方程组有库塔公式求解方程组有1123411234()6(22)6nnnnhyykkkkhzzllll以及以及1213243 22nnnnkzhkzlhkzlkzhl1211322233(,) (,)222(,)222(,)nnnnnnn

57、nnnnnlf xy zhhhlf xyk zlhhhlf xyk zllf xh yhk zhl，消去消去 1234,k kkk，则上述格式可表示为，则上述格式可表示为2112311234()6(22)6nnnnnhyyhzlllhzzllll其中，其中，12123122223(,) (,) 222(,)2242(,)2nnnnnnnnnnnnnnnlf xyzhhhlf xyzzlhhhhlf xyzl zlhlf xh yhzlzhlLorenz方程组求解方程组求解 使用四阶使用四阶Ronge-kutta法求解洛伦兹方程，其中参数值取法法求解洛伦兹方程，其中参数值取法如下：如下： 102

58、.666667b ,和和 28r 0005xyz。，积分区间，积分区间 0,20t。1010 28 2.666667 dxxydtdyxyxzdtdzzxydt 解：将初始条件代入四阶解：将初始条件代入四阶Ronge-kutta求解公式，步长取求解公式，步长取 0.5h 11234112341123411112121(22)61(22)61(22)6( 1010), (28),( 2.666667) 10()10(),22 28()()()(222nnnnnnnnnnnnnnnnnnnnnxxkkkkyysssszzmmmmkhxyshxyx zmhzx ykhkhxykmhhshxyxz1112232232234343)2 2.666667()