第四章 常微分方程数值解法

1四常微分方程数值解法2常微分方程数值解法引言（常微分方程数值解法概述）显式欧拉法、隐式欧拉法、二步欧拉法局部截断误差与精度改进的欧拉方法龙格-库塔方法收敛性与稳定性简述一阶常微分方程组与高阶常微分方程3引言一阶常微分方程初值问题：00(,)()yfxyyxy微分方程初始条件定理：若f(x,y)在某闭区域R：00||,||(0,0)xxayybab上连续，且在R域内满足李普希兹(Lipschitz)条件，即存在正数L，使得对于R域内的任意两值y1,y2，下列不等式成立：1221|(,)(,)|||fxyfxyLyy则上述初值问题的连续可微的解y(x)存在并且唯一。4引言（续）实际生产与科研中，除少数简单情况能获得初值问题的初等解（用初等函数表示的解）外，绝大多数情况下是求不出初等解的。有些初值问题即便有初等解，也往往由于形式过于复杂而不便处理。实用的方法是在计算机上进行数值求解：即不直接求y(x)的显式解，而是在解所存在的区间上，求得一系列点xn(n0,1,2,…)上解的近似值。5欧拉（Euler）方法方法一化导数为差商的方法10()()()()()limnnnnnhyxhyxyxyxyxhh由于在逐步求解的过程中，y(xn)的准确值无法求解出来，因此用其近似值代替。为避免混淆，以下学习简记：y(xn)：待求函数y(x)在xn处的精确函数值yn：待求函数y(x)在xn处的近似函数值6代入初值问题表达式可得：根据y0可以一步步计算出函数yy(x)在x1,x2,x3x4,…上的近似值y1,y2,y3,y4,…常微分方程数值解是一组离散的函数值数据，它的精确表达式很难求解得到，但可以进行插值计算后用插值函数逼近y(x)10()()()()()limnnnnnhyxhyxyxyxyxhh100(,)0,1,2,()nnnnyyhfxynyyx1(,)nnnnyyfxyh7欧拉方法（续）方法二数值积分法1111d()()[,()]d()[,()][,()]nnnnxxnnxxnnnnnnyxyxyxfxyxxxxfxyxhfxyx同样以近似值yn代替精确值y(xn)可得：100(,)()nnnnyyhfxyyyx将微分方程yf(x,y)在区间[xn,xn+1]上积分：8欧拉方法的几何意义xy00x1x2x3x4x5x6x0P1P2P3P4P5P6P()yyx9隐式欧拉法在数值积分法推导中，积分的近似值取为积分区间宽度与右端点处的函数值乘积，即：这样便得到了隐式欧拉法：11100(,)()nnnnyyhfxyyyx含有未知的函数值隐式欧拉法没有显式欧拉法方便11111111d()()[,()]d()[,()][,()]nnnnxxnnxxnnnnnnyxyxyxfxyxxxxfxyxhfxyx10二步欧拉法在数值积分法推导中，积分区间宽度选为两步步长，即积分区间为：[xn1,xn1]，则：11111111d()()[,()]d()[,()]2[,()]nnnnxxnnxxnnnnnnyxyxyxfxyxxxxfxyxhfxyx以y(x)在xn1,xn上的近似值代替精确值可得：11002(,)()nnnnyyhfxyyyx需要前两步的计算结果中矩形公式11梯形公式欧拉法在数值积分法中，如果用梯形公式近似计算f(x,y)在区间[xn,xn+1]上的积分，即：用近似值代替精确值可得梯形公式欧拉法：111[(,)(,)]2nnnnnnhyyfxyfxy上式右端出现了未知项，可见梯形法是隐式欧拉法的一种；实际上，梯形公式欧拉法是显式欧拉法与隐式欧拉法的算术平均。111111[,()][,()][,()]d()2[,()][,()]2nnxnnnnnnxnnnnfxyxfxyxfxyxxxxhfxyxfxyx12例用显式欧拉法、隐式欧拉法、梯形法求解初值问题：1(0)1yyxy取h0.1，计算到x0.5，并与精确解进行比较解：由已知条件可得：h0.1，x00，y01，f(x,y)yx1显式欧拉法：1(1)(1)0.90.10.1nnnnnnnnyyhyxhyhxhyx13例：（续）隐式欧拉法：111(1)nnnnyyhyx化简得：11(1)(1)(1)nnnnnhyyhxyhxh11[(1)](0.10.11)1.1(1)nnnnnyyhxhyxh111[(1)(1)]2nnnnnnhyyyxyx梯形公式欧拉法：11[1(2)][1(2)](2)(2)nnnnhyhyhxx11[(2)(2)](2)(1.90.20.21)2.1nnnnnnyhyhxxhyx14计算结果：xn显式法yn隐式法yn梯形法yn精确解y(xn)0.011110.11.0000001.0090911.0047621.0048370.21.0100001.0264461.0185941.0197310.31.0290001.0513151.0406331.0408180.41.0561001.0830141.0700971.0703200.51.0904901.1209221.1062781.106531()exyxx本题的精确解为：10.90.10.1nnnyyx1(0.10.11)1.1nnnyyx1(1.90.20.21)2.1nnnyyx显式法隐式法梯形法15局部截断误差为了简化分析某常微分方程数值算法的误差，现假设yny(xn)，即在前一步yn准确的前提下，估计：111()nnnTyxy称上述误差Tn1为该常微分方程数值算法的局部截断误差如果某个常微分方程数值算法的局部截断误差可表示为O(hp1)，则称该数值算法的精度是p阶欧拉法的精度为一阶；二步欧拉法的精度为二阶；梯形公式欧拉法的精度为二阶。16泰勒展开法如果初值问题中的f(x,y)充分可微，则可将y(xn1)在点xn处展开：21()()()()()2!nnnnnhyxyxhyxhyxyx如果只保留线性项，忽略h2及以上各项，则：11()()()nnnnyxyxhyxy()[,()]nnnyxfxyxny1(,)nnnnyyhfxy显式欧拉公式17局部截断误差的分析利用泰勒公式展开，比较各算法与展开式的前几项将y(xn1)在xn点处用泰勒公式展开：211()()()()()2!(,)nnnnnnyyxyxhyxhyxhxx22111()()()2!nnnyTyxyhOh显式欧拉法的局部截断误差：1()[,()]()()nnnnnnyyxhfxyxyxhyx1(,)nnnnyyhfxy欧拉法()nnyyx令：1阶精度18补充：二元函数微分中值定理1100(,)(,)fxyfxy1001010010()[()]()[()]xyxxfxxxyyfyyy010000(,)(,)fxhykfxy00()()xyhfxhkfyk000000(,)(,)()()xyfxhykfxyhfxhkfyk0000(,)(,)fxyhkfxhykxy191111111()(,)[,](,)[()]nnnynnnnyfxyfxxxfxyy111()[()()(,)]nnnynnyyxhyxOhfxT21()()()()nnnyxyxhyxOhy(xn1)在xn点处展开：211111()(,)()nnnynnTyxyhfxTOh111(,)nnnnyyhfxy隐式欧拉法：1阶精度()nyx2111()1(,)nynTOhhfx11[(),]nnyxy111()(,)nynnyxfxT11()()(,)nynnyxOhfxT211()()(,)()nnynnyxhyxhfxTOh20111()2[,()]()2()nnnnnnyyxhfxyxyxhyx分别将y(xn1),y(xn1)在xn点处用泰勒公式展开：231231()()()()()2!()()()()()()2!nnnnnnnnyxyxyxhyxhOhyxyxyxhyxhOh111()nnnTyxy二步欧拉法的局部截断误差：112(,)nnnnyyhfxy二步欧拉法：2阶精度1()nyx()nyx311()()2()()nnnyxyxhyxOh211111111111(,)[,](,)[()]()()(,)nnnynnnnynnnfxyfxfxyyxyxfTyxx111[(,)(,)]2nnnnnnhyyfxyfxy梯形公式欧拉法：()nyx2111()[2()()()(,)]2nnnnynnhyyxyxhyxOhfxTy(xn1)在xn点处展开：231()()()()()2!nnnnhyxyxhyxyxOh3111(,)()2nynnhTfxTOh31211()1(,)nhynTOhfx2阶精度11[(),]nnyxy2311()()()(,)()22nnnynnhhyxhyxyxfxTOh211()()()(,)nnynnyxhyxOhfxT22各种欧拉法的比较方法精度评述显式欧拉法1最简单，精度低隐式欧拉法1不便计算，稳定性好二步欧拉法2需要两步初值，且第2个初值只能由其它方法给出，可能对后面的递推精度有影响梯形公式欧拉法2精度有所提高，但为隐式，需要迭代求解，计算量大23改进的欧拉法从上述例子可以看到，梯形法由于具有二阶精度，其局部截断误差比显式欧拉法和隐式欧拉法小，但梯形法实质上是一种隐式算法显式欧拉法是一个显式算法，虽然计算量较小，但是精度不高综合两种方法的长处，可以先用显式欧拉法求出y(xn+1)的一个粗略近似值，然后用它代入梯形法公式的右端，用梯形法计算y(xn+1)的较为精确的近似值。1ny预报值1ny校正值24改进的欧拉法（续）按照上述思想，可以建立如下预报-校正系统：1(,)nnnnyyhfxy预报:按以上两式求解常微分方程的算法称为改进的欧拉法，它还可以表示为：11(,)[,(,)]2nnnnnnnnhyyfxyfxyhfxy嵌套形式11(,)(,)()2pnnncnnpnpcyyhfxyyyhfxyyyy平均化形式2阶精度111[(,)(,)]2nnnnnnhyyfxyfxy校正:25用改进欧拉法求上例所述的初值问题并与欧拉法和梯形法比较误差的大小。解：采用改进欧拉法的嵌套形式：11(,)[,(,)]2nnnnnnnnhyyfxyfxyhfxy10.1,0.5(0)1yyxhxy计算到1(1)[(,)]12nnnnnnnhyyxyhfxyx(1)[(1)]12nnnnnnnhyyxyhyxxh(2)(2)122nnhhhhyxh0.9050.0950.1nnyx