应用数值分析常微分方程求解

1、第九章 常微分方程数值解9.4 线性多步法线性多步法9.5 收敛性与稳定性分析收敛性与稳定性分析9.4 线性多步法线性多步法在前面所讨论的方法中 ,在计算yn ?1时只用到前一步yn的信息(单步法)，为提高截断误差的阶，每个时间步必须增加计算右端函数f(x, y)的次数。当f(x, y)的结构比较复杂时，计算量较大。现在指出另一个提高截断误差阶的办法，即构造这样的方法 : 在计算yn?1公式中,充分利用前几步得到的信息yn, yn?1,?及fn? f(xn, yn),fn?1? f(xn?1, yn?1)，但每进一步，只计算一次f(x, y)的值。这样的方法称为多步方法，若函数值fn, fn?

2、1,?以线性组合的形式出现于公式中，则称方法为线性多步方法。初值问题:?y? f(x,y)a? x? b?y(a) ? y0中心差商 y(xf(xn? h)? f(xn? h)n) ?2hyn?1? yn?1? 2hf(xn, yn)称为Euler二步法.计算yn?1时,用到了yn, yn?1.-(1)离散方程得 1. (l步)线性多步法的一般格式为：yn?1?0yn?1yn?1?L ?l?1yn?(l?1)? h(?1fn?1?0fn?L ?1fn?(l?1)-(2)fi? f(xi,yi)当当?1?0 时，为隐式公时，为隐式公; ?1=0 则为显式公式。则为显式公式。其中?0,L,?l?1

3、,?1,?0,L,?l?1是独立常数, ?l?1,?l?1不全为零.局部截断误差:Tn?1? y(xn?1) ?iy(xn?i) ? h?if(xn?i, y(xn?i)i?0i?1l?1l?1? y(xn?1) ?iy(xn?i)? h?iy(xn?i)i?0i?1l?1l?1? y(xn?1) ?iy(xn?ih) ? h?iy(xn?ih)i?0i?1l?1l?1-(3)2. 线性多步方法的构造构造多步法有多种途径，常用的有基于 Taylor展开的构造方法(待定系数法)和基于数值积分的构造方法。根据公式具体的精度 (p阶),即Tn?1? O(h),确定p?1?i,?i(2l ?1个未知数

4、)将(3)中所有项全在x? xn进行Taylor 展开,观察含h幂次前的系数,令h 的低于p?1的幂次项系数为零.3. 几个重要的线性多步法Adams(阿当姆斯)方法, Milne(米纳)方法, Hamming(哈明)方法, Simpson(辛普生)方法.通过指定l和p确定?i,?i来构造.(1) Adams 4步4阶显式公式hyn?1? yn?(55fn?59 fn?1? 37 fn?2? 9 fn?3)242515(5)6Tn?1?h y (xn)?O(h )270(2) Adams 3步4阶隐式公式yyhn?1?n?24(9 fn?1?19 fn? 5fn?1? fn?2)T?195(5

5、)6n?1?270h y (xn)? O(h )Milne4步4阶显式公式yy4hn?1?n?3?3(2 fn? fn?1? 2 fn?2)T145(5)6n?1?45h y (xn) ?O(h )(3)(4) Ham ming 3步4阶隐式公式y?18(9y3hn?1n?yn?2)?8( fn?1? 2 fn? fn?1)T15(5)6n?1? ?40h y (xn)?O(h )Simpson 2步4阶隐式公式yhn?1? yn?1?3( fn?1? 4 fn? 2 fn?1)T? ?15(5)6n?190h y (xn)?O(h )(5)多步方法的特点：多步方法的特点：(1)、 因初始条件

6、只有一个，多步方法的启动要借助因初始条件只有一个，多步方法的启动要借助高阶的单步方法来开始高阶的单步方法来开始 . (2)、多步方法比较简单，只要在这几个点的函数、多步方法比较简单，只要在这几个点的函数值的线性组合值的线性组合, 而且每步中所用函数值而且每步中所用函数值 , 有些下一有些下一步还可使用。步还可使用。9.5 收敛性与稳定性分析收敛性与稳定性分析1. 单步法(显式)的收敛性yn?1? yn? h?(xn, yn,h)yn?1? yn-(1)?(xn, yn,h)h用上式的差分方程来逼近微分方程的初值问题是否合理,就要看差分方程的解是否收敛到初值问题的精确解.?定义定义1由初值问题的

7、显式单步法 (1)产生的近似解 yn,如果对任一固定的 xk? x0? kh均有 yk?y(xk) as h? 0.则称单步法是收敛的.(即整体截断误差 Ek? y(xk) ? yk趋于零)由定理可得:1. Euler公式,改进的Euler公式和Runge ? Kutta公式均收敛.2. 一个方法的整体截断误差比局部截断误差低一阶 .若某些引入的误差,在以后的传播中被压缩,衰减或增长可以控制,就认为数值方法 (1)是数值稳定的,反之,若在传播中被放大而无法控制,就认为是数值不稳定.其中,若误差的传播可以被压缩,衰减, 则称绝对稳定.?y=f(x,y),x? D定义2 对初值问题?对于固定的?y

8、(x0) ? y0,步长 h,在数值计算中, 节点值 yi产生一扰动?i(包括初值y0),而仅由这一个扰动引起的以后各节点值 yj( j ? i)的变化?j都不超过?i,即|?j|?|?i|, 就称这个数值方法是稳定的.讨论具体初值问题的数值方法(1) 的数值绝对稳定性往往比较困难,原因在于f(x, y)常常比较复杂的函数.在具体的实际应用中, 通常用试验方程 y ?y 来检验,其中?为复常数(称模型方程) .选择这一试验方程的理由 :(1)该试验方程比较简单 ,若对它方法已不稳定,则对其他方程也就靠不住了 .(2) 对于一般的初值问题 y ? f(x, y)可局部化成这一形式.y ? f(x

9、, y)? f(x0, y0)? fx(x0, y0)(x? x0) ? fy(x0, y0)(y? y0)?O(| x? x0|? | y? y0|)? f(x0, y0)? fx(x0, y0)(x? x0)? y0fy(x0, y0)? yfy(x0, y0)?O(| x? x0|? | y? y0|)略去O(| x? x0| ? | y? y0|)项,并%? y(x) ?令 ? fy(x0, y0), yfx(x0, y0)%y ?y? f(x0, y0) ? x0fx(x0, y0) ? y0fy(x0, y0) ?x,则fx(x0, y0)?在平移一下,即化成检验方程形式.?y ?

10、y?y(x0) ? y0-(2)?y? ye0?(x?x0)(y0? 0)?当Re? 0时,?当Re? 0时,其关系式为y(x)|? ? (asx? ? );y(x) |? 0 (as x? ? ),此时,试验方程是稳定的.一般地, 设单步法(1)用于试验方程(2), 从 yn计算一步得 yn?1,yn?1? E(?h)yn.其中记号E(?h)与所选数值方法有关.若在节点值yi产生一扰动?i,而以后的计算全是正确的,则在计算yj时将有误差 E(?h)?i.j?i定义定义就是绝对稳定的.(1) 当|E(?h)|?1时,数值方法(1) 在复平面上,变量h ?h,满足|E(h)|?1的区域称为方程的

11、绝对稳定区域,绝对稳定域与实轴之交称为绝对稳定区间.(2) 当方法A的绝对稳定域比B的绝对稳定域大时,称方法A比方法B稳定.1. Euler公式yn?1? yn? hf(xn, yn) ? yn?hyn? (1?h)yn要使 |1 ?h|?1,即 |1? h|?1给出了绝对稳定区域z | z?1|?1|,这是复平面上以 (?1,0)为圆心的单位圆,绝对稳定区间为(-2,0).2. 隐式Euler公式yn?1? yn? hf(xn?1, yn?1) ? yn?hyn?1?1yn?1?yn1?h11E(h) ?,1?h1? h1|?1? |1? h |?11? h稳定区域 :Re(h) ? Re(

12、?h) ? 0,绝对稳定区间(? ,0所以,隐式Euler方法比显示Euler方法的稳定性好.3. 梯形方法hhyn?1? yn?f(xn, yn)? f(xn?1, yn?1)? yn?(yn? yn?1),22?h1?1?h2yn?1?yn?yn?h1?h1?2? h? 0, Re(h) ? 0,1? h| E(h) |?|? 1成立.1? h绝对稳定区域 : 左半平面;绝对稳定区间 :(? ,04. 经典 Runge-Kutta方法K1? f(xn, yn) ?ynhhh2K2? f(xn?, yn?K1) ?yn?yn222hhhh2K3? f(xn?, yn?K2) ?yn?(?yn?yn)22222h2h3?yn?yn?yn)24K4? f(xn? h, yn? hK3)h2h4?(yn?hyn?yn?yn)24223于是,hyn?1? yn?(K1? 2K2? 2K3? K4)6111234? y1?h?(?h) ?(?h) ?(