第九章 常微分方程的数值解法 ppt课件

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

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

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

4、令等等份份分分成成将将求求解解区区间间nnxNnnhaxNabhNba. ), 1 , 0(,(*)(,()(nnnxyxfxy ,)()()()(1hxyxyxyxynnnn即即用用一一阶阶差差商商近近似似代代替替，将将令令yn为为y(xn)的近似值，将上式代入的近似值，将上式代入(*)式可得式可得1, 2 , 1 , 0 ),(1Nnyxhfyynnnn此式称为欧拉此式称为欧拉(Euler)公式公式.,210Nyyyy值值问问题题的的数数值值解解即即可可由由此此式式依依次次求求出出初初已已知知nnnyxyT)( 为Euler方法的局部截断误差.上一页上一页 下一页下一页 返回返回 例例1

5、用欧拉公式解初值问题用欧拉公式解初值问题 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，y6，y7，y8，y9，y10上一页上一页 下一页下一页 返回返回 其部分结果见下表其部分

6、结果见下表 可见可见Euler方法的计算结果精度不太高。方法的计算结果精度不太高。 xy21 准准确确解解上一页上一页 下一页下一页 返回返回 出出发发的的一一条条曲曲线线表表示示从从的的解解方方程程),()()(),(0000yxPxyyaybxayxfy 欧拉公式的几何意义：欧拉公式的几何意义：x0P0 x1P1),)(,()(),(0000000000 xxyxfyxxxyyyyxP 的的切切线线方方程程过过点点,),(000101hyxfyyhxxx的的交交点点的的纵纵坐坐标标为为与与直直线线 x2P2xnPn2111212111111111,),(),)(,(),(),(yhyxfy

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

8、似值，则可得1, 2 , 1 , 0 ),(111Nnyxhfyynnnn称为后退称为后退Euler公式公式 知知 yn时时,必须通过解方程才能求出必须通过解方程才能求出yn1 ,这样的公式称为这样的公式称为隐式公式隐式公式, 而而Euler公式为显式公式公式为显式公式. Euler公式和后退公式和后退Euler公式都是由公式都是由yn去计算去计算yn+1 ，因而，称，因而，称它们为单步法。它们为单步法。上一页上一页 下一页下一页 返回返回 定义定义 在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提下，考步计算是精确的前提下，考虑的截断误差虑的截断误差 Ti+1 = y(

9、xi+1) yi+1 称为局部截断误差。称为局部截断误差。定义定义 若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。显然显然, p越大越大, 精度越高精度越高.三、三、 局部截断误差与方法的阶局部截断误差与方法的阶（将准确解（将准确解)(xy代入公式的左、右两端，其左端与右端之差）代入公式的左、右两端，其左端与右端之差） Euler方法的精度方法的精度 )(,()()()(111nnnnnnnxyxhfxyhxyyxyT其中：其中：)()(,(nnnxyxyxf 上一页上一页 下一页下一页 返回返回 所以，所以， Euler方法具有

10、方法具有 1 阶精度。阶精度。),( , ! 2)()()()(12 nnnnnxxhyxyhxyhxy )(hxynnx将将在点在点处一阶处一阶TaylorTaylor展开展开)()()(2hOxyhxynn)()()()()()(221hOxyhxyhOxyhxyTnnnnn上一页上一页 下一页下一页 返回返回 所以，所以， 后退的后退的Euler方法也具有方法也具有 1 阶精度。阶精度。 , ! 2)()()()()(2111hyxyhxyhxyxynnnn )(nxy1nx将将在点在点处一阶处一阶TaylorTaylor展开展开 )()()(211hOxyhxynn 隐式隐式Euler

11、方法的精度方法的精度 )( )()()()()(2121111hOxyhhOxyhxyxyTnnnnn上一页上一页 下一页下一页 返回返回 显、隐式两种算法的平均显、隐式两种算法的平均),(),(2111 nnnnnnyxfyxfhyy 欧拉公式的改进欧拉公式的改进)()(1233hOxyhTn 其局部误差为：其局部误差为：此公式具有此公式具有2 2阶精度阶精度. .称平均公式或梯形公式称平均公式或梯形公式梯形公式可由下迭代式计算：其中迭代初值是梯形公式可由下迭代式计算：其中迭代初值是Euler公式提供公式提供. ),1 ,0(),(),(2),(111101kyxfyxfhyyyxhfyyk

12、nnnnnknnnnn上一页上一页 下一页下一页 返回返回 四、改进的欧拉公式四、改进的欧拉公式Step 1: 先用显式欧拉公式作预测，算出先用显式欧拉公式作预测，算出),(1iiiiyxfhyy Step 2: 再将再将 代入隐式梯形公式的右边作校正，得到代入隐式梯形公式的右边作校正，得到1 iy),(),(2111 iiiiiiyxfyxfhyy注：此法亦称为预测注：此法亦称为预测-校正法校正法 。 可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个单阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。它的精步递推格式，比隐式公式的迭代求解过程简单。它的

13、精度高于显式欧拉法。度高于显式欧拉法。 )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii上一页上一页 下一页下一页 返回返回 为了便于编程为了便于编程, 常将改进的欧拉公式写为：常将改进的欧拉公式写为：cpnpnncnnnpyyyyxhfyyyxhfyy21),(),(11上一页上一页 下一页下一页 返回返回 例例2 用改进的欧拉法解例用改进的欧拉法解例1中的初值问题中的初值问题.1)0(10,2yxyxyy解：取步长解：取步长 h=0.1, 改进欧拉法的具体改进欧拉法的具体 形式为形式为2/)()2()2(11cpnpnpncnnnnpyyyyxyhyyy

14、xyhyy 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具体计算过程如下具体计算过程如下上一页上一页 下一页下一页 返回返回 xn改进的欧拉法改进的欧拉法误差误差xn改进的欧拉法改进的欧拉法误差误差0100.61.4859560.0027160.21.1840960.0000880.81.6164760.0040240.41.3433600.0017191.01.7378690.005818

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

16、.1, 改进的欧拉方法的具体公式为改进的欧拉方法的具体公式为取步长取步长h=0.1，分别用，分别用Euler方法、改进的方法、改进的Euler方法求数值解。方法求数值解。上一页上一页 下一页下一页 返回返回 计算结果见下表计算结果见下表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 2970.430 4670.387 4210.348 6790.905 0000.819 0250.741 2180.670

17、 8020.607 0760.549 4040.497 2100.449 9750.407 2280.368 541上一页上一页 下一页下一页 返回返回 .,),(),(11拉拉公公式式高高精精度度比比欧欧值值代代替替平平均均斜斜率率两两点点的的斜斜率率的的平平均均相相当当于于取取点点hkyxyxnnnn第二节第二节 龙格龙格 - 库塔法库塔法基本思想基本思想考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：2/)( ),( ),(2111121kkhyyhkyxfkyxfknnnnnn斜率斜率一定取一定取k1 k2 的平均值吗？的平均值吗？步长一定是一个步长一定是一个h 吗

18、？吗？根根据据微微分分中中值值定定理理的的准准确确解解为为初初值值问问题题设设,)(),()(0 yaybxayxfyxy),(),(,()()()()(11 nnnnnnnnnxxyhfxyhyxyxy .,)()(,(1上上的的平平均均斜斜率率在在区区间间为为解解曲曲线线 nnnnxxxyyf 只要能对平均斜率提供一种近似算法只要能对平均斜率提供一种近似算法,就能得到一种对应的差分格式就能得到一种对应的差分格式.,),(),(精精度度较较低低平平均均斜斜率率近近似似代代替替处处的的斜斜率率点点欧欧拉拉公公式式相相当当于于取取一一个个nnnnyxfyx上一页上一页 下一页下一页 返回返回 .

19、,1公公式式可可能能构构造造出出精精度度更更高高的的代代替替平平均均斜斜率率用用他他们们的的加加权权平平均均值值上上多多预预测测几几个个点点的的斜斜率率在在区区间间 nnxx例如取例如取 m 个点的斜率构造如下形式的公式个点的斜率构造如下形式的公式 )(),(),(),(),(2211111,1123213133121221mmnnmmmmnmnmnnnnnnkckckchyyhkbhkbyhaxfkhkbhkbyhaxfkhkbyhaxfkyxfk,的的常常数数无无关关是是与与其其中中nfcbaiiji该公式称为该公式称为m级龙格库塔级龙格库塔(Runge-Kutta)公式公式,简称简称R-

20、K公式公式.求解：只需将公式的局部截断误差在求解：只需将公式的局部截断误差在xn点进行点进行Taylor展开展开,令其令其前面尽可能多的项为前面尽可能多的项为0, 便可导出便可导出ai，bij，ci所满足的方程组所满足的方程组,即即可从中求出这些系数可从中求出这些系数. 上一页上一页 下一页下一页 返回返回 以以 m=2 的情形为例说明建立的情形为例说明建立R-K公式的方法公式的方法.公公式式级级建建立立令令KRaba 2,212 )(),(),(22111121kckchyyahkyahxfkyxfknnnnnn.,21待待定定其其中中cca其局部截断误差为：其局部截断误差为：)()()()

21、(2211111kckchxyxyyxyTnnnn点点进进行行泰泰勒勒展展开开得得在在nx)()(,(1nnnxyxyxfk )()()(,()()(,(22hOxyyfxfahxyxfxyahxyahxfknnnnnn )(6)(2)()()(321nnnnnxyhxyhxyhxyxy),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx )()()(2hOxyahxynn 上一页上一页 下一页下一页 返回返回 因此有：因此有：)()(21)()(132221hOxyhacxyhccTnn .,2项项系系数数为为零零与与必必须须的的阶阶数数尽尽

22、可可能能高高要要使使hhT 211221accc即满足：即满足： 而对于而对于h3h3，若将，若将k2k2的的TaylorTaylor展开式多取一项展开式多取一项, ,会发现会发现h3h3项的系项的系数不可能为数不可能为0. 0. 而对于上式有无穷多个解而对于上式有无穷多个解, ,它的每一组解都给出了一个局部截它的每一组解都给出了一个局部截断误差为断误差为 的二级的二级R-KR-K公式公式, ,即二阶即二阶R-KR-K公式公式. .)(3hOT 当取当取1,2121acc时时, ,二阶二阶R-KR-K公式就是改进的公式就是改进的EulerEuler公式公式 这里有这里有 个未知个未知数，数，

23、个方程。个方程。32上一页上一页 下一页下一页 返回返回 )22( ),( ),( ),( ),(43216134222312221kkkkyyhkyhxfkkyxfkkyxfkyxfkhnnnnhnhnhnhnnn 常用的标准四阶常用的标准四阶RK公式经典公式经典R-K方法）方法）最常用的四阶标准最常用的四阶标准RK公式经典公式经典R-K方法为：方法为：上一页上一页 下一页下一页 返回返回 例例 用四阶标准用四阶标准R-K公式解初值问题公式解初值问题 1)0(10 ,2yxyxyy解：取解：取 h=0.2 h=0.2 ，四阶标准，四阶标准R-KR-K法的具体格式如下法的具体格式如下: : )

24、22(6 )(2)(),(2)2(2)2()2,2(2)2(2)2()2,2( /2),(432113334222311121kkkkhyyhkyhxhkyhkyhxfkkhyhxkhykhyhxfkkhyhxkhykhyhxfkyxyyxfknnnnnnnnnnnnnnnnnnnnnn上一页上一页 下一页下一页 返回返回 10y5 , 1 , 0 , 2 . 0nnxn知知 183229. 1)22(6 84324. 033849. 0181728. 1)(2)(90864. 018318. 0091818. 12)2(2)2( 91818. 01 . 12 . 01 . 12)2(2- )

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

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

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

28、yy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 上一页上一页 下一页下一页 返回返回 稳定性稳定性 /* Stability */例：考察初值问题例：考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉法、隐式欧拉法和改进的欧拉格式计算数值解。分别用欧拉法、隐式欧拉法和改进的欧拉格式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法隐式欧拉法隐式欧拉法欧拉法欧拉法 节点 xixey30 1.00002.0000 4.00008.0000 1.6000101 3.2000101 1.0000

29、2.5000101 6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107What is wrong ?!上一页上一页 下一页下一页 返回返回 定义定义 若某算法在计算过程中任一步产生的误差在以后的计若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是绝对稳定的算中都逐步衰减，则称该算法是绝对稳定的. 一般分析时为简单起见，只考虑试验方程一般分析时为简单起见，只考虑试验方程yy 常数常数l0，可以是复数可以是复数 当步长取为当步长取为 h 时，将某算法应用于上式