第7章节常微分方程数值解法.ppt

略去高阶项，再做变换即可得到的形式 u?=λu. 对于m个方程的方程组, 可线性化为y?=Ay, 这里A为m×m雅可比矩阵(?fi/?yj)，若A有m个特征值λ1,λ2,?,λm，其中λi可能是复数，所以，为了使模型方程结果能推广到方程组，方程(4.8)中λ为复数. 为保证微分方程本身的稳定性，还应假定Re(λ)0. 　　下面先研究欧拉方法的稳定性. 模型方程y?=λy的欧拉公式为 设在节点yn上有一扰动值εn，它的传播使节点值yn+1产生大小为的扰动值εn+1，假设用yn*=yn+εn，按欧拉公式得出yn+1*=yn+1+εn+1的计算过程不再有新的误差，则扰动值满足 可见扰动值满足原来的差分方程(4.9). 这样，如果差分方程的解是不增长的，即有 则它就是稳定的. 这一论断对于下面将要研究的其它方法同样适用. 　　显然，为要保证差分方程(4.9)的解是不增长的，只要选取h充分小，使 　　在μ=hλ的复平面上，这是以(-1,0)为圆心，1为半径的单位圆. 称为欧拉法的绝对稳定域，一般情形可由下面定义. 定义6 单步法(4.1)用于解模型方程y?=λy，若得到的解yn+1=E(hλ)yn，满足|E(hλ)|1，则称方法(4.1)是绝对稳定的. 在μ=hλ的平面上, 使|E(hλ)|1的变量围成的区域，称为绝对稳定区域，它与实轴的交称为绝对稳定区间. 　　对欧拉法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方法计算其误差见下表 xn 0.2 0.4 0.6 0.8 1.0 h=0.1 h=0.2 0.93×10-1 4.98 0.12×10-1 25.0 0.14×10-2 125.0 0.15×10-3 625.0 0.17×10-4 3125.0 　　从以上结果看到，如果步长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∞时隐式梯形法均是稳定的. 7.5 线性多步法 在逐步推进的求解过程中，计算yn+1之前事实上已经求出了一系列的近似值y0,y1,?,yn，如果充分利用前面多步的信息来预测yn+1，则可以期望会获得较高的精度. 这就是构造所得线性多步法的基本思想. 构造多步法的主要途径基于数值积分方法和基于泰勒展开方法，前者可直接由方程(1.1)两端积分后利用插值求积公式得到. 本节主要介绍基于泰勒展开的构造方法. 7.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)相同, 要用迭代