第八章 常微分方程数值解法 - 副本

1、1/47第八章 常微分方程数值解 1 引言 2 欧拉方法 3 龙格库塔方法 4 阿达姆斯方法 小结 作业与实验2/47本章要求 1. 熟悉Euler显格式,梯形法及Euler预校法; 2. 熟悉局部截断误差及绝对稳定性; 3. 掌握龙格库塔法。3/47 一. 问题提出 有一个或多个导数及其函数的方程式称为微分方程，在工程中常遇到求解微分方程的问题。1 引言)(1)(,),(00 xyyyxybaxyxfdxdy 解解函函数数式式初初值值问问题题考考虑虑一一阶阶常常微微分分方方程程的的自自变变量量）偏偏微微分分方方程程（一一个个以以上上量量）常常微微分分方方程程（一一个个自自变变微微分分方方程程

2、4/47 很多微分方程的解不能用初等函数来表示，有时很多微分方程的解不能用初等函数来表示，有时即使能够用解析式表示其解，但计算量太大而不实即使能够用解析式表示其解，但计算量太大而不实用用(表达式过于复杂表达式过于复杂)。 需要用数值方法来求解，一般只要求得到若干个需要用数值方法来求解，一般只要求得到若干个点上的近似值或者解的简单的近似表达式点上的近似值或者解的简单的近似表达式(精度要精度要求满足即可求满足即可)。 计计算算）数数值值解解（适适合合于于计计算算机机近近似似解解析析解解求求近近似似解解求求精精确确解解一一般般极极困困难难，1 引言5/471 引言 二. 两类求定解问题 实际中求解常

3、微分方程的所谓定解问题有两类:初值问题和边值问题 定解指已知因变量和/或其导数在某些点上是已知的(约束条件)。 1. 边解问题 约束条件为已知，在自变量的任一非初值上，已知函数值和/或其导数值，如 常常可以将边解问题转化为初值问题求解。ybyayyyxfy求求解解 )(,)(),(6/471 引言 2. 初值问题得得到到。，这这可可以以通通过过递递推推公公式式，求求出出知知的的进进，由由已已点点：按按节节点点顺顺序序依依次次推推初初值值问问题题的的数数值值解解法法特特，有有通通常常取取等等距距节节点点，即即。的的近近似似值值上上的的在在，以以满满足足上上述述两两式式，即即求求解解：初初值值上上

4、已已知知函函数数值值，如如约约束束条条件件为为在在自自变变量量的的110011000000,), 2 , 1 , 0(), 2 , 1 , 0()()()(/),(/)(),( iiiiiiinyyyyniihxxxxhniyxybxxxaxyyxyxxyxfdxdyyxyyxfy7/47初值初值问题问题的的常见常见解法解法单步法：单步法： 利用前一个单步的信息利用前一个单步的信息(一个点一个点)，在，在y=f(x) 上找下一点上找下一点yi， 有欧拉法，龙格库格法。有欧拉法，龙格库格法。预测校正法：预测校正法： 多步法，利用一个以上的前点信息求多步法，利用一个以上的前点信息求f(x)上的下一

5、个上的下一个yi， 常用迭代法，如改进欧拉法，阿当姆常用迭代法，如改进欧拉法，阿当姆斯法。斯法。1 引言8/47yx( )yf x 0y0 x0hxxxxx)()()()()()()1(1211211210nnnnnnnnxyyyyyyxyxyxyxyxxxxxxy 的的近近似似解解上上的的值值在在一一些些离离散散点点的的数数值值解解，即即解解函函数数式式本本章章：求求1 引言9/47 3)1(42),(3)1(42),(4242yxyyxfdxdyyxyyxfxyyxyyx的的初初值值问问题题：就就得得到到一一阶阶常常微微分分方方程程，且且给给出出初初值值令令：例例：方方程程1 引言10/4

6、7唯唯一一。初初值值问问题题的的解解必必存存在在且且由由常常微微分分方方程程理理论论知知：，使使得得即即存存在在常常数数条条件件，满满足足李李普普希希兹兹且且关关于于适适当当光光滑滑连连续续，只只要要函函数数yyLyxfyxfLLipschitzyyxf ),(),()(),(1 引言11/47 一. 欧拉（Euler）格式展展开开法法（略略）实实质质还还是是方方法法二二：数数值值微微分分法法显显格格式式，可可循循环环求求解解。叫叫式式式式或或展展开开法法方方法法一一：式式中中设设节节点点为为TaylorEulerniyxhfyyyhxyxhfxyxxyxxxyxyxyTaylornabhni

7、ihaxiiiiiiiiiiiiiiiii22)1,2 , 1 , 0(),()(2)(,()()(!2)()()()(:),2 , 1 , 0(122111 可忽略可忽略2 欧拉方法12/47),(),()()(3)(,()()()(,(,),(1111111iiiiiiiixxiixxxxiiyxhfyyyxhfxyxydxxyxfxyxydxxyxfdxyxxyxfdxdyiiiiii 求积公式，得求积公式，得对右端积分采用左矩形对右端积分采用左矩形式式积分积分两端于两端于将方程将方程方法三：数值积分法方法三：数值积分法左矩形左矩形2 欧拉方法13/47),()(,(),(),(),()

8、(,(),()(,()()(),()(,(11122211111111000111000000000000000yxhfyyyxxxxyxfyyyxfyxyxhfyyyxxxxyxfyyyxfkxyxfxykxxkyyyxxyxEuler ，则则纵纵坐坐标标记记为为求求交交点点，与与为为斜斜率率作作直直线线以以过过，则则纵纵坐坐标标记记为为求求交交点点，与与斜斜率率作作切切线线过过法法的的几几何何意意义义方方法法四四：几几何何方方法法2 欧拉方法14/47( )yy x 0yx012 xxx012 y yy2 欧拉方法15/47 欧拉公式：欧拉公式：/* Eulers Method */向前差

9、商近似导数向前差商近似导数hxyxyxy)()()(010 ),()()()(000001yxfhyxyhxyxy 1y记为记为)1,., 0(),(1 niyxfhyyiiii2 欧拉方法16/472 欧拉方法 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商近似导数向后差商近似导数hxyxyxy)()()(011 )(,()(1101xyxfhyxy )1,., 0(),(111 niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得同时出现在等式的两边，不能直接得到，故称为到，故称为隐式隐式 /* implicit *

10、/ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。17/47)5 , 4 , 3 , 2 , 1 , 0()1(2 . 0),()(2 . 01),(42 . 00)0()10(111 iyxyyxhfyyxyhxyyxfhyxxydxdyiiiiiiii欧拉公式为：欧拉公式为：，解：因为解：因为位小数。位小数。计算过程保留。计算过程保留取步长取步长题题例：用欧拉法解初值问例：用欧拉法解初值问2 欧拉方法18/477.2 欧拉方法当当i=0，即，即x0=0时，时，

11、y1=0+0.2(1-00)=0.2000当当i=1 ，即，即x1=0.2时，时， y2=0.2+0.2(1-0.20.2)=0.3920当当i=2 ，即，即x2=0.4时，时， y3=0.392+0.2(1-0.40.392)=0.56064当当i=3 ，即，即x3=0.6时，时， y4=0.56064+0.2(1-0.60.56064)=0.6933632当当i=4 ，即，即x4=0.8时，时， y4=0.6933632+0.2(10.80.6933632)=0.782425088列表计算如下：列表计算如下：i012345xi00.20.40.60.81yi0.00000.2000.392

12、00.56060.69340.782419/47，计计算算结结果果如如下下表表：分分别别取取步步长长公公式式为为解解：此此时时的的。确确解解是是的的数数值值解解，此此问问题题的的精精，方方法法求求初初值值问问题题例例：利利用用05. 0 , 1 . 0 , 2 . 02 , 1 , 0,0)211()1/()(0)0(202110221222 hiyyxhyyEulerxxxyyxyxyEuleriiii自己看看是不是自己看看是不是2 欧拉方法20/477.2 欧拉方法hxiyiy(xi)y(xi)-yih=0.20.000.400.801.201.602.000.000000.376310.

13、542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.201.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.6

14、02.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.0022721/47显显格格式式确确定定。，初初值值用用、用用迭迭代代法法求求情情况况）、化化为为显显格格式式（仅仅小小数数是是隐隐格格式式。实实用用上上式式式式（形形求求积积公公式式得得：中中，对对右右端端积积分分采采用用梯梯式式在在Eulerynyxfyxfhyyxyxfxyxfhxyxyiiiiiiiiiiiii1111111 2

15、14)(4)( )1,0,1,2,(i ),),2)(,)(,2)()()3( 0)的近似值的近似值yk 。abx1x2x3积分曲线积分曲线bxxxxann 11040/477.3 龙格 库塔方法折线折线法法法法的几的几何构何构造造 Iabx1x2x3Euler法法(折线法折线法)斜率近似斜率近似 y(a)y(x1)y(x2)y(x3)前进前进Euler法的计算公式法的计算公式)1,2 , 1 , 0(),(1 niyxhfyyiiii显式公式显式公式41/477.3 龙格 库塔方法折线折线法法法法的几的几何构何构造造IIabx1x2x3斜率近似斜率近似y(x1)y(x2)y(x3)y(x4)

16、后退后退Euler法的公式法的公式)1,2 , 1 , 0(),(111 niyxhfyyiiii隐式公式隐式公式42/47单步法单步法：即利用前一个节点的函数值：即利用前一个节点的函数值yi，计算后计算后一个节点的函数值一个节点的函数值yi+1。目的：建立高精度的单步递推格式。目的：建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜率某一斜率沿直线达到沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其点。欧拉法及其各种变形所能达到的最高精度为各种变形所能达到的最高精度为2 阶。阶。3 龙格 库塔方法 二. 2

17、阶龙格 库塔格式43/47斜率斜率一定取一定取K1 K2的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？ )(),(),(2121),(),(2),(00121211111111xyyhKyhxfKyxfKKKhyyyxfyxfhyyyxhfyyiiiiiikkkkkkkkkk可可将将其其改改写写为为：考考察察改改进进的的欧欧拉拉法法，3 龙格 库塔方法44/477.3 龙格 库塔方法首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hOyxy

18、Riii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx 45/47Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphx

19、yhyhOxyphxyxyhyyiiiiiiii 3 龙格 库塔方法46/47Step 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii 21,1221 p 这里有这里有3个未知个未知数，数，2个方程。个方程。存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙阶龙格格 - 库塔格式库塔格式。21, 121 p注意到，注意到，

20、就是就是改进的欧拉法改进的欧拉法。 3 龙格 库塔方法47/47其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待定系数，确定这些系数的步骤均为待定系数，确定这些系数的步骤与前面相似。与前面相似。 ).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 问题问题: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？3 龙格 库

21、塔方法 三. 高阶龙格 库塔格式48/473阶阶龙格龙格-库塔法库塔法)2,(),(),()4(2131222132161hKhKyhxfKKyxfKyxfKKKKyyiihihiiihii 3 龙格 库塔方法49/473 龙格 库塔方法最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 50/47注：注： 龙格龙格-库塔法的主要运算在于计算库塔法的主要运算在于计

22、算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2 nhO8 n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。3 龙格 库塔方法51/477.

23、3 龙格 库塔方法例：例：隐式龙格隐式龙格-库塔法库塔法 ),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii 而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为 )2,2(1111KhyhxfKhKyyiiii其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定52/477.3 龙格 库塔方法例：使用高阶例：使用高阶R-K方法计算初值问题方法计算初值问题解解: （1）使用三阶）使用三阶R-K方法方法(公式公式) 1)0(5 . 002yxyy. 1

24、. 0 h取取1111. 1)4(61 . 02555. 1)2( 1 . 0(1025. 1)21 . 0(1132101212032102201 KKKyyKKyKKyKyKn时时53/477.3 龙格 库塔方法其余结果如下其余结果如下: i xi k1 k2 k3 yi 1.0000 0.1000 1.0000 1.1025 1.2555 1.1111 2.0000 0.2000 1.2345 1.3755 1.5945 1.2499 3.0000 0.3000 1.5624 1.7637 2.0922 1.4284 4.0000 0.4000 2.0404 2.3423 2.8658

25、1.6664 5.0000 0.5000 2.7768 3.2587 4.1634 1.9993 54/477.3 龙格 库塔方法（2）如果使用四阶）如果使用四阶R-K方法方法(公式公式)1111. 1)22(61 . 02351. 1)1 . 0(1133. 1)21 . 0(1025. 1)21 . 0(11432101230422032102201 KKKKyyKyKKyKKyKyKn时时55/477.3 龙格 库塔方法其余结果如下其余结果如下: i xi k1 k2 k3 k4 yi 1.0000 0.1000 1.0000 1.1025 1.1133 1.2351 1.1111 2.

26、0000 0.2000 1.2346 1.3756 1.3921 1.5633 1.2500 3.0000 0.3000 1.5625 1.7639 1.7908 2.0423 1.4286 4.0000 0.4000 2.0408 2.3428 2.3892 2.7805 1.6667 5.0000 0.5000 2.7777 3.2600 3.3476 4.0057 2.000056/4700.050.10.150.20.250.30.350.40.450.511.21.41.61.822.2精确 解 00.050.10.150.20.250.30.350.40.450.511.21.41

27、.61.822.2精确 解 Euler法 00.050.10.150.20.250.30.350.40.450.511.21.41.61.822.2精确 解 Euler法 预测 -校 正 系 统 00.050.10.150.20.250.30.350.40.450.511.21.41.61.822.2精确 解 Euler法 预测 -校 正 系 统 改 进 Eul er法 00.050.10.150.20.250.30.350.40.450.511.21.41.61.822.2精确 解 Euler法 预测 -校 正 系 统 改 进 Eul er法 3阶RK法 00.050.10.150.20.2

28、50.30.350.40.450.511.21.41.61.822.2精确 解 Euler法 预测 -校 正 系 统 改 进 Eul er法 3阶RK法 4阶RK法0.1280.130.1320.1340.1360.1380.140.1420.1440.1460.1481.131.141.151.161.171.181.19精确 解 Euler法 预测 -校 正 系 统 改 进 Eul er法 3阶RK法 4阶RK法0.29650.2970.29750.2980.29850.2990.29950.30.30051.4181.421.4221.4241.4261.428精确 解 改 进 Eul er法 3阶RK法 4阶RK法57/4700.050.10.150.20.250.3