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

第五章常微分方程的数值解法主要内容：1、引言2、欧拉方法3、龙格－库塔方法4、单步法的收敛性和稳定性5、线性多步法6、一阶方程组与高阶方程第一节引言●在常微分方程课程里面讨论的是一些典型方程求解解析解的基本方法。●然而在生产实践和科学研究中遇到的微分方程往往比较复杂，在很多情况下，不能给出解的解析表达式；有时候即时能用解析表达式来表示，又因为计算量太大而不实用，有时候一些是已经有了求解的基本方法的典型方程，但实际使用时也是有困难的。●以上情况说明用求解解析解的基本方法来求微分方程的解往往是不适宜的，甚至很难办到。●实际问题中，对于求解微分方程，一般只要求得到解的若干个点上的近似值或者解的便于计算的近似表达式。●本章研究微分方程的数值解法，而且着重讨论微分方程中最简单的一类问题——一阶方程的初值问题。第一节引言1、一阶方程的初值问题假定上式在区间[a,b]上存在唯一且足够光滑的解y(x)。●所谓数值解法就是寻求解y(x)在一系列离散点，也称为节点处的值：要计算出解函数y(x)在一系列节点a=x0x1…xn=b处的近似值0)(],[),(yaybaxyxfdxdyy),...,1()(nixyyii第一节引言●节点间距，即步长为：通常采用等距节点，即hi=h(常数)●等间距节点●在这些节点上采用离散化方法（通常用数值积分、微分、泰勒展开等）将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解yn作为y(xn)的近似值。这样求得的yn就是上述初值问题在节点xn上的数值解。一般说来，不同的离散化导致不同的方法。,...)2,1,0(nnhaxn记,...)2,1,0()(,)(nyxyyxxynnnn上的近似值在点)1,...,0(1nixxhiii第二节欧拉方法一、欧拉法Euler1、向前差商近似导数hxyxyxy)()()(010),()(0001yxfhyxy1y记为),(1112yxfhyy)1,...,0(),(1niyxfhyyiiii0)(],[),(yaybaxyxfdxdyy第二节欧拉方法2、举例例1用欧拉法求初值问题当h=0.02时在区间[0,0.10]上的数值解。解：根据欧拉公式可以得到：此外，可以得到方程的真解：01)(219.0'00xxyyxy1)219.0(*02.0),(01yyxyyxhfyyiiiiiii45.0)21()(xxy第二节欧拉方法求解过程如下：nxnyny(xn)n=y(xn)-yn001.00001.0000010.020.98200.98250.000520.040.96500.96600.000530.060.94890.95030.001440.080.93360.93540.001850.100.91920.9230.0021第二节欧拉方法3、欧拉方法的几何意义()yyxxy0x1x2x0P1P2P根据已知条件：曲线y(x)上的点（x0,y0）及该点处曲线的导数f(x0,y0)，则可以得到过该点的直线：))(,(0000xxyxfyy该直线与x＝x1的交点P1，则P1的纵坐标y1为：))(,(010001xxyxfyy),(000yxhfy就用y1作为y(x1)的近似值…逐次进行后可以得到一条折线p0P1…Pn，该折线看作是初值问题的积分曲线的近似，因此欧拉方法也称为欧拉折线法第二节欧拉方法从上述几何意义上得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙即误差太大。4、欧拉法的局部截断误差（1）截断误差定义在假设yi=y(xi)，即第i步计算是精确的前提下，考虑的截断误差Ri+1=y(xi+1)yi+1，称为局部截断误差如图所示：APi+1即为欧拉方法在xi+1点的截断误差（2）如果某种方法的局部截断误差是则称该方法具有p阶精度ix1ix1iPA)(1phO第二节欧拉方法（3）则截断误差的大小？写出y(xn+1)的泰勒展开式：由欧拉方法可以得到：则上面两个公式相减得到：)(!2)(')()()(21nnnnnxyhxhyxyhxyxy)(')(),(1nnnnnnxhyxyyxhfyy)(2)(2111nnnnxyhyxyR)()(222hoxyhn具有1阶精度第二节欧拉方法二、改进的欧拉法一阶方程的初值问题与如下积分方程是等价的：当x=x1时可以借助于数值积分，求y(x1)的值1、用矩形公式xxdttytfyxy0))(,()(010))(,()(01xxdttytfyxy10)())(,()(,(0100xxxxxyxfdttytf第二节欧拉方法可以推导出：用矩形法计算右端的积分与用欧拉法计出的结果完全相同2、用梯形公式则可以推导出：)))((,()(010001xxxyxfyxy),(000yxhfy10)())}(,())(,({21))(,(011100xxxxxyxfxyxfdttytf)],(),([21110001yxfyxfhyy),(),(2111nnnnnnyxfyxfhyy第二节欧拉方法梯形公式的截断误差：111111()()[(,)(,)]2iiiiiiiiihyxyfxyTyyxyfx231()()()[()()]23!2iiiiihhhhyxyxyxyxyx2324()()()23![()()()()]()22iiiiiiihhhyxyxyxhhyxyxhyxyxOh梯形公式具有二阶精度，比欧拉方法有了进步)()(1233hOxyhi第二节欧拉方法和欧拉公式相比较，梯形公式在计算yi+1时候也只用到前一步的值yi，但是若yi已知，将yi带入公式求解时候，一般不能直接得到yi+1，而需要通过其他方法（比如迭代法）求解，所以梯形公式被称为隐式公式。3、改进的欧拉方法梯形公式是隐式的，一般用迭代法求解，计算量较大。实际中常将欧拉公式和梯形公式联合使用，先用欧拉公式得出一个y(xi+1)的近似值称为预估值，然后对预估值使用梯形公式对它进行精确化，得到较为精确的近似值yi+1，称之为校正值，计算公式为：这样的预估校正系统称为改进的欧拉方法。1iy,...2,1,0),(),(2),(1111iyxfyxfhyyyxhfyyiiiiiiiiii第二节欧拉方法为了便于编写程序，常将上面的公式改写为如下式：)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy第二节欧拉方法4、举例在区间[0,1.5]上，取h=0.1，求解。解：（1）用欧拉法计算公式如下：（2）用改进欧拉法计算公式如下：1)0(2yyxydxdyy1.0,1201hyyxyhyynnnnnnnnnnyxyhyy21111111222),(),(2nnnnnnnnnnnnnyxyyxyhyyxfyxfhyy1.0,10hy本题的精确解为，可用来检验近似解的精确程度。计算结果如表：xxy21)(nnxxy21)(xn欧拉法yn迭代一次改进欧拉法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.7958221.7888541.21.9174641.8522421.8439091.31.9840461.9073231.8973671.42.0514041.9612531.9493591.52.1200522.0142072.000000第二节欧拉方法P93例5-2第二节欧拉方法5、欧拉两步公式中心差商00110011)(),(2)(2)()(),(yxyyxhfyyyxyhxyxyyxfynnnnnn第三节龙格－库塔方法一、龙格库塔法的基本思想1、平均斜率考察差商：根据微分中值定理：在闭区间[a,b]上连续，开区间(a,b)上可导，则至少存在(a,b)上的一点，使得下式成立：根据上面公式可以得到：hxyxyii)()(1abafbff)()()(),(ba)](,[)()(1hxyhxhfxyxyiiii称为区间[xi,xi+1]上的平均斜率*K)1,0(第三节龙格－库塔方法因此只要对K*提供一种算法，就可以求得数值解，根据该观点对欧拉法及改进的欧拉进行分析。2、基于平均斜率对欧拉法和改进的欧拉法进行分析（1）在欧拉公式中，是取了一个点xi上的斜率值f(xi,yi)作为平均斜率K*的近似值的，已经知道其精度较低。（2）对于改进的欧拉公式：)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy),(),()(21121211hKxfKyxfKKKhyyiiiii可以看出它用两个点xi和xi+1上的斜率K1和K2的算术平均值作为平均斜率K*的近似值的，而xi+1处的斜率是由已知信息预测得到的。第三节龙格－库塔方法根据上面的分析得到：如果在区间[xi,xi+1]上多测几个点的斜率值，然后取其加权平均作为平均斜率的近似值，有可能构造出具有更高精度的计算公式，此即为龙格库塔算法的基本思想。二、二阶龙格库塔方法1、推广改进的欧拉方法，考察区间[xi,xi+1]上的一点：用xi和xp两个点上的斜率值K1和K2的加权平均作为平均斜率K*的近似值：即取：与改进的欧拉法类似，有：如何得到xp的斜率？phxxipi10p2211*KKK)(22111KKhyyii为待定常数21,),(1iiyxfK第三节龙格－库塔方法如何得到xp的斜率？根据改进的欧拉法，可以利用欧拉法预测的值：则可以得到点xi+p斜率K2：则可以得到算法的具体表达式：)(ypix1yphKyipi),(2pipiyxfK),(),()(12122111phKyxfKyxfKKKhyyipiiiii为参数p,,21第三节龙格－库塔方法2、选择参数使得算法具有2阶精度计算上面公式的局部截断误差：根据泰勒公式有：另外：)(O)(!2)(')()(3i2ii1ihxyhxhyxyxy),(iiiyxhfy)(O)],(),(),([!232hyxfyxfyxfhiiiiyiix))],(,(),([)(2122111iiipiiiiiiyxphfyxfyxfhyKKhyy第三节龙格－库塔方法考虑到：则有：),()),(,(iiiiipiyxfyxphfyxf)(O)],(),(),([2hyxfyxfyxfphiiiiyiix二元泰勒展开)]},(),(),(),([),({21iiiiyiixiiiiiyxfyxfphyxfphyxfyxfhy