1、 考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题/*Initial-Value Problem*/:00(,)()yfx yaxby xy 其中其中f(x,y)为为x,y的已知函数，的已知函数，y0为给定的初始值为给定的初始值 第六章第六章 常微分方程数值解常微分方程数值解/*Numerical Methods for Ordinary Differential Equations*/例如：例如：,()0101y=x+y xy其解析解为：其解析解为：1 2xyxe 只有一些特殊类型的微分方程问题能够得到用解只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方
2、程问题很难析表达式表示的函数解，而大量的微分方程问题很难得到其解析解。得到其解析解。这是微积分的发明者之一Leibniz在1686年曾经让当时数学界人士求解的一阶微分方程式，吸引了许多数学家的注意，大约经过150年的探索到1838年,刘维尔（Liouville)在理论上证明了这个微分方程不能用初等积分法求解，得借助于数值方法dyxydx22 只能依赖于只能依赖于数值方法数值方法去获得微分方程的去获得微分方程的数值解数值解。bxxxxan210iiixxh1nabh数值方法的基本思想是：在解的存在区间上取数值方法的基本思想是：在解的存在区间上取n+1个节点个节点这里这里hi可以不相等，但一般取成
3、相等的，这时可以不相等，但一般取成相等的，这时在这些节点上采用离散化方法，（通常用数值积分、微分、在这些节点上采用离散化方法，（通常用数值积分、微分、泰勒展开等）将上述初值问题化成关于离散变量的相应问泰勒展开等）将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解题。把这个相应问题的解yn作为作为y(xn)的近似值。这样求得的近似值。这样求得的的yn就是上述初值问题在节点就是上述初值问题在节点xn上的数值解。一般说来，上的数值解。一般说来，不同的离散化导致不同的方法。不同的离散化导致不同的方法。，i=0,1,n-1称为由称为由xi到到xi+1的步长。的步长。步进式步进式：根据已知的或已
4、求出的节点上的函数值计算当前节点根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。上的函数值，一步一步向前推进。因此只需建立由已知的或已求出的节点上的函数值求因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。当前节点函数值的递推公式即可。1 欧拉方法欧拉方法/*Eulers Method*/欧拉公式：欧拉公式：x0 x1向前差商近似导数向前差商近似导数hxyxyxy)()()(010 ),()()()(000001yxfhyxyhxyxy 1y记为记为)1,.,0(),(1 niyxfhyyiiii亦称为亦称为欧拉折线法欧拉折线法/*Eule
5、rs polygonal arc method*/也称欧拉折线法欧拉折线法.用这条折线近似地代替曲线()y x()yy xxynx1nxny1ny1ny x由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙。非常粗糙。例例 用欧拉法求初值问题当h=0.02时在区间0,0.1上的数值解。000.912()10yyxy xx 解解：把 代入欧拉法计算公式，得yxyxf219.0),(5,1,021018.01219.01nyxyxhyynnnnnn具体计算结果nxnyny(xn)n=y(xn)-yn001.00001.0000010.020.98200.98250.000520.040
6、.96500.96600.000530.060.94890.95030.001440.080.93360.93540.001850.100.91920.9230.0021 表中y(xn)，是初值问题的真解 在xn上的值。为近似值yn的误差。从表中可以看出，随着n的增大，误差也在增大，所以说，欧拉法计算简便，对一些问题有较大的使用价值，但是，它的误差较大，所得的数值解精确度不高。45.0)21()(xxyn定义定义在假设在假设 yi=y(xi)，即第，即第 i 步计算是精确的前提下，考步计算是精确的前提下，考虑的截断误差虑的截断误差 Ri=y(xi+1)yi+1 称为称为局部截断误差局部截断误差
7、/*local truncation error*/。定义定义若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。欧拉法的局部截断误差：欧拉法的局部截断误差：),()()()()()(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR )()(322hOxyih 欧拉法具有欧拉法具有 1 阶精度。阶精度。Ri 的的主项主项/*leading term*/衡量求解公式好坏的一个主要标准是求解公式的精度衡量求解公式好坏的一个主要标准是求解公式的精度,因此引入局部截断误差和阶数的概念。因此引入局部截断误差和阶数的概念。Eule
8、rs Method 欧拉公式的改进：欧拉公式的改进：隐式欧拉法隐式欧拉法/*implicit Euler method*/向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,.,0(),(111 niyxfhyyiiii Eulers Method由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐式隐式/*implicit*/欧拉公式，欧拉公式，而前者称为而前者称为显式显式/*explicit*/欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算
9、一个初值，再迭代迭代求解。求解。隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11)(iiiyxyR)()(322hOxyih 即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。Hey!Isnt the leading term of the local truncation error of Eulers method?Seems that we can make a good use of it)(22ihxy Eulers Method 梯形公式梯形公式/*trapezoid formula*/显、隐式两种算法的显、隐式两种算法的平均平均)1,.,0(),(),(2111
10、niyxfyxfhyyiiiiii注：注：局部截断误差局部截断误差 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。)()(311hOyxyRiii例例 1)0(2yyxydxdy在区间在区间0,1.5上上，取取h=0.1，求解。求解。本题的精确解为，本题的精确解为，可用来检验近似解可用来检验近似解xxy21)(的精确程度。计算结果如下表：的精确程度。计算结果如下表：xn 欧拉法欧拉法yn迭代一次迭代
11、一次梯形公式梯形公式yn准确解准确解01110.11.11.0959091.0954450.21.1918181.1840961.1832160.31.2774381.2602011.2649110.41.3582131.3433601.3416410.51.4351331.4161021.4142140.61.5089661.4829561.4832400.71.5803381.5525151.5491930.81.6497831.6164761.6124520.91.7177791.6781681.6733201.01.7847701.7378691.7320511.11.851181.7
12、958221.7888541.21.9174641.8522421.8439091.31.9840461.9073231.8973671.42.0514041.9612531.9493591.52.1200522.0142072.000000 中点欧拉公式中点欧拉公式/*midpoint formula*/中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,.,1),(211 niyxfhyyiiii假设假设 ，则可以导出，则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。)(),(11iiiixyyxy
13、y )()(311hOyxyRiii 需要需要2个初值个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双步法/*double-step method*/，而前面的三种算法都是，而前面的三种算法都是单步法单步法/*single-step method*/。方方 法法 Eulers Method显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低,计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高,显式显式多一个初值多一个初值,可能影响精度可能影响精度 Cant you give
14、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:先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy Step 2:再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxf
15、yxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法/*predictor-corrector method*/。可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将。后面将看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。)1,.,0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii Eulers Method2 龙格龙格-库塔法库塔法/*Runge-Kutta Method*/建立高精度的单步递推格式。建立高精度的
16、单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从(xi,yi)点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？2 Runge-Kutta Method首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有，使
17、得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hOyxyRiii Step 1:将将 K2 在在(xi,yi)点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2:将将 K2
18、代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii 2 Runge-Kutta MethodStep 3:将将 yi+1 与与 y(xi+1)在在 xi 点的点的泰勒泰勒展开作比较展开作比较)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii21,1221 p 这里有这里有 个未知个未知数，数，个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为
19、。所有满足上式的格式统称为2阶龙格阶龙格-库库塔格式塔格式。21,121 p注意到，注意到，就是改进的欧拉法。就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？其中其中 i (i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待定均为待定系数，确定这些系数的系数，确定这些系数的步骤与前面相似。步骤与前面相似。2 Runge-Kutta Method).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法/*Classical Runge-Kutta Method*/：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii
