数值分析第9章-常微分方程初值问题数值解法.pptVIP

第9章 常微分方程初值问题数值解法 9.2 简单的数值方法与基本概念 9.2.2 梯形方法 9.2.3 单步法的局部截断误差与阶 9.3 龙格—库塔方法 9.3.2 二阶显式R-K方法 9.3.3 三阶与四阶显式R-K方法 　　对欧拉法E(hλ)=1+hλ，其绝对稳定域为|1+hλ|1,绝对稳定区间为-2λ0,在例5中λ=-100,-2-100h0,即0h2/100=0.02为稳定区间，在例4中取h=0.025，故它是不稳定的，当取h=0.005时它是稳定的. 对二阶R-K方法，解模型方程(4.1)可得到 故 绝对稳定域由|E(hλ)|1得到，于是可得绝对稳定区间为-2hλ0，即0h2/λ. 类似可得三阶及四阶R-K方法的E(hλ)分别为 由|1+hλ|1可得到相应的绝对稳定域. 当λ为实数时，则得绝对稳定区间，它们分别为 三阶显式R-K方法： 四阶显式R-K方法： 　　从以上讨论可知显式R-K方法的绝对稳定域均为有限域，都对步长h有限制. 如果h不在所给的绝对稳定区间内，方法就不稳定. 例4 分别取h=0.1及h=0.2，用经典的四阶R-K方法(3.1)计算初值问题 解 本例λ=-20,hλ分别为-2及-4，前者在绝对稳定区间内，后者则不在，用四阶R-K方法计算其误差见下表 0.17×10-4 3125.0 0.15×10-3 625.0 0.14×10-2 125.0 0.12×10-1 25.0 0.93×10-1 4.98 h=0.1 h=0.2 1.0 0.8 0.6 0.4 0.2 xn 　　从以上结果看到，如果步长h不满足绝对稳定条件，误差增长很快. 　　对隐式单步法,可以同样讨论方法的绝对稳定性,例如对后退欧拉法，用他解模型方程可得 故 由|E(hλ)|1可得绝对稳定域为|1-hλ|1，这是以(1,0)为圆心，1为半径的单位圆外部. 故方法的绝对稳定区间为-∞hλ0. 当λ0时，则0h∞，即对任何步长均为稳定的. 对隐式梯形法，它用于解模型方程(4.8)得 故 对Re(λ)0有|E(hλ)|1，故绝对稳定域为μ=hλ的左半平面，绝对稳定区间为-∞hλ0，即0h∞时隐式梯形法均是稳定的. 9.5 线性多步法 在逐步推进的求解过程中，计算yn+1之前事实上已经求出了一系列的近似值y0,y1,?,yn，如果充分利用前面多步的信息来预测yn+1，则可以期望会获得较高的精度. 这就是构造所得线性多步法的基本思想. 构造多步法的主要途径基于数值积分方法和基于泰勒展开方法，前者可直接由方程(1.1)两端积分后利用插值求积公式得到. 本节主要介绍基于泰勒展开的构造方法. 9.5.1 线性多步法的一般公式 如果计算yn+k时，除用yn+k-1的值，还要用到yn+i (i=0,1,?,k-2)的值，则称此方法为线性多步法. 一般的线性多步法公式可表示为 其中yn+1为y(xn+1)的近似，fn+i=f(xn+i, yn+i), 这里xn+i=xn+ih，αi, βi为常数， α0及β0不全为零，则称(5.1)为线性k步法，计算时需先给出前面k个近似值y0,y1,?,yk-1，再由(5.1)逐次求出yk,yk+1,?. 如果βk=0，则(5.1)称为显式k步法，这时yn+k可直接由(5.1)算出；如果βk≠0, 则(5.1)称为隐式k步法,求解时与梯形法(2.7)相同, 要用迭代法方可算出yn+k. (5.1)中系数αi及βi可根据方法的局部截断误差及阶确定，其定义为 定义7 设y(x)是初值问题(1.1), (1.2)的准确解，线性多步法(5.1)在xn+k上局部截断误差为 若Tn+k=O(hp+1)，则称方法(5.1)是p阶的，p≥1则称方法(5.1)与方程(1.1)是相容的. 当r=1, ?(xn, yn, h)=f(xn, yn)时，就是欧拉法，此时方法的阶为p=1. 当r=2时，改进欧拉法(3.1),(3.1)是其中一种，下面将证明其阶p=2. 要使公式(3.4),(3.5)具有更高的阶p，就要增加点数r. 下面我们只就r=2推导R-K方法. 并给出 r=3,4 时的常用公式，其推导方法与r=2时类似，只是计算较复杂. 对r=2的R-K方法，由(3.4),(3.5)式可得如下计算公式 这里 c1, c2, λ2, μ21 均为待定常数，我们希望适当选取这些系数，使公式阶数 p 尽量高. 根据局部截断误差定义，推导出(3.6)的局部截断误差为 其中