数值分析第9章常微分方程数值解法

1、（numerical methods for ordinary differential equations ）问题驱动：蝴蝶效应问题驱动：蝴蝶效应 洛伦兹吸引子洛伦兹吸引子(lorenz attractor)是由是由mit大学的气象学家大学的气象学家edward lorenz在在1963年给出的，他给出第一个混沌现象年给出的，他给出第一个混沌现象蝴蝴蝶效应。蝶效应。 图图10.1.1蝴蝶效应示意图蝴蝶效应示意图 洛伦兹方程是大气流体动力学模型的一个简化的常微分方程组：洛伦兹方程是大气流体动力学模型的一个简化的常微分方程组： dxxydtdyrxyxzdtdzbzxydt 该方程组来源于模拟大

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

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

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

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

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

7、ipschitz条件条件是指是指: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 */。定义定义若某算法的局部截断误差为若某算法的局部截断误差为o(hp+1

11、)，则称该，则称该 算法有算法有p 阶精度。阶精度。定义定义 欧拉法的局部截断误差：欧拉法的局部截断误差：),()()()()()(32112iiiihiiiiiyxhfyhoxyxyhxyyxyr )()(322hoxyih ri 的的主项主项/* leading term */欧拉法具有欧拉法具有 1 阶精阶精度。度。例例1:1: 用欧拉公式求解初值问题用欧拉公式求解初值问题 2201.201yxyxy （）取步长取步长 。 0.1h 解解: : 应用应用eulereuler公式于题给初值问题的具体形式为：公式于题给初值问题的具体形式为： 2120,1,.,1101iiiiyyhx yiy

12、 其中其中 。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.8000000.7352940.6711410.6097560.5524860.5000000.4524890.4098360.009901

13、0.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.002053可用来检验近似解的准确程度。可用来检验近似解的准确程度。 进行计算，数值解已达到了一定的精度。进行计算，数值解已达到了一定的精度。这个初值问题的准确解为这个初值问题的准确解为 ， 21 1y xx从上表最后一列，我们看到取步长从上表最后一列，我们看到取步长0.1h 欧拉公式的改进：欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit euler method */向后差商近似导数向后差商近似导数hxyxyxy)

14、()()(011 x0 x1)(,()(1101xyxfhyxy 由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接同时出现在等式的两边，不能直接得到，故称为得到，故称为隐式隐式 /* implicit */ 欧拉公式，而前者欧拉公式，而前者称为称为显式显式 /* explicit */ 欧拉公式。欧拉公式。111,0,1iiiiyhf xiyyn 一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11()iiiry xy23( )()2ihy xo h即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。

15、阶精度。 梯形公式梯形公式 / /* *trapezoid formula trapezoid formula * */ / 显、隐式两种算法的显、隐式两种算法的平均平均) 1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：的确有局部截断误差的确有局部截断误差 ，311iiiry xyo h即梯形公式即梯形公式具有具有2 阶精度阶精度，比欧拉方法有了进步。，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式公式隐式公式，计算时不得不用到，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。 中点欧拉公式中点欧拉公式 /* m

16、idpoint 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 1: 先用显式欧拉公式作预测，算出先用显式欧拉公式作预测，算出step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 i

18、y )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii1,iiiiyyhf x y111,2iiiiiiyyyyhf xf x注注: :此法亦称为此法亦称为预测预测- -校正法校正法 / /* * predictor-corrector method predictor-corrector method * */ /可以证明该算法可以证明该算法具有具有 2 阶精度阶精度，同时可以看到它，同时可以看到它是个是个单步单步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将看到，它的。后面将看到，它的稳定性高于稳定性高于显式欧拉法。显

19、式欧拉法。改进的欧拉法改进的欧拉法 ),(),(2111 nnnnnnyxfyxfhyy在实际计算时，可将欧拉法与梯形法则相结合，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为计算公式为 ,.2, 1, 0),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyyknnnnnknnnnn应用改进欧拉法应用改进欧拉法, ,如果序列如果序列 收敛收敛, ,)1(1)0(1 nnyy它的极限便满足方程它的极限便满足方程 ),(),(2111 nnnnnnyxfyxfhyy改进欧拉法的截断误差改进欧拉法的截断误差)(0)(311hyxynn 因此，改进欧拉法公式具有因此，改进欧

20、拉法公式具有 2 2 阶精度阶精度例例2:2: 用改进用改进euler公式求解例公式求解例1中的初值问题，中的初值问题， 取步长取步长 。0.1h 解：解：对此初值问题采用改进对此初值问题采用改进eulereuler公式，公式， 其具体形式为其具体形式为 0212111111122()0,1,.,1112piiiicpiiiipciiiyyyhx yyyhxyiyyy 计算结果列于下表：计算结果列于下表：iixiy 1piy 1ciyiiy xy改进的改进的euler法法iiy xyeuler法法01234567891011120.00.10.20.30.40.50.60.70.80.91.0

21、1.11.21.0000000.9900000.9613660.9172460.8619540.8000340.7355270.6175870.6103990.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.4134

22、810.0000000.0000990.0001730.0001850.0001150.0000340.0002330.0004460.0006430.0008030.0009190.0009900.0010230.0000000.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.002053通过计算结果得比较可以看出，改进的通过计算结果得比较可以看出，改进的euler方法方法的计算精度比的计算精度比euler方法要高。方法要高。建立高精度的单步递推格式。建立高精度的单

23、步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，点出发，欧拉法及其各种变形所能达到的最高精度为欧拉法及其各种变形所能达到的最高精度为2 2阶。阶。以以某一斜率某一斜率沿直线达到沿直线达到 点。点。11,iixy 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(2121121211hkyhxfkyxfkkkhyyiiiiii 斜率斜率一定取一定取k1 k2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phkyphxfkyxfkkkh

24、yyiiiiii l ll l首先希望能确定系数首先希望能确定系数 l l1、l l2、p，使得到的算法格，使得到的算法格式有式有2阶精度，即在阶精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hoyxyriii step 1: 将将 k2 在在 ( xi , yi ) 点作点作 taylor 展开展开)(),(),(),(),(2112hoyxfphkyxphfyxfphkyphxfkiiyiixiiii )()()(2hoxyphxyii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx step 2:

25、将将 k2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hoxyphxyhyhoxyphxyxyhyyiiiiiiii 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存在存在

26、无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶阶龙格龙格 - 库塔格式库塔格式。注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 21, 121 l ll lp 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihkhkhkyhxfkhkhkyhxfkhkyhxfkyxfkkkkhyyb bb bb ba ab bb ba ab ba al ll ll l其中其中l li ( i = 1, , m )，a ai (

27、 i = 2, , m ) 和和 b bij ( i = 2, , m; j = 1, , i 1 ) 均为待定系数，确定这些均为待定系数，确定这些系数的步骤与前面相似。系数的步骤与前面相似。 v考虑一阶常微分方程初值问题考虑一阶常微分方程初值问题0( , ),( ),yf x yaxby ay 11()( )( , ( )( )( , ( )iixiixiy xy xf x y x dxy xhfy1( )(, ()nijijijjy xhb f xch y xch11njjb将区域将区域a,ba,b进行分划：进行分划：0,1, 2,.,.jxa xajhjnbahn称 为 步 长 若若11

28、1222(, ()(, ().(, ().iiiinininkf xch y xchkf xc h y xc hkf xc h y xc h(, ()(, ( )( )(, ( )( , ( )jijijijijiijijiikf xc h y xc hf xc h y xhc y xf xc h y xhc f x y x-11jacjmjm则则11(,)jijijmmmf xc h yha k1iy x11111,( ,),(,),2,3,.,niijjjiijjijijmmmyyhb kkf x ykf xc h yha kjnn n级显式级显式runge-kuttarunge-kutt

29、a方法方法11jjmjmac二阶二阶runge-kuttarunge-kutta方法方法取取n=2n=211 12 212221(),( ,),(,).iiiiiiyyh bkb kkf x ykf xc h yc hk( , ),( , ),( , ),( , ),( , ),( , ),iixxiiyyiixxxxiiyyyyiixyxyiiff x yff x yff x yffx yffx yffx y记记由此得由此得12221221231122 21,(,)( , )()( ),()()( )iiiixyiixykfkf xc h yc hkf x yc h fk fo hyybb

30、hfcbh fk fo h另一方面另一方面23123( , ),()( )( )( )( )2()( )2iixyiiiiixyyf x yf yff fhy xy xhy xy xo hhyhfff fo h1222112bbcb为使局部截断误差为为使局部截断误差为 , ,应取应取3()o h改进的改进的eulereuler方法方法v取取112121(/ 2/ 2),(,),(,),0,1,.,1.iiiiiiyyh kkkfxykfxh yhkin2121,12bbc中点方法中点方法取取 12121,(,),(,),220,1,.,1.iiiiiiyyhkkfxyhhkfxykin2121

31、1,0,2bbc二阶二阶heunheun方法方法取取 11212113(),44(,),22(,),330,1,.,1.iiiiiiyyhkkkfxyhhkfxykin212312,443bbc二级二级runge-kuttarunge-kutta方法不超过二阶方法不超过二阶v记记2,2xyxxxyyyfff f gff ff f2232221()2kfc hfc h go h则则22 341122 22 21()( ),2iiyybb hfbc h fbc hg o h2341234()( )( )( )( )()2!3!()().26iiiiiiyhhy xy xhy xy xyxo hhh

32、yhffgf fo h因此局部截断误差只能达到因此局部截断误差只能达到3().o h三级三级runge-kuttarunge-kutta方法方法取取n=3n=311 12233122213331132233132(),( ,),(,),(,),.iiiiiiiiyyh bkb kb kkf x ykf xc h yc hkkf xc h ya hka hkcaa记记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 h33

33、311322331132222333113222331132222332323(,)()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 2 322 23 3()()12() ().2iiyyybbb hfbcbc h fhbc a ffbcbcgo h2341()()( ),26iiyhhy xyhffgf fo h因此要使局部截断误差为因此要使局部截断误差为 ，必须，必须 12322332222333

34、23 21 ,1 / 2 ,1 / 3 ,1 / 6 .bbbbcbcbcbcbca4o hkuttakutta方法方法取取1232332311/6,2/3,1/6,1/2,1,2,1bbbccaa 1123121312121(),636(,),11(,),22(,2).iiiiiiiiyyhkkkkfxykfxh yhkkfxh yhkhk三阶三阶heunheun方法方法v取取1232332311/4,0,3/4,1/3,2/3,2/3,0bbbccaa1131213213(),44(,),11(,),3322(,).33iiiiiiiiyyhkkkfxykfxhyh kkfxhyh k三级

35、三级runge-kuttarunge-kutta方法不超过三阶方法不超过三阶v完全类似于二级完全类似于二级runge-kuttarunge-kutta方法的分析方法的分析只能达到只能达到三级三级runge-kuttarunge-kutta方法的局部截断误差方法的局部截断误差将将 和和 都展开到都展开到 项项, ,易证易证1iy1()iy x4h4()o h 四级四级r-kr-k方法方法取取n=4n=411 1223344122213331132244411422433331324414243(),(,),(,),(,),(,),.iiiiiiiiiiyyh b kb kb kb kkf x y

36、kf xc h yc hkkf 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,aacaaacbbbbb cb cb cb cb cb cb cb cb cb c ab c ac ab c ab c ac ab c c ab c c ac ab c a32431/24.a112341213243(22) / 6,(,),(

37、/ 2,/ 2),(/ 2,/ 2),(,).iiiiiiiiiiyyh kkkkkfxykfxhyhkkfxhyhkkfxh yhk局部截断误差为o(h5)四阶四阶经典龙格经典龙格-库塔法库塔法 /* classical runge-kutta method */ ：附注：附注：二阶二阶runge-kuttarunge-kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到3();o h 五阶五阶runge-kuttarunge-kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 四阶四阶runge-kuttarunge-kutta方法的局部截断误差方法的局部截断误差

38、 只能达到只能达到 三阶三阶runge-kuttarunge-kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 4();o h5();o h6().o h注：注： 龙格龙格-库塔法的主要运算在于计算库塔法的主要运算在于计算 的值，即计的值，即计算算 的值。的值。butcher 于于1965年给出了计算量与可达年给出了计算量与可达到的最高精度阶数的关系：到的最高精度阶数的关系：ikf753可达到的最高精度可达到的最高精度642每步须算每步须算ki 的个数的个数3()o h4()o h5()o h7()o h8()o h6()o h)(2nho8n 由于龙格由于龙格-库塔法的导出基

39、于泰勒展开，故精库塔法的导出基于泰勒展开，故精太好的解，最好采用太好的解，最好采用低阶算法低阶算法而将步长而将步长h 取小取小。度主要受解函数的光滑性影响。对于光滑性不度主要受解函数的光滑性影响。对于光滑性不 收敛性收敛性 /* convergency */定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xi = x0 + i h，当，当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是，则称该算法是收敛收敛的。的。 例：例：就初值问题就初值问题 考察欧拉显式格式的考察欧拉显式格式的收敛性。收敛性。 0)0(yyyyl l解：解：该问题的精确解为该问题的精确解为 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上的解上的解.分别用欧拉显、隐式格式和改进的欧拉格分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。式计算数值解。 1)0()(30)(yxyxy0