Euler法与修正的Euler法局部截断误差Range-Kutta公式

1、1Euler法与修正的法与修正的Euler法法局部截断误差局部截断误差Range-Kutta公式公式常微分方程数值解2一阶常微分方程初值问题一阶常微分方程初值问题: 000)(),(yxyxxyxfdxdy数值方法数值方法取定离散点取定离散点: x0 x1 x2 xN 其中其中, y = y(x) 是未知函数是未知函数,右端函数右端函数 f(x, y )是已知函数是已知函数, 初值初值 y0 是已知数据。是已知数据。求未知函数求未知函数 y(x) 在离散点处的近似值在离散点处的近似值y1, y2, y3, , yN),(1nnnnyxfhyy ),(yxfdxdy Euler法与修正的法与修正

2、的Euler法法3求解常微分方程初值问题的求解常微分方程初值问题的Euler方法方法 取定步长取定步长: h,记记 xn = x0 + nh, ( n = 1,2, , N )称计算格式称计算格式: yn+1 = yn + h f( xn, yn ) 为为Euler公式公式。对应的求初值问题数值解的方法称为对应的求初值问题数值解的方法称为Euler方法方法。例例2 用用Euler法求初值问题法求初值问题 1)0(20,2yxxyydxdy的数值解。的数值解。解解: 记记 f (x, y) = y x y2, xn= nh (n = 0, 1, 2, N ) 由由Euler公式得公式得: yn+

3、1 = yn + h( yn xn yn2 ) (n = 0, 1, ,N)4取步长取步长 h = 2/10, 2/20, 2/30, 2/40, 用用Euler法求解法求解的数值实验结果如下的数值实验结果如下. N 10 20 30 40 h 0.2 0.1 0.0667 0.05误差误差 0.1059 0.0521 0.0342 0.0256xexxy 211)(解析解解析解:-101234500.511.52 o 数值解数值解- 准确解准确解 Comparison with exact resultsFigure 4. Comparison of Eulers method with e

4、xact solution for different step sizes 56y = f (x, y) 1)(,()()(1nnxxnndxxyxfxyxy),(),(2)(,(111 nnnnxxyxfyxfhdxxyxfnn梯形公式梯形公式: ),(1nnnnyxhfyy 左矩形公式左矩形公式),()(,(1nnxxyxhfdxxyxfnn 用数值积分方法离散化常微分方程用数值积分方法离散化常微分方程 11)(,()(nnnnxxxxdxxyxfdxxy),(),(2111 nnnnnnyxfyxfhyy7),(nnnnyxhfyy1),(),(2111 nnnnnnyxfyxfhyy

5、),(),(2111 nnnnnnyxfyxfhyy预预- -校方法又称为修正的校方法又称为修正的Euler法法,算法如下算法如下 k1 = f(xn , yn) , k2 = f( xn+1 , yn+ h k1) ,2211kkhyynn 由梯形公式推出的预由梯形公式推出的预- -校方法校方法:8-101234500.511.5预预- -校方法校方法, h=0.2时时误差最大值误差最大值: 0.0123-101234500.511.52 n 10 20 30 40 h 0.2 0.1 0.0667 0.05误差误差2 0.0123 0.0026 0.0011 5.9612e-004误差误差

6、1 0.1059 0.0521 0.0342 0.0256欧拉方法欧拉方法, h=0.2时时误差最大值误差最大值: 0.10599设设 yn= y(xn), 称称 Rn+1=y(xn+1) - yn+1为为局部截断误差局部截断误差.)(2)()()()()(2111 yxxxyxxxyxynnnnnnn )(2),()()(21 yhyxhfxyxynnnn 即即由泰勒公式由泰勒公式Euler公式公式: yn+1 = yn+ hf (xn, yn) 的局部截断误差的局部截断误差y(xn+1) yn+1=y(xn) yn+ O(h2) = O(h2)Euler公式的局部截断误差记为公式的局部截断

7、误差记为: O(h2) 称称Euler公式具有公式具有1阶精度。阶精度。局部截断误差局部截断误差10 若局部截断误差为若局部截断误差为: O(h p +1)则称显式单步法具有则称显式单步法具有 p 阶精度阶精度 。 例例 3 证明修正的证明修正的Euler法具有法具有2阶精度阶精度),(),(2111 nnnnnnyxfyxfhyy),(1nnnnyxhfyy 将预测公式将预测公式代入代入得得 yn+1 = yn + 0.5hf(xn,yn)+f(xn+1, yn+hf(xn, yn)11yn+1 = yn + 0.5hf(xn,yn)+f(xn+1, yn+hf(xn, yn) f(xn+1

8、, yn+hf(xn, yn)= f(xn+h, yn+hf(xn, yn) = f (xn, yn)+hfxn + hf(xn,yn) fyn + O(h2),()(nnnyxfxy ),()(nnnyxfdxdxy 0.5hf(xn,yn)+f(xn+1, yn+hf(xn, yn) =hy(xn)+0.5h2y”(xn)+0.5h2y(xn) fyn+ O(h3)yn+1 = yn + hy(xn)+0.5h2y”(xn)+ O(h3)y = f (x, y)局部截断误差局部截断误差: y(xn+1) yn+1 = y(xn) yn =O(h3)故修正的故修正的Euler法法具有具有2阶

9、精度。阶精度。 000)(),(yxyxxyxfdxdy12三阶三阶Range-Kutta公式一般形式公式一般形式yn+1= yn+ hk1+4k2+k3/6k1=f(xn,yn), k2=f(xn+0.5h, yn+0.5hk1)k3=f(xn+h, yn hk1+2hk2)四阶四阶Range-Kutta公式一般形式公式一般形式yn+1= yn+ hk1+2k2+2k3+k4/6k1=f(xn,yn), k2=f(xn+0.5h, yn+0.5hk1)k3=f(xn+0.5h, yn+0.5hk2), k4=f(xn+h, yn+hk3)Range-Kutta公式公式13 1)0(20,2y

10、xxyydxdyxexxy 211)(例例4数值实验数值实验:几种不同求数值解公式的误差比较几种不同求数值解公式的误差比较 n 10 20 30 40 h 0.2 0.1 0.0667 0.05RK4 6.862e-005 3.747e-006 7.071e-007 2.186e-007RK3 0.0012 1.529e-004 4.517e-005 1.906e-005RK2 0.0123 0.0026 0.0011 5.9612e-004Euler 0.1059 0.0521 0.0342 0.025614差分格式的定性分析差分格式的定性分析 差分算法的差分性质精度精度：差分方程对源方程的

11、逼近误差：差分方程对源方程的逼近误差相容性相容性：时空步长趋于：时空步长趋于0时，差分方程的极限为源方程时，差分方程的极限为源方程相容性、稳定性、收敛性、耗散性、色散性、相容性、稳定性、收敛性、耗散性、色散性、和守恒性和守恒性等等15稳定性稳定性：任何初值扰动对差分数值解的影响随时间推移 不再增加（强稳定）或在一段时间内有界（弱稳定）收敛性收敛性：当步长趋于0，差分数值解收敛于源问题的真解 Lax等价定理等价定理：对一个适定的初值问题，在满足相容条件的前提下，稳定性是收敛性的充要条件。耗散性耗散性(diffusion):差分余项对解产生的耗散效应色散性色散性(dispersive)：差分余项对

12、解产生的色散效应守恒性守恒性：数值解保持真解所固有的守恒性的程度16收敛性收敛性 数值解法的基本思想是，通过某种离散化手段将数值解法的基本思想是，通过某种离散化手段将微分方程转化为差分方程，如单步法，即微分方程转化为差分方程，如单步法，即 ),(1hyxhyynnnn它在它在 处的解为处的解为 ，而初值问题在，而初值问题在 的精确解为的精确解为,记记 称为整体截断称为整体截断 误差误差. . nxnynx)(nxynnnyxye)( 收敛性就是讨论当收敛性就是讨论当 固定且固定且 时时 的问题的问题. . nxx 00nxxhn0ne17 定义定义 若一种数值方法对于固定的若一种数值方法对于固

13、定的 , ,当当 时有时有 ，其中，其中 是原问题的确解，是原问题的确解，则称该方法是则称该方法是收敛收敛的的. . nhxxn00h)(nnxyy)(xy 定理定理 假设单步法具有假设单步法具有p阶精度，且增量函数阶精度，且增量函数 关于关于y 满足利普希茨条件满足利普希茨条件 ),(hyx,),(),(yyLhyxhyx又设初值又设初值 是准确的，即是准确的，即 ，则其，则其整体截断误差整体截断误差 0y)(00 xyy).()(pnnhOyxy18绝对稳定性与绝对稳定域绝对稳定性与绝对稳定域 定义定义 若一种数值方法在节点值若一种数值方法在节点值 上大小为上大小为 的扰动，的扰动，于以后

14、各节点值于以后各节点值 上产生的偏差均不超过上产生的偏差均不超过 ，则称，则称该方法是该方法是稳定稳定的的. . ny)(nmym 以欧拉法为例考察计算稳定性以欧拉法为例考察计算稳定性. . 例例 考察初值问题考察初值问题 .1)0(,100yyy其准确解其准确解 是一个按指数曲线衰减得很快的函数，是一个按指数曲线衰减得很快的函数，如图如图9-39-3所示所示. . xexy100)( 用欧拉法解方程用欧拉法解方程 得得 yy10019.)1001(1nnyhy若取若取 ，则欧拉公式的，则欧拉公式的具体形式为具体形式为 025.0h,5.11nnyy 可以看到，欧拉方法的解可以看到，欧拉方法的

15、解 （图中用（图中用号标出）在准确值号标出）在准确值 的上下波动，计算过程明的上下波动，计算过程明显地不稳定显地不稳定. . ny)(nxy 但若取但若取 则计算过程稳定则计算过程稳定. . nnyyh5.0,005.0120再考察后退的欧拉方法，取再考察后退的欧拉方法，取 时计算公式为时计算公式为 025.0h.5.311nnyy计算结果如下，这时计算过程是稳定的计算结果如下，这时计算过程是稳定的. . 0067.00625.5100.00233.0375.3075.00816.025.2050.02857.05.1025.0 后退欧拉方法欧拉方法节点计算结果对比21其中其中 为为 的近似，

16、的近似， 为常数为常数, , 及及 不全为零，则称为不全为零，则称为线性线性 步法步法. . iny)(inxy,),(0ihxxyxffininininii,00k 线性多步法的一般公式线性多步法的一般公式 如果计算如果计算 时，除用时，除用 的值，还用到的值，还用到 的值，则称此方法为的值，则称此方法为线性多步法线性多步法. . kny1 knyiyin()2, 1 ,0k 一般的线性多步法公式可表示为一般的线性多步法公式可表示为 ,010kiinikiiniknfhyy 计算时需先给出前面计算时需先给出前面 个近似值个近似值 , , 再由再由逐次求出逐次求出 . . k110,kyyy,1kkyy22 阿当姆斯显式与隐式公式阿当姆斯显式与隐式公式 考虑形如考虑形如 kiiniknknfhyy01的的 步法，称为步法，称为阿当姆斯阿当姆斯（Adams)Adams)方法方法. . k 为显式方法，为显式方法， 为隐式方法，通常称为阿为隐式方法，通常称为阿当姆斯显式与隐式公式，也称当姆斯显式与隐式公式，也称Adams-BashforthAdams-Bashforth公式与公式与AdamsAdams-Monlton-Monlton公式公式.