第九章常微分方程的数值解法53180

1、1第九章第九章 常微分方程数值解法常微分方程数值解法第一节第一节 euler方法方法第三节第三节 单步法的收敛性和稳定性单步法的收敛性和稳定性第二节第二节 runge-kutta方法方法上一页上一页 下一页下一页 返回返回 2上一页上一页 下一页下一页 返回返回 本章介绍求解微分方程数值解的基本思想和方法本章介绍求解微分方程数值解的基本思想和方法. 含有自变量、未知函数和它的一阶导数和高阶含有自变量、未知函数和它的一阶导数和高阶导数的方程导数的方程. 常微分方程常微分方程 它是描述运动、变化规律的重要数学方法之一，它是描述运动、变化规律的重要数学方法之一，分为两类：分为两类：1. 初值问题初值

2、问题 即给出未知函数及导数在初始点的值；即给出未知函数及导数在初始点的值；2. 边值问题边值问题 即给出未知函数及（或）它的某些导数即给出未知函数及（或）它的某些导数在区间两个端点的值在区间两个端点的值 。3 考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题 :0)(,),(yaybaxyxfy只要只要 f (x, y) 在在a, b r1 上连续，且关于上连续，且关于 y 满足满足 lipschitz 条条件件，即存在与，即存在与 x, y 无关的常数无关的常数 l 使使对任意定义在对任意定义在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，则上述问题解都成立，则上述问题

3、解存在唯一解存在唯一解。| ),(),(|2121yylyxfyxf 所谓所谓数值解法数值解法就是要计算出初值问题的解函数就是要计算出初值问题的解函数 y(x) 在一系列在一系列离散点离散点 a = x0 x1 xn= b上的近似值：上的近似值： y0 , y1 ,yn .节点间距节点间距 为步长，为步长，) 1,., 0(1 nixxhiii通常采用通常采用等距节点等距节点，即取，即取 hi = h (常数常数)。 yn 称为问题的称为问题的数值解数值解.数值解所满足的离散方程统称为数值解所满足的离散方程统称为差分格式差分格式. 上一页上一页 下一页下一页 返回返回 4第一节第一节 欧拉方法

4、欧拉方法一、一、 欧拉公式欧拉公式点点列列出出方方程程式式在在令令等等份份分分成成将将求求解解区区间间nnxnnnhaxnabhnba. ), 1 , 0(,(*)(,()(nnnxyxfxy ,)()()()(1hxyxyxyxynnnn即即用用一一阶阶差差商商近近似似代代替替，将将令令yn为为y(xn)的近似值，将上式代入的近似值，将上式代入(*)式可得式可得1, 2 , 1 , 0 ),(1nnyxhfyynnnn此式称为此式称为欧拉欧拉(euler)公式公式.,210nyyyy值值问问题题的的数数值值解解即即可可由由此此式式依依次次求求出出初初已已知知nnnyxyt)( 为为euler

5、方法的方法的局部截断误差局部截断误差. .上一页上一页 下一页下一页 返回返回 5例例1 用欧拉公式解初值问题用欧拉公式解初值问题 1)0(10,2yxyxyy解解： 取步长取步长 h=0.1, 欧拉公式的具体形式为：欧拉公式的具体形式为：)2(1 . 0),(1nnnnnnnnyxyyyxfhyy由由此此式式可可得得已已知知其其中中, 1),10, 1 , 0(1 . 00ynnnhaxn 1 . 11 . 01)2(00001 yxyhyy191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(11112 yxyhyy依次计算可得依次计算可得 y3，y4，y5，y

6、6，y7，y8，y9，y10上一页上一页 下一页下一页 返回返回 6其部分结果见下表其部分结果见下表 可见可见euler方法的计算结果精度不太高。方法的计算结果精度不太高。 xy21 准准确确解解上一页上一页 下一页下一页 返回返回 7出出发发的的一一条条曲曲线线表表示示从从的的解解方方程程),()()(),(0000yxpxyyaybxayxfy 欧拉公式的几何意义：欧拉公式的几何意义：x0p0 x1p1),)(,()(),(0000000000 xxyxfyxxxyyyyxp 的的切切线线方方程程过过点点,),(000101hyxfyyhxxx的的交交点点的的纵纵坐坐标标为为与与直直线线

7、x2p2xnpn2111212111111111,),(),)(,(),(),(yhyxfyyhxxxxxyxfyyyxfyxp是是欧欧拉拉公公式式求求出出的的正正好好的的交交点点的的纵纵坐坐标标为为与与的的直直线线为为且且斜斜率率为为过过点点 .1y的的正正好好是是由由欧欧拉拉公公式式求求出出几何意义：几何意义：用折线近似代替方程的解曲线用折线近似代替方程的解曲线,因而也称因而也称euler方法为方法为折线法折线法.上一页上一页 下一页下一页 返回返回 8二、二、 后退的后退的欧拉公式欧拉公式)(,()(1111点点列列出出方方程程在在nnnnxyxfxyx 也用一阶差商逼近导数也用一阶差商

8、逼近导数hxyxyxynnn)()()(11令令yn+1为为y(xn+1)的近似值，则可得的近似值，则可得1, 2 , 1 , 0 ),(111nnyxhfyynnnn称为称为后退后退euler公式公式 已知已知 yn时时,必须通过解方程才能求出必须通过解方程才能求出yn1 ,这样的公式称为这样的公式称为隐式公式隐式公式, 而而euler公式为公式为显式公式显式公式. euler公式和后退公式和后退euler公式都是由公式都是由yn去计算去计算yn+1 ，因此，称，因此，称它们为它们为单步法。单步法。上一页上一页 下一页下一页 返回返回 9定义定义 在假设在假设 yi = y(xi)，即第，即

9、第 i 步计算是精确的前提下，考步计算是精确的前提下，考虑的截断误差虑的截断误差 ti+1 = y(xi+1) yi+1 称为称为局部截断误差。局部截断误差。定义定义 若某算法的局部截断误差为若某算法的局部截断误差为o(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。显然显然, p越大越大, 精度越高精度越高.三、三、 局部截断误差与方法的阶局部截断误差与方法的阶（将准确解（将准确解)(xy代入公式的左、右两端，其左端与右端之差）代入公式的左、右两端，其左端与右端之差） euler方法的精度方法的精度 )(,()()()(111nnnnnnnxyxhfxyhxyyxyt其中：其中：)

10、()(,(nnnxyxyxf 上一页上一页 下一页下一页 返回返回 10所以，所以， euler方法具有方法具有 1 阶精度。阶精度。),( , ! 2)()()()(12 nnnnnxxhyxyhxyhxy )(hxynnx将将在点在点处一阶处一阶taylor展开展开)()()(2hoxyhxynn)()()()()()(221hoxyhxyhoxyhxytnnnnn上一页上一页 下一页下一页 返回返回 11所以，所以， 后退的后退的euler方法也具有方法也具有 1 阶精度。阶精度。 , ! 2)()()()()(2111hyxyhxyhxyxynnnn )(nxy1nx将将在点在点处一阶

11、处一阶taylor展开展开 )()()(211hoxyhxynn 隐式隐式euler方法的精度方法的精度 )( )()()()()(2121111hoxyhhoxyhxyxytnnnnn上一页上一页 下一页下一页 返回返回 12 显、隐式两种算法的显、隐式两种算法的平均平均),(),(2111 nnnnnnyxfyxfhyy 欧拉公式的改进欧拉公式的改进)()(1233hoxyhtn 其局部误差为：其局部误差为：此公式具有此公式具有2 2阶精度阶精度. .称称平均公式或梯形公式平均公式或梯形公式梯形公式可由下迭代式计算：其中迭代初值是梯形公式可由下迭代式计算：其中迭代初值是euler公式提供公

12、式提供. ),1 ,0(),(),(2),(111101kyxfyxfhyyyxhfyyknnnnnknnnnn上一页上一页 下一页下一页 返回返回 13四、四、改进的欧拉公式改进的欧拉公式step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法 。 可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单

13、步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。它的。它的精精度高度高于显式欧拉法。于显式欧拉法。 )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii上一页上一页 下一页下一页 返回返回 14为了便于编程为了便于编程, 常将常将改进的欧拉公式改进的欧拉公式写为：写为：cpnpnncnnnpyyyyxhfyyyxhfyy21),(),(11上一页上一页 下一页下一页 返回返回 15例例2 用用改进的欧拉法解例改进的欧拉法解例1中的初值问题中的初值问题.1)0(10,2yxyxyy解：解：取步长取步长 h=0.1, 改进欧拉法的具

14、体改进欧拉法的具体 形式为形式为2/)()2()2(11cpnpnpncnnnnpyyyyxyhyyyxyhyy 095909. 1)091818. 11 . 1(21)(21091818. 1)1 . 12 . 01 . 1( 1 . 01)2(1 . 0 1 . 1)01( 1 . 01)2(1 . 01100000 yyyyxyyyyxyyycpppcp具体计算过程如下具体计算过程如下上一页上一页 下一页下一页 返回返回 16xn改进的欧拉法改进的欧拉法误差误差xn改进的欧拉法改进的欧拉法误差误差0100.61.4859560.0027160.21.1840960.0000880.81.

15、6164760.0040240.41.3433600.0017191.01.7378690.005818 184096.1)(21180942.1)2(1 .0187250.1)2(1 .02211111cpppcpyyyyxyyyyxyyy依次计算可得依次计算可得 y3，y4，y5，y6，y7，y8，y9，y10其部分结果见下表其部分结果见下表 上一页上一页 下一页下一页 返回返回 17例例3 对下面的对下面的初值问题初值问题1)0(10,yxyy解解 （1）取步长取步长 h=0.1, 欧拉方法的具体公式为欧拉方法的具体公式为nnnnyyhyy9.0)(1 905. 02/ )(91. 0)

16、(1 . 0 9 . 0)(1 . 01ncpnnpncnnnpyyyyyyyyyyyy（2）取步长取步长 h=0.1, 改进的欧拉方法的具体公式为改进的欧拉方法的具体公式为取步长取步长h=0.1，分别用，分别用euler方法、改进的方法、改进的euler方法求数值解。方法求数值解。上一页上一页 下一页下一页 返回返回 18计算结果见下表计算结果见下表euler方法方法改进的改进的euler方法方法xnynyn0.10.20.30.40.50.60.70.80.91.00.900 0000.810 0000.729 0000.656 1000.590 4900.531 4410.478 297

17、0.430 4670.387 4210.348 6790.905 0000.819 0250.741 2180.670 8020.607 0760.549 4040.497 2100.449 9750.407 2280.368 541上一页上一页 下一页下一页 返回返回 19.,),(),(11拉拉公公式式高高精精度度比比欧欧值值代代替替平平均均斜斜率率两两点点的的斜斜率率的的平平均均相相当当于于取取点点hkyxyxnnnn第二节第二节 龙格龙格 - 库塔法库塔法基本思想基本思想考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：2/)( ),( ),(2111121kkhyy

18、hkyxfkyxfknnnnnn斜率斜率一定取一定取k1 k2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？根根据据微微分分中中值值定定理理的的准准确确解解为为初初值值问问题题设设,)(),()(0 yaybxayxfyxy),(),(,()()()()(11 nnnnnnnnnxxyhfxyhyxyxy .,)()(,(1上上的的平平均均斜斜率率在在区区间间为为解解曲曲线线 nnnnxxxyyf 只要能对平均斜率提供一种近似算法只要能对平均斜率提供一种近似算法,就能得到一种对应的差分格式就能得到一种对应的差分格式.,),(),(精精度度较较低低平平均均斜斜率率近近似似代代

19、替替处处的的斜斜率率点点欧欧拉拉公公式式相相当当于于取取一一个个nnnnyxfyx上一页上一页 下一页下一页 返回返回 20.,1公公式式可可能能构构造造出出精精度度更更高高的的代代替替平平均均斜斜率率用用他他们们的的加加权权平平均均值值上上多多预预测测几几个个点点的的斜斜率率在在区区间间 nnxx例如取例如取 m 个点的斜率构造如下形式的公式个点的斜率构造如下形式的公式 )(),(),(),(),(2211111,1123213133121221mmnnmmmmnmnmnnnnnnkckckchyyhkbhkbyhaxfkhkbhkbyhaxfkhkbyhaxfkyxfk,的的常常数数无无关

20、关是是与与其其中中nfcbaiiji该公式称为该公式称为m级级龙格库塔龙格库塔(runge-kutta)公式公式,简称简称r-k公式公式.求解：求解：只需将公式的局部截断误差在只需将公式的局部截断误差在xn点进行点进行taylor展开展开,令其令其前面尽可能多的项为前面尽可能多的项为0, 便可导出便可导出ai，bij，ci所满足的方程组所满足的方程组,即可即可从中求出这些系数从中求出这些系数. 上一页上一页 下一页下一页 返回返回 21以以 m=2 的情形为例说明建立的情形为例说明建立r-k公式的方法公式的方法.公公式式级级建建立立令令kraba 2,212 )(),(),(22111121k

21、ckchyyahkyahxfkyxfknnnnnn.,21待待定定其其中中cca其局部截断误差为：其局部截断误差为：)()()()(2211111kckchxyxyyxytnnnn点点进进行行泰泰勒勒展展开开得得在在nx)()(,(1nnnxyxyxfk )()()(,()()(,(22hoxyyfxfahxyxfxyahxyahxfknnnnnn )(6)(2)()()(321nnnnnxyhxyhxyhxyxy),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx )()()(2hoxyahxynn 上一页上一页 下一页下一页 返回返回 22

22、因此有：因此有：)()(21)()(132221hoxyhacxyhcctnn .,2项项系系数数为为零零与与必必须须的的阶阶数数尽尽可可能能高高要要使使hht 211221accc即满足：即满足： 而对于而对于h3，若将若将k2的的taylor展开式多取一项展开式多取一项,会发现会发现h3项的系数项的系数不可能为不可能为0. 而对于上式有无穷多个解而对于上式有无穷多个解, ,它的每一组解都给出了一个局部截断误它的每一组解都给出了一个局部截断误差为差为 的二级的二级r-k公式公式,即即二阶二阶r-k公式公式.)(3hot 当取当取1,2121acc时时, ,二阶二阶r-k公式就是改进的公式就是

23、改进的euler公式公式 这里有这里有 个未知个未知数，数， 个方程。个方程。32上一页上一页 下一页下一页 返回返回 23)22( ),( ),( ),( ),(43216134222312221kkkkyyhkyhxfkkyxfkkyxfkyxfkhnnnnhnhnhnhnnn 常用的标准四阶常用的标准四阶rk公式公式（经典（经典r-k方法）方法）最常用的四阶标准最常用的四阶标准rk公式（经典公式（经典r-k方法）为：方法）为：上一页上一页 下一页下一页 返回返回 24例例 用四阶标准用四阶标准r-k公式解初值问题公式解初值问题 1)0(10 ,2yxyxyy解：取解：取 h=0.2 ，四

24、阶标准，四阶标准r-k法的具体格式如下法的具体格式如下: )22(6 )(2)(),(2)2(2)2()2,2(2)2(2)2()2,2( /2),(432113334222311121kkkkhyyhkyhxhkyhkyhxfkkhyhxkhykhyhxfkkhyhxkhykhyhxfkyxyyxfknnnnnnnnnnnnnnnnnnnnnn上一页上一页 下一页下一页 返回返回 2510y5 , 1 , 0 , 2 . 0nnxn已知已知 183229. 1)22(6 84324. 033849. 0181728. 1)(2)(90864. 018318. 0091818. 12)2(2)

25、2( 91818. 01 . 12 . 01 . 12)2(2- )2( 1012-4321013003042002031001020001kkkkhyyhkyhxhkykkhyhxkhykkhyhxkhykyxyk上一页上一页 下一页下一页 返回返回 26同理可计算得同理可计算得， , ,5432yyyy具体结果见下表具体结果见下表 至少具有四位有效数字至少具有四位有效数字. 比较比较: :上节用改进的上节用改进的euler公式计算公式计算, ,取取h= =0.1，最多具有四位有，最多具有四位有效数字效数字 。 上一页上一页 下一页下一页 返回返回 27 改进的改进的euler公式每前进一步

26、只要计算两次公式每前进一步只要计算两次f 值值, ,而而4 4阶阶r-k公式每前进一步要计算四次公式每前进一步要计算四次f 值值,但改进的但改进的euler法的步长比法的步长比4阶阶r-k法的小一半法的小一半,两者计算总量差不多两者计算总量差不多. 而而4阶阶r-k法的效果要比改进的法的效果要比改进的euler法好法好. 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。上一页上一页 下一页下一页 返回返回

27、 28第三节第三节 单步法的单步法的收敛性与稳定性收敛性与稳定性 收敛性收敛性 /* convergency */定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xi = x0 + i h，当，当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是，则称该算法是收敛收敛的。的。 例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。 0)0(yyyy 解：解：该问题的精确解为该问题的精确解为 xeyxy 0)( 欧拉公式为欧拉公式为iiiiyhyhyy)1 (1 0)1 (yhyii 对任意固定的对任意固定的 x = xi =

28、i h ，有，有iixhhxihyhyy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 上一页上一页 下一页下一页 返回返回 29 稳定性稳定性 /* stability */例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉法、隐式欧拉法和改进的欧拉格式计算数值解。分别用欧拉法、隐式欧拉法和改进的欧拉格式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法隐式欧拉法隐式欧拉法欧拉法欧拉法 节点节点 xixey30 1.0000 2.0000 4.0000 8.0000 1

29、.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.4788 10 31.2341 10 46.1442 10 63.0590 10 7what is wrong ?!上一页上一页 下一页下一页 返回返回 30定义定义 若某算法在计算过程中任一步产生的误差在以后的计若某算法在计算过程中任一步产生的误差在以后的计算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的. 一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程yy 常数常数 00，可以是复