第五章：常微分方程数值解法第一节欧拉法

1、第五章第五章 常微分方程数值解常微分方程数值解/* Numerical Methods for Ordinary Differential Equations */ 待求解的问题待求解的问题：一阶一阶常微分方程的常微分方程的初值问题初值问题 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy 解的存在唯一性解的存在唯一性（“常微分方程常微分方程”理论）：只要理论）：只要 f (x, y) 在在a, b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条件条件，即存在与，即存在与 x, y 无关的常数无关的常数 L 使使对任意

2、定义在对任意定义在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，则上述都成立，则上述IVP存存在唯一解在唯一解。| ),(),(|2121yyLyxfyxf 解析解法解析解法：（常微分方程理论）：（常微分方程理论）只能求解极少一类常微分方程；实际中给定的问题不一只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法。定是解析表达式，而是函数表，无法用解析解法。如何求解如何求解计算解函数计算解函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值),., 1()(nixyyii 节点间距节点间距 为步长，通

3、常采用为步长，通常采用等距节点等距节点，即取即取 hi = h (常数常数)。) 1,., 0(1 nixxhiii数值解法数值解法: 求解所有的常微分方程求解所有的常微分方程步进式步进式：根据已知的或已求出的节点上的函数值计算：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。因此只需当前节点上的函数值，一步一步向前推进。因此只需建立由已知的或已求出的节点上的函数值求当前节点建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。函数值的递推公式即可。111()()() ()()(,)nnnnnnnnnny xy xhy xy xyy xyyh f xy

4、1(,) 0, 1,.nnnnyyh f xyn-Eulers Method1 欧拉方法欧拉方法 /* Eulers Method */1 Eulers MethodTaylor展开法展开法几何意义几何意义亦称为亦称为欧拉折线法欧拉折线法 /* Eulers polygonal arc method*/ 几何直观是帮助我们寻几何直观是帮助我们寻找解决一个问题的思路找解决一个问题的思路的好办法哦的好办法哦定义定义在假设在假设 yn = y(xn)，即第，即第 n 步计算是精确的前提下，考步计算是精确的前提下，考虑公式或方法本身带来的误差虑公式或方法本身带来的误差: Rn = y(xn+1) yn

5、+1 , 称为称为局部局部截断误差截断误差 /* local truncation error */。说明 显然，这种近似有一定误差，显然，这种近似有一定误差，而且步长越大，误差越大，而且步长越大，误差越大，如何估计这种误差如何估计这种误差y(xn+1) yn+1 ？1 Eulers Method截断误差截断误差: 实际上，实际上，y(xn) yn， yn 也有误差，它对也有误差，它对yn+1的误差的误差也有影响，见下图。但这里不考虑此误差的影响，仅考虑方也有影响，见下图。但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此称为方法误差或截断误差。法或公式本身带来的误差，因此称为方法

6、误差或截断误差。局部截断误差的分析局部截断误差的分析：由于假设：由于假设yn = y(xn) ，即，即yn准确，因此准确，因此分析局部截断误差时将分析局部截断误差时将y(xn+1) 和和 yn+1都用点都用点xn上的信息来表上的信息来表示，工具：示，工具：Taylor展开。展开。 欧拉法的局部截断误差：欧拉法的局部截断误差：223111232() ( )( )( )( ) ( ,) ( )( )hnnnnnnnnnhnRy xyy xhy xy xO hyhf x yy xO hRn+1 的的主项主项/* leading term */1 Eulers Method定义定义若某算法的局部截断误

7、差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。 欧拉法具有欧拉法具有 1 阶精度。阶精度。在在xn点用一阶向前差点用一阶向前差商近似一阶导数商近似一阶导数1()()()nnny xy xyxh 在第在第2章讨论牛顿插值公式时章讨论牛顿插值公式时介绍了介绍了差商差商的概念和性质，的概念和性质，各阶差商可以近似各阶导数，具有不同的精度，各阶差商可以近似各阶导数，具有不同的精度，且可以用函数值来表示。且可以用函数值来表示。上一章中数值微分的方法之一上一章中数值微分的方法之一就是用差商近似导数就是用差商近似导数111()()() ()()(,)nnnnnnn

8、nnny xy xhy xy xyy xyyh f xyEulers method1 Eulers Method1 Eulers Method 欧拉公式的改进欧拉公式的改进：隐式欧拉法或后退欧拉法隐式欧拉法或后退欧拉法 /* implicit Euler method or backward Euler method*/11()()()nnny xy xy xhxn+1点向后差商近似导数点向后差商近似导数111111()()() ()()(,)nnnnnnnnnny xy xhy xy xyy xyyh f xy隐式或后退欧拉公式隐式或后退欧拉公式由于未知数由于未知数 yn+1 同时出现在等式

9、的两边，故称为同时出现在等式的两边，故称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公欧拉公式。隐式公式不能直接求解，一般需要用式。隐式公式不能直接求解，一般需要用Euler显式公式显式公式得到初值，然后用得到初值，然后用Euler隐式公式迭代求解。因此隐式公隐式公式迭代求解。因此隐式公式较显式公式计算复杂，但稳定性好式较显式公式计算复杂，但稳定性好01(1)( )111(,)(,)nnnnkknnnnyyh f xyyyh f xy(1)()1111111()(0 )1111(1)11111()1(,)(,)

10、1, ()(,)kknnnnnnkknnnnknnnnnnknyyh fxyfxyhL yyhLyyhLyykyyh fxyy 在 迭 代 公 式 中 取 极 限 ， 有因 此的 极 限 就 是 隐 式 方 程 的 解收敛性收敛性1 Eulers Method 见上图，见上图， 显然，这种近似也有一定误差，显然，这种近似也有一定误差，如何估计这种误差如何估计这种误差y(xn+1) yn+1 ？方法同上，基于方法同上，基于Taylor展开估计局部截断误差。展开估计局部截断误差。但是注意，隐式公式中右边含有但是注意，隐式公式中右边含有f(xn+1 , yn +1 ) ,由于由于yn +1不准确，所

11、以不能直接用不准确，所以不能直接用y (xn+1)代替代替f(xn+1 , yn +1 ) 设已知曲线上一点设已知曲线上一点 Pn (xn , yn ),过该过该点作弦线，斜率为点作弦线，斜率为(xn+1 , yn +1 ) 点的点的方向场方向场f(x,y)方向方向,若步长若步长h充分小，充分小，可用弦线和垂线可用弦线和垂线x=xn+1的交点近似的交点近似曲线与垂线的交点。曲线与垂线的交点。几何意义xnxn+1PnPn+1xyy(x)1 Eulers Method 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11111111121111111321, ,2 , 2nnnnynnnnnn

12、nnnnnnynnnnnnnnfxyfxy xfxyy xyy xhfxy xyxyxhyxyxyhfxyy xy xhhyxh yxyxy xy x由微分中值定理，得在，之间；又而 2326nnnnhhhyxyxyx1 Eulers Method111111232311(), 231,23nnnynnnnnynnnnRy xyhfxy xyhhyxyxhhhfxRyxyx 从而即2211121,1,(1)1ynynynhfxhfxhfxxxx 111 Eulers Method 2322111231221 1,23 3,226 ()2nynynnnnynnnnnhhRhfxhfxy xyxh

13、hy xfxy xy xhRy xo h 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：111()nnnRy xy232()()hny xO h即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。1 Eulers Method1(,) 0, 1,.nnnnyyh f xyn比较尤拉显式公式和隐式公式及其局部截断误差231112()()()hnnnnRy xyy xO h显式公式111(,)nnnnyyh f xy隐式公式231112()()()hnnnnRy xyy xO h1 Eulers Method 若将这两种方法进行算术平均，即可消除误差若将这两种方法进行算术平均，即可消除

14、误差的主要部分的主要部分/*leading term*/而获得更高的精度而获得更高的精度,称为梯形法称为梯形法1 Eulers Method 梯形公式梯形公式 /* trapezoid formula */ 显、隐式两种算法的显、隐式两种算法的平均平均111 (,)(,)2nnnnnnhyyf xyf xy注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。但阶精度，比欧拉方法有了进步。但注意到该公式是注意到该公式是隐式隐式公式，计算时不得不用到迭代公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。法，其迭代收敛性与欧拉公式相似

15、。3111()()nnnRy xyO h梯形法的迭代计算和收敛性梯形法的迭代计算和收敛性01(1)( )111(,)(,)(,)2nnnnkknnnnnnyyh f xyhyyf xyf xy(1)()1111111()(0 )1111(1)11111()1(,)(,)2222 1, ()2(,)(,)2kknnnnnnkknnnnknnnnnnnnknhyyfxyfxyhhLL yyyyhLhyykLhyyfxyfxyy 当 时 ，在 迭 代 公 式 中 取 极 限 ， 有因 此的 极 限 就 是 隐 式 方 程 的 解收敛性收敛性1 Eulers Method梯形法的简化计算梯形法的简化计

16、算 迭代计算量大，且难以预测迭代次数。为了控制计算量，通常只迭代迭代计算量大，且难以预测迭代次数。为了控制计算量，通常只迭代一次就转入下一点的计算。用显式公式作预测，梯形公式作校正，得到如下一次就转入下一点的计算。用显式公式作预测，梯形公式作校正，得到如下预测校正系统，也称为改进尤拉法预测校正系统，也称为改进尤拉法： 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1nnnnyxfhyy Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 n

17、y),(),(2111 nnnnnnyxfyxfhyy11( ,),( ,)2nnnnnnnnhyyf x yf xyh f x y1 Eulers Method注：注：此法亦称为此法亦称为预测预测-校正法校正法 /* predictor-corrector method */。可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将。后面将看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。1 Eulers Method+1+11+1 (,)(,(,)

18、21 (,)(,)2nnnnnnnnnnnnnnnhyyf xyf xh yh f xyyyhf xyyhf xy或其它形式其它形式1+1(,) (,)12pnnncnnpnpcyyh f xyyyh f xyyyy或几何解释几何解释xnxn+1ABPn+1=(A+B)/2尤拉法尤拉法后退尤拉法后退尤拉法梯形法梯形法1 Eulers Method 0)(,),(yaybaxyxfdxdy 00( ),xxy xyft y tdt令令x=x1,得得 1010(),xxy xyft y tdtAnother point of view对右端积分采用左矩形、右矩形、梯形积分公式，即对右端积分采用左矩

19、形、右矩形、梯形积分公式，即可得尤拉显式、隐式、梯形公式可得尤拉显式、隐式、梯形公式1 Eulers Method 中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假设假设 ，则可以导出，则可以导出即中点公式也具有即中点公式也具有 2 阶精度，且是显式的。阶精度，且是显式的。)(),(11iiiixyyxyy )()(311hOyxyRiii 需要需要2个初值个初值 y0和和 y1来启动递推来启动

20、递推过程，这样的算法称为过程，这样的算法称为双步法双步法 /* double-step method */，而前面的三种算法都是，而前面的三种算法都是单步法单步法 /* single-step method */。1 Eulers Method几何解释几何解释xnxn+1尤拉法尤拉法后退尤拉法后退尤拉法中点法中点法xn-1 0)(,),(yaybaxyxfdxdy 00( ),xxy xyft y tdt令令x=x2,得得 2020(),xxy xyft y tdtAnother point of view对右端积分采用中矩形公式即得中点公式对右端积分采用中矩形公式即得中点公式1 Eulers

21、 Method公式公式局部截断误差局部截断误差精精度度显显隐隐稳稳定定性性步数步数尤拉显尤拉显式公式式公式1 1阶阶显显差差单步单步尤拉隐尤拉隐式公式式公式1 1阶阶隐隐好好单步单步梯形公梯形公式式2 2阶阶隐隐差差单步单步中点法中点法2 2阶阶显显好好二步二步3(3)3nhyx3(3)12nhyx2(2)2nhyx2(2)2nhyxsummary 算例算例 分别用显式分别用显式Euler方法，梯形方法和预估校正方法，梯形方法和预估校正Euler方法方法解初值问题解初值问题解：解：取取 h =0.1，1),( xyyxf, 1 xydxdy10 x1)0( y 梯形方法为：梯形方法为：),()

22、,(2111 kkkkkkkyxfyxfhyy)1() 1(211 kkkkkkxyxyhy).1( 1 . 0) 1(,1 . 0, 1 . 01 khkxkkhxhkkkkk10110521191 kyykk续续 算例算例 分别用分别用显式显式Euler方法，梯形方法和预估校正方法，梯形方法和预估校正Euler方法方法解初值问题解初值问题解：解：取取 h =0.1，1),( xyyxf, 1 xydxdy10 x1)0( y 梯形方法为：梯形方法为：),(),(2111 kkkkkkkyxfyxfhyy)1() 1(211 kkkkkkxyxyhy).1( 1 . 0) 1(,1 . 0,

23、 1 . 01 khkxkkhxhkkkkk10110521191 kyykk1 . 00095. 0905. 01 kyykk预估校正预估校正EulerEuler方法：方法：续续 EulerEuler方法方法 梯形方法梯形方法 预估校正方法预估校正方法0.01.0000000.01.0000000.01.0000000.00.11.0000004.81031.0047627.51051.0050001.61040.21.0100008.71031.0185941.41041.0190252.91040.31.0290001.21021.0406331.91041.0412184.01040.41.0561001.41021.0700962.21041.0708004.81040.51.0904901.61021.1062782.51041.1070765.51040.61.1314411.71021.1485372.71041.1494045.91040.71.1782971.81021.1962952.91041.1972106.21040.81.2304671.91021.2490193.01041.2499756.51040.91.28