8.4 线性多步法

线性多步法matlab程序,数学实验“线性多步法（数值积分法,Taylor展开法）”实验报。

数学实验“线性多步法(数值积分法,Taylor展开法)”实验报告(内含matlab程序).西京学院数学软件实验任务书课程名称数学软件实验班级数0901学号0912020107姓名李亚强实验课题线性多步法(数值积分法，Taylor展开法)实验⽬的熟悉线性多步法(数值积分法，Taylor展开法)实验要求运⽤Matlab/C/C++/Java/Maple/Mathematica等其中⼀种语⾔完成实验内容线性多步法(数值积分法，Taylor展开法)成绩教师实验⼆⼗五实验报告实验名称：线性多步法(数值积分法，Taylor展开法)。

实验⽬的：进⼀步熟悉线性多步法(数值积分法，Taylor展开法)。

实验要求：运⽤Matlab/C/C++/Java/Maple/Mathematica等其中⼀种语⾔完成程序设计。

实验原理：数值积分法：常微分⽅程初值问题(1)数值解法中，某⼀步的公式不仅与前⼀步解的值有关，⽽且与前若⼲步解的值有关，利⽤前⾯多步的信息预测下⼀步的值，这就是多步法的基本思想，可以期望获得较⾼的精度。

将(1)中的⽅程在区间上积分，可以得到：⽤等距节点的插值多项式来替代被积函数，再对插值多项式积分，这样就得到⼀系列求积公式。

⽤梯形⽅法计算积分项代⼊(1)中得：设由个数据点 构造插值多项式，这⾥，运⽤插值公式有：得到下列计算公式：(2)其中，由此可得(2)中的系数。

公式(2)是⼀个r+1步的显式公式，称为Adams显式公式。

Taylor展开法：基于数值积分可以构造出⼀系列求解常微分⽅程的计算公式，下⾯介绍基于 Taylor 展开的待定系数法，它可灵活地构造出线性多步法。

对固定的系数，可以选取待定系数使线性多步法的阶尽可能⾼。

还可以根据需要，确定显式还是隐式。

设构造如下具有 p 阶精度的线性多步公式：(4)它们的局部截断误差为:即： (5)对(5)式的右端各项在点处作Taylor展开有:代⼊(5)中得：使的系数为零，得到关于和的线性⽅程组：⽽且得到线性多步法的局部截断误差：实验内容：%数值积分法function [k,X,Y,wucha,P]= Adams4x(funfcn,x0,b,y0,h)x=x0; y=y0;p=128; n=fix((b-x0)/h);if n<5return;endX=zeros(p,1);Y=zeros(p,length(y)); f=zeros(p,1);k=1; X(k)=x; Y(k,:)=y';for k=2:4c1=1/6;c2=2/6;c3=2/6;c4=1/6;a2=1/2; a3=1/2;a4=1;b21=1/2;b31=0;b32=1/2; b41=0;b42=0;b43=1;x1=x+a2*h; x2=x+a3*h;x3=x+a4*h; k1=feval(funfcn,x,y);y1=y+b21*h*k1; x=x+h;k2=feval(funfcn,x1,y1);y2=y+b31*h*k1+b32*h*k2;k3=feval(funfcn,x2,y2);y3=y+b41*h*k1+b42*h*k2+b43*h*k3;k4=feval(funfcn,x3,y3);y=y+h*(c1*k1+c2*k2+c3*k3+c4*k4); X(k)=x; Y(k,:)=y; endX;Y;f(1:4)=feval(funfcn,X(1:4),Y(1:4));for k=4:nf(k)=feval(funfcn,X(k),Y(k));X(k+1)=X(1)+h*k;Y(k+1)=Y(k)+(h/24)*((f(k-3:k))'*[-9 37 -59 55]');f(k+1)= feval(funfcn,X(k+1),Y(k+1));f(k)=f(k+1); k=k+1; endfor k=2:n+1wucha(k)=norm(Y(k)-Y(k-1)); k=k+1; endX=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1, wucha=wucha(1:n,:); P=[n',X,Y,wucha']; - 0 -。

1.1 求下列各数的具有四位有效数字的近似值, 并指出其绝对误差限和相对误差限)1.0ln(,121,1011,1014321====x x x x1.2 下列各数都是对准确值进行四舍五入得到的近似值, 指出它们的绝对误差限、相对误差限和有效数字的位数。

3*5*4*3*2*1100.5,5000,50.31,3015.0,0315.0⨯=====x x x x x1.3 为了使31的近似值的相对误差不超过0.1%, 问应取几位有效数字?1.4 怎样计算下列各题才能使得结果比较精确?(1) x x sin )sin(-+ε,其中ε充分小 (2) ⎰++121N Nx dx,其中N 是充分大的正数(3)xxsin cos 1-,其中x 充分小(4) o 1cos 1- (5) 1001.0-e(6) )11010ln(84--1.5 求方程01562=+-x x 的两个根, 使至少具有四位有效数字。

2.1 证明方程043=-+x x 在区间[1,2]内有且仅有一个根。

如果用二分法求它具有五位有效数字的根，试问需对分多少次？(不必求根)2.2 用二分法求方程0134=+-x x 在[0.3, 0.4]内的一个根, 精度要求21021-⨯=ε。

2.3 找出下列方程的有根区间，选择适当的初始点用二分法求方程的根，精度要求210-=ε。

(1) 02=--x x ；(2) 06cos 2=-++-x e x x ； (3) 01tan =--x x ； (4) 0sin 2=--x e x 。

2.4 考虑方程032=-x e x ，将其改写为3xex ±=，取00=x ，用两种迭代公式迭代，分别收敛到1.0和-0.5附近的两个根(取精度要求310-=ε)。

2.5 为求方程0123=--x x 在5.1=x 附近的一个根，建立下列形式的迭代公式：(1) 2121111kk x x xx +=⇒+=+，；(2) 3212311k k x x x x +=⇒+=+，；(3) 111112-=⇒-=+k k x x x x ，。

a
k 0
r
k
yn k h
k 1

r
k f nk
其中, f nk f ( xnk , ynk )(k 1,0,1,..., r ).
当 10 时，为隐式公式；1=0 则为显式公式。
2019/1/30 2
构造线性多步法的主要方法：数值积分法和泰 勒展开法。 基于数值积分的构造法 将 y f ( x, y ) 在 [ xn r , xn 1 ] 上积分，得到
2019/1/30
5
若积分
xn1

xn1
xnr
y' ( x)dx 用节点
xn , xn1,, xnq
作为积分点，则有

xnr
y' ( x)dx h[a0 y' ( xn ) a1 y' ( xn1) aq y' ( xnq )] hRn1
局部截断误差
q
y' ( xn ) f ( xn , y( xn ))
y( xn1) y( xnr )

xn1
xnr
f ( x, y( x))dx
f ( x, y( x))dx ， 只要近似地算出右边的积分 I r xnr 则可通过 yn 1 yn r I r 近似y(xn+1) 。而选用不同近 似式 Ir，可得到不同的计算公式。
2019/1/30 3

xn1
3 9 y( xn1 ) y( xn2 ) hy( xn1 ) hy( xn1 ) 4 4
2019/1/30 4
(n 0,1,2,...)
(2) 对右端定积分采用梯形 公式, 略去误差项有 h yn 1 yn （f ( xn , yn ) f ( xn 1 , yn 1 )) 2 这就是熟知的梯形公式 。 有望提高精度。 (n 0,1,2,...)

方程组
(6.1)

uukk

1
ak 1uk 1 ak 1uk


a1u1 a1u2

a0u0 a0u1
n0
n1


2 方程的解法
首先，令E表示位移算子，即Eun un1, E iun uni ,i 0,1,2, ,若
L Ek ak 1E k1 a1E a0I
un=c12n+c23n,
由u0=2,u1=5,得
c1 + c2 =2 2c1 + 3c2 =5
2=un=c12n+c23n ( n=0) 5=un=c12n+c23n (n=1)
解得c1=1,c2=1,所以差分方程的通解为un=2n +3n,n=0,1,…。
6.2 线性多步法的方法稳定性
设有线性多步法 引进多项式
它的特征多项式是 p() 2 4 5 ，其根为 1 1,2 5 ，不满足 根条件，因此，该方法虽是三阶方法，但不具有收敛性，不能使用. §5的其它方法都满足根条件。
说明：满足根条件的线性多步法，当 h 0 时，数值解收敛到 准确解。不满足根条件的线性多步法，舍入误差会随计算步数增加 而呈指数增长，截断误差也会随步数增加呈指数形式影响解。把 满足根条件的线性多步法，称为稳定的方法。
替代复根a+bi及a-bi，
因此差分方程(6.1)的通解中含有项
c1,c（2为3）任若意常r~ 是数特，n(c征1c多osa项n2式bp2c(，2rs)in的nma重)rct根g ，ab 。则Lr n (r r~)m q(r) rn ,
其中q(r)是某个次数为k-m的多项式。于是对任意的 0 l m,

Created with SmartPrinter trail version
2. l 步隐式Adams（亚当姆斯）方法(内插法)
公式: 用过点 ( xk i , fk i ),i = 1,0,L,l 1 的插值多项式
~p(x) = l 1 fk
i= 1
l1 x x i
fk 5
fk
fk 2) 1 + fk
2)
yk+1 = yk + 7h20(251fk+1 + 646fk 264fk 1 + 106fk 2 19fk 3)
表 8.3
显式Adams方法
步数(l ) 方法阶 ( p)
公式
1
1
2
2
3
3
4
4
yk+1 = yk + hfk
yk+1
=
yk
+
h(3 2
fk
fk 1)
= xk+1 [ f (x, y(x))
i=0
p( x)]dx
xk
t=x
xk
= xk+1 y(l+1)( (x))(x
xk
l!
xk )(x
xk 1)L(x
xk l+1)dx
h
=
ቤተ መጻሕፍቲ ባይዱ
y(l +1)(
( x )) {
1
t(t + 1)(t + 2)L(t + l
1)dt }hl +1
l!
0
Created with SmartPrinter trail version
xx xxkk i

(
xk
1
,
y(xk1)) 4 f (xk ,
y(xk ))
f ( xk1,
y( xk 1 ))]
1 90
h5
y(5)
(
k
)
用 yk1, yk1 近似 y( xk1), y( xk1) ，得Simpson方法
yk 1

yk 1

h( 3
fk1

4
fk

fk 1 )
局部截断误差为
注（2）预估－校正－改进方法 同预估－校正－改进Adams方法一样，利用两个同阶的显式方 法和隐式方法，即四阶Milne方法（显式）和四阶Hamming方法 （隐式）可以合成以下预估－校正－改进方法.
则(5.12)确定的线性多步法至少为2l-1阶的.隐式类似. 但是线性 多步法与单步法有很大的不同，并不是所有的线性多步法都可以
应用.这是因为由(5.16)得到的方程组的解的情况不定.另外，常 用四步四阶Adams显式公式、三步四阶Adams隐式公式、 Milne(米 尔尼)公式、四阶Hamming（哈明）公式也可由Taylor展开得到.
解得 0 4,1 5, 0 4, 1 2,
且
C4

1 (1 4!
1

4
1)

1 (1 24

5

4
2)

1 6

0
因此得阶数最高的二步显式线性多步法为
yk1 4 yk 5 yk1 h(4 fk 2 fk 1 )
(5.17)
其主局部截断误差为
y( xk1) yk1
l 1
l 1
y( xk1 ) i y( xki ) h i f ( xki , y( xki ))

第七章 常微分方程数值解法常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。

另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。

因此，有必要探讨常微分方程初值问题的数值解法。

本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。

第一节 欧拉法求解常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy（1）的数值解，就是寻求准确解)(x y 在一系列离散节点 <<<<<n x x x x 210 上的近似值 ,,,,,210n y y y y{}n y 称为问题的数值解，数值解所满足的离散方程统称为差分格式，1--=i i ix x h 称为步长，实用中常取定步长。

显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。

定理 设函数()y x f ,在区域+∞≤≤-∞≤≤y b x a D ,:上连续，且在区域D 内满足李普希兹(Lipschitz)条件，即存在正数L ，使得对于R 内任意两点()1,y x 与()2,y x ,恒有()()2121,,y y L y x f y x f -≤-则初值问题(1)的解()x y 存在并且唯一。

一、欧拉法(欧拉折线法)若将函数)xy 在点nx处的导数()n x y '用两点式代替， 即()hx y x y x y n n n )()(1-≈'+，再用n y 近似地代替()n x y ,则初值问题(1)变为⎩⎨⎧==++=+ ,2,1,0),()(001n x y y y x hf y y n n n n(2)(2)式就是著名的欧拉(Euler)公式。

以上方法称 为欧 拉法或欧拉折线法。

河北联合大学第9章常微分方程初值问题数值解法§9.3 龙格-库塔法 §9.4 线性多步法1. 什么是p 阶Runge-Kutta 法？写出经典的四阶Runge-Kutta 法。

答：考虑增加计算),(y x f 的次数，就有可能构造出精度较高的近似公式。

构造如下形式的计算公式：),,3,2( ),( ),(11111p i K b h y h a x f K y x f K K c h y y i j j ij n i n i n n p i i i n n （1）其中i ij i c b a ,,都是参数，确定它们的原则是使近似公式在),(n n y x 处的Taylor 展开式与)(x y 在n x 处的Taylor 展开式的前面的项尽可能多地重合。

将式（1）称为p 阶显式Runge-Kutta 方法，简称R-K 方法。

取4 p ，通过更复杂的计算，式（1）成为),( )2,2( )2,2( ),( )22(6342312143211hk y h x f K K h y h x f K K h y h x f K y x f K K K K K h y y n n n n n n n n n n 上式称为经典形式的四阶R-K 公式，是常用的四阶R-K 公式。

局部截断误差为5()O h 。

2. 如何导出线性多步法的公式？它与单步法有何区别？答：（1）利用Taylor 展开构造线性多步方法线性多步法(Linear-Multistep Method)的一般形式为ri i n i r i i n i n f h y y 101 （2）其中i i ,均为常数，),(k k k y x f f ),,,1(r n n n k 。

若01 ，公式（2）为显式法；01 ，公式为隐式法。

利用Taylor 展开导出线性多步公式的基本方法是：将线性多步公式（2）在n x 处进行Taylor 展开，然后与1() n y x 在n x 处的Taylor 展开式相比较，要求它们前面的项重合，由此确定参数i i ,。

第四节 解线性方程组的 超松弛迭代法
1
第八章 解线性方程组的迭代法
SOR迭代法是Gauss—Seidel 迭代法的一种修正， 可由下述思想得到.
设已知x(k)及已计算x(k+1)的分量xj(k+1) (j=1,2,,i-1).
(1) 首先用Gauss—Seidel 迭代法定义辅助量 x~i(k1),

1
1

1
4 1 1
1 4 1
1 1 4





x2 x3 x4






1 1 1

.
Hale Waihona Puke 它的精确解为x*=(-1, -1, -1, -1 )T.
解 取初始向量x(0)=0，迭代公式为

x1( k x2( k x3( k
x(11) (0.99999646, 1.00000310, 0.99999953, 0.99999912)T ,
(11) 0.46 105. 2
满足误差

(k) 105.
2
迭代次数k
对取其它值，迭代次数如表. 1.0
22
1.1
17
从此例看到，松弛因子选择
1.2
aii
.
j1
ji
i 1,2,, n 也可写作：


x(0) x(k1)
i
( x1(0) ,, xn(0)
xi(k) xi ,
)T
,

i 1
n
xi
(bi

aij
x
(k j

考虑引例1 考虑引例1，对于微分方程 dy = 2x dx 函数y=x2+C是其通解，y=x2+1是其满足初始 是其通解， 函数 是其通解 是其满足初始 条件x 的特解。 条件 0=1，y(x0)=2的特解。 ， 的特解
3
微分方程的(部分 积 微分方程的 部分)积 部分 分曲线(0≤C≤2)见右 分曲线 见右 图，其中红色曲线 为满足初始条件的 特解曲线。 特解曲线。
例8.2 求微分方程 x2 y″+xy′+(x2 -1/4)y=0 的通解 s=dsolve('x^2*D2y+x*Dy+(x^2-1/4)*y=0','x') s = (C1*cos(x)+C2*sin(x))/x^(1/2) 注意 (1)若不加自变量 ，则将把 作为常数求解。 若不加自变量x，则将把x作为常数求解 作为常数求解。 若不加自变量 (2)所得结果可以用命令 所得结果可以用命令pretty进行简化。 进行简化。 所得结果可以用命令 进行简化 例 pretty(s) C1 cos(x) + C2 sin(x) --------------------1/2 x
应满足条件x=1时y=2，因此可解得 又y(x)应满足条件 应满足条件 时 ，因此可解得C=1 故所求曲线方程为y=x2+1。 故所求曲线方程为 。
列车在平直路线上以20米 秒的速度行驶 秒的速度行驶， 例2 列车在平直路线上以 米/秒的速度行驶， 当制动时，列车获得加速度-0.4米/秒2，问开始 当制动时，列车获得加速度 米秒 制动后多长时间能停住， 制动后多长时间能停住，以及列车在这段时间 行驶了多少路程？ 行驶了多少路程？ 设列车开始制动后t秒内行驶了 解:设列车开始制动后 秒内行驶了 米，则 设列车开始制动后 秒内行驶了s米

数值分析
Hale Waihona Puke 8.48.4.1单步法的收敛性与稳定性
收敛性与相容性
数值解法的基本思想是通过某种离散化手段将微分方
程转化为差分方程，如单步法(2.11)，即
y n 1 y n h ( x n , y n , h ),
（4.1）
它在 xn 处的解为 yn ，而初值问题（1.1）,（1.2）在 xn 处的精确解为 y ( x n ) ， (真值减去近似值!)
一般的线性多步法公式可表示为
yn k

i 0
k 1
i
yni h i f ni ,
i 0
k
（5.1）
其中 y n i 为 y ( x n i ) 的近似， f n i f ( x n i , y ni ), x n i x 0 ih,
i , i 为常数. 若 0及 0不全为零，则称为线性 k 步法.
h (55 f n 3 59 f n 2 37 f n 1 9 f n ), 24
p 求值E： f np f ( x , y 4 n4 n 4 ),
校正C： y n 4
h y n 3 ( 9 f np 4 19 f n 3 5 f n 2 f n 1 ), 24
E ( h ) 1 h , 其绝对稳定域由
1 h 1
(参阅教材415页!)
绝对稳定区间为 2 h 0 .
郑勋烨
11
数值分析
用二阶R-K方法解模型方程 y y 可得到
y n 1 ( h ) 2 [1 h ] yn , 2
故
(参阅教材415页!)
相容定理 p 阶方法(4.1)与初值问题(1.1)，(1.2)相容

第8章常微分方程的数值解法8.4单步法的收敛性与稳定性8.4.1相容性与收敛性上面所介绍的方法都是用离散化的方法，将微分方程初值问题化为差分方程初值问题求解的.这些转化是否合理？即当h →∞时，差分方程是否能无限逼近微分方程，差分方程的解n y 是否能无限逼近微分方程初值问题的准确解()n y x ，这就是相容性与收敛性问题.用单步法(8.3.14)求解初值问题(8.1.1)，即用差分方程初值问题100(,,)()n n n n y y h x y h y x y ϕ+=+⎧⎨=⎩(8.4.1)的解作为问题(8.1.1)的近似解，如果近似是合理的，则应有()()(,(),)0 (0)y x h y x x y x h h hϕ+--→→(8.4.2)其中()y x 为问题(8.1.1)的精确解.因为0()()lim ()(,)h y x h y x y x f x y h→+-'==故由(8.4.2)得lim (,,)(,)h x y h f x y ϕ→=如果增量函数(,(),)x y x h ϕ关于h 连续，则有(,,0)(,)x y f x y ϕ=(8.4.3)定义8.3如果单步法的增量函数(,,)x y h ϕ满足条件(8.4.3)，则称单步法(8.3.14)与初值问题(8.1.1)相容.通常称(8.4.3)为单步法的相容条件.满足相容条件(8.4.3)是可以用单步法求解初值问题(8.1.1)的必要条件.容易验证欧拉法和改进欧拉法均满足相容性条件.一般地，如果单步法有p 阶精度（1p ≥），则其局部截断误差为[]1()()(,(),)()p y x h y x h x y x h O h ϕ++-+=上式两端同除以h ，得()()(,,)()p y x h y x x y h O h hϕ+--=令0h →，如果(,(),)x y x h ϕ连续，则有()(,,0)0y x x y ϕ'-=所以1p ≥的单步法均与问题(8.1.1)相容.由此即得各阶龙格－库塔法与初值问题(8.1.1)相容.定义8.4一种数值方法称为是收敛的，如果对于任意初值0y 及任意固定的(,]x a b ∈，都有lim () ()n h y y x x a nh →==+其中()y x 为初值问题(8.1.1)的精确解.如果我们取消局部化假定，使用某单步法公式，从0x 出发，一步一步地推算到1n x +处的近似值1n y +.若不计各步的舍入误差，而每一步都有局部截断误差，这些局部截断误差的积累就是整体截断误差.定义8.5称111()n n n e y x y +++=-为某数值方法的整体截断误差.其中()y x 为初值问题(8.1.1)的精确解，1n y +为不计舍入误差时用某数值方法从0x 开始，逐步得到的在1n x +处的近似值（不考虑舍入误差的情况下，局部截断误差的积累）.定理8.1设单步法(8.3.14)具有p 阶精度，其增量函数(,,)x y h ϕ关于y 满足利普希茨条件，问题(8.1.1)的初值是精确的，即00()y x y =，则单步法的整体截断误差为111()()p n n n e y x y O h +++=-=证明由已知，(,,)x y h ϕ关于y 满足利普希茨条件，故存在0L >，使得对任意的12,y y 及[,]x a b ∈，00h h <≤，都有1212(,,)(,,)x y h x y h L y y ϕϕ-≤-记1()(,(),)n n n n y y x h x y x h ϕ+=+，因为单步法具有p 阶精度，故存在0M >，使得1111()p n n n R y x y Mh ++++=-≤从而有111111111()()()(,(),)(,,)()(,(),)(,,)n n n n n n n p n n n n n n p n n n n n n e y x y y x y y y Mh y x h x y x h y h x y h Mh y x y h x y x h x y h ϕϕϕϕ+++++++++=-≤-+-≤++--≤+-+-1(1)p nMh hL e +≤++反复递推得11111101110(1)(1)1(1)(1)(1)(1)1(1)p p n n n p n n p n e Mh hL Mh hL e hL hL Mh hL e hL Mh hL e hL+++-+++++⎡⎤≤++++⎣⎦⎡⎤≤+++++++⎣⎦+-≤++因为00()y x y =，即00e =，又(1)n h b a +≤-，于是ln(1)1()(1)(1)b a b a hL n L b a h h hL hL e e --++-+≤+=≤所以()11()p L b a p n M e h e O h L -+⎡⎤≤-=⎣⎦推论设单步法具有p (1p ≥)阶精度，增量函数(,,)x y h ϕ在区域G ：, , 0a x b y h h ≤≤-∞<<+∞≤≤上连续，且关于y 满足利普希茨条件，则单步法是收敛的.当(,)f x y 在区域:,D a x b y ≤≤-∞<<+∞上连续，且关于y 满足利普希茨条件时，改进欧拉法，各阶龙格－库塔法的增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，因而它们都是收敛的.关于单步法收敛的一般结果是：定理8.2设增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，则单步法收敛的充分必要条件是相容性条件(8.4.3).8.4.2稳定性稳定性与收敛性是两个不同的概念，收敛性是在假定每一步计算都准确的前提下，讨论当步长0h →时，方法的整体截断误差是否趋于零的问题.而稳定性则是讨论舍入误差的积累能否对计算结果有严重影响的问题.定义8.6若一种数值方法在节点值n y 上有一个大小为δ的扰动，于以后各节点()m y m n >上产生的偏差均不超过δ，则称该方法是稳定的.我们以欧拉法为例进行讨论.假设由于舍入误差，实际得到的不是n y 而是n n n y y δ=+，其中n δ是误差.由此再计算一步，得到1(,)n n n n y y hf x y +=+把它与不考虑舍入误差的欧拉公式相减，并记111n n n y y δ+++=-，就有[]1(,)(,)1(,)n n n n n n y n nh f x y f x y hf x δδηδ+⎡⎤=+-=+⎣⎦其中y f f y∂=∂.如果满足条件1(,)1y n hf x η+≤，(8.4.4)则从n y 到1n y +的计算，误差是不增的，可以认为计算是稳定的.如果条件(8.4.4)不满足，则每步误差将增大.当0y f >时，显然条件(8.4.4)不可能满足，我们认为问题本身具有先天的不稳定性.当0y f <时，为了满足稳定性要求(8.4.4)，有时h 要很小.一般的，稳定性与方法有关，也与步长h 的大小有关，当然也与方程中的(,)f x y 有关.为简单起见，通常只考虑数值方法用于求解模型方程的稳定性，模型方程为y y λ'=(8.4.5)其中λ为复数.一般的方程可以通过局部线性化转化为模型方程，例如在(,)x y 的邻域(,)(,)(,)()(,)()x y y f x y f x y f x y x x f x y y y '==+-+-+略去高阶项，再作变量替换就得到u u λ'=的形式.对于模型方程(8.4.5)，若Re 0λ>，类似以上分析，可以认为方程是不稳定的.所以我们只考虑Re 0λ<的情形，这时不同的数值方法可能是数值稳定的或者是数值不稳定的.当一个单步法用于试验方程y y λ'=，从n y 计算一步得到1()n n y E h y λ+=(8.4.6)其中()E h λ依赖于所选的方法.因为通过点(,)n n x y 试验方程的解曲线（它满足,()n n y y y x y λ'==）为[]exp ()n n y y x x λ=-，而一个p 阶单步法的局部截断误差在()n n y x y =时有1111()()p n n n T y x y O h ++++=-=，所以有1exp()()()p n n y h E h y O h λλ+-=(8.4.7)这样可以看出()E h λ是h e λ的一个近似值.由(8.4.6)可以看到，若n y 计算中有误差ε，则计算1n y +时将产生误差()E h λε，所以有下面定义.定义8.7如果(8.4.6)式中，()1E h λ<，则称单步法(8.3.14)是绝对稳定的.在复平面上复变量h λ满足()1E h λ<的区域，称为方法(8.3.14)的绝对稳定区域，它与实轴的交称为绝对稳定区间.在上述定义中，规定严格不等式成立，是为了和线性多步法的绝对稳定性定义一致.事实上，()1E h λ=时也可以认为误差不增长.(1)欧拉法的稳定性欧拉法用于模型方程(8.4.5)，得1(1)n n y h y λ+=+，所以有()1E h h λλ=+.所以绝对稳定条件是11h λ+<，它的绝对稳定区域是h λ复平面上以(1,0)-为中心的单位圆，见图8.3.而λ为实数时，绝对稳定区间是(2,0)-.Im()h λRe()h λ2-1-O 图8.3欧拉法的绝对稳定区域(2)梯形公式的稳定性对模型方程，梯形公式的具体表达式为11()2n n n n h y y y y λλ++=++，即11212n nh y y h λλ++=-，所以梯形公式的绝对稳定区域为12112h h λλ+<-.化简得Re()0h λ<，因此梯形公式的绝对稳定区域为h λ平面的左半平面，见图8.4.特别地，当λ为负实数时，对任意的0h >，梯形公式都是稳定的.Im()h λRe()h λO 图8.4梯形公式的绝对稳定区域(3)龙格－库塔法的稳定性与前面的讨论相仿，将龙格－库塔法用于模型方程(8.4.5)，可得二、三、四阶龙格－库塔法的绝对稳定区域分别为211()12h h λλ++<23111()()126h h h λλλ+++<2341111()()()12624h h h h λλλλ++++<当λ为实数时，二、三、四阶显式龙格－库塔法的绝对稳定区域分别为20h λ-<<、2.510h λ-<<、 2.780h λ-<<.例8.5设有初值问题21010101(0)0xy y x x y ⎧'=-≤≤⎪+⎨⎪=⎩用四阶经典龙格－库塔公式求解时，从绝对稳定性考虑，对步长h 有何限制？解对于所给的微分方程有2100,(010)1f x x y xλ∂==-<≤≤∂+在区间[0,10]上，有201010max ||max51t x x λ<<==+由于四阶经典龙格－库塔公式的绝对稳定区间为 2.7850h λ-<<，则步长h 应满足00.557h <<.。

同的计算公式。
§4 Multistep Method
亚当姆斯显式公式 /* Adams explicit formulae */
利用k+1 个节点上的被积函数值 fn , fn-1, ... , fn-k 构造 k 阶牛顿 后插多项式 Nk ( xn + t h) , t [0, 1], 有
xn+1 f (x, y(x))dx = xn
… … … … … … …
常用的是 k = 3 的4阶亚当姆斯显式公式
h yn+1 = yn + 24 (55 fn - 59 fn-1 + 37 fn-2 - 9 fn-3 )
§4 Multistep Method
亚当姆斯隐式公式 /* Adams implicit formulae */
利用k+1 个节点上的被积函数值 fn+1 , fn , …, fn-k+1 构造 k 阶
注：上式也可分通项过，数亦值可积得分此导S出im，ps即on将公式y =。f (x, y) 在区间
[xn-3 , xn+1]
再过 fn ,
上fn-积1 , 分fn-，2 做得f到的y插( x值n+1多) =项y(式xn即-3 )可+。xxnn-+31
f
( x,
y( x))dx
,
辛甫生 /* Simpson */ 公式
将通式中的右端各项 yn-1, … , yn-k ; fn+1, fn-1, … , fn-k 分别在 xn 点作泰勒展开，与精确解 y(xn+1) 在 xn 点的泰勒展开作比较。通过令
同类项系数相等，得到足以确定待定系数a0, … , ak ; b-1, b0, … , bk 的等式，则可构造出线性多步法

h yn 1 yn 2 ( K1 K 2 ) K1 f ( x n , y n ) K f ( x h, y hK ) n n 1 2
2013/12/11 18
五、单步法的局部截断误差和阶
设一般的单步法为： 显式公式：yn 1 yn h（xn，yn，h） 隐式公式：yn 1 yn h（xn，yn，yn 1，h） 设 y n 为数值方法的精确值，y（ xn）为微分方程的 精确解。
欧拉（Euler）方法的几何意义
Euler方法有明显的几何意义。 如右图所示，一阶常微分方 程初值问题的解曲线y(x)过 点P0(x0,y0)。从P0出发以 f(x0,y0)为斜率作一直线段， 与x=x1相交于点P1(x1,y1)， 显然有y1=y0+hf(x0,y0)。 同理，再由P1出发以f(x1,y1)为斜率作一直线段，与x=x2 相交于点P2(x2,y2)，显然有y2=y1+hf(x1,y1)。这样一直做 下去，得到一条折线P0P1P2……，作为y(x)的近似曲线。 因此，显示Euler公式又称为Euler折线法。
其中第一式称为预测算式，第二式称为校正算式。
2013/12/11
17
Euler预测校正公式（续） 若将Euler预-校方法中的第一式带入第二式，得
h yn 1 yn [ f ( xn , y n ) f ( xn 1 , yn hf ( xn , yn ))] 2
Remark：这是一种显示的单步方法。有时为了计 算方便，常将上式改写成：
用两点差商公式 代替导数y( xn ) ,再用 y n 表示 y ( xn ) 的近似值,则得到
2013/12/11 8
y ( xn 1 ) y ( xn ) xn 1 xn

i = n, n −1, L, n − k
插值节点的不同取法就导致不同的多步法。
（1）Adams外插法（显式多步法）
取k+1个节点tn-k,…，tn-1,tn及函数值f(tn-i,u(tn-i)) i=k,…,1,0
构造区间[tn,tn+1]上逼近f(t,u(t))的k次Lagrange插值多项式Ln,k(t)
其中
k +1
∑ ( ) u = u + n+1
n h bk +1i f tn−i+1, un−i+1
i=0
=∫ ∏ bk+1i
0 k+1 τ + j dt
−1
j=0 j≠i
j−i
且 t = tn+1 +τ h , τ ∈[−1, 0]。
注： t − tn− j+1 = tn+1 +τ h − tn + ( j −1) h = (τ + j ) h
第1章 常微分方程初值问题数值解法
§2 线性多步法
§2 线性多步法
前节所讨论的方法如Euler方法、改进Euler方法都称为单步法 （单步长法）。 因为它们只利用前一个点的信息来计算下一个点，
即，只用初始点u0计算u1; 一般说来，只用un来计算un+1。
线性单步法一般说来，精度是较低的。 为提高精度，我们考虑
3）内插法是隐式格式（稳定性好），外插法是 显式格式。
2.2 待定系数法（基于Taylor展开式的求解公式） 用数值积分法只能构造一类特殊的多步法,其系数 一般只满足：
ak=1,ak-m=-1 al=0，当l≠k-m, k。
本节我们将基于Taylor展开式来构造出更一般的求 解公式。

讯 ＊ 一ｋ ｈＩ －
展 式得 ｕ＝ ０ｈ（ ，０ ，１ 是 ｕｔ 的近 似值 。Ｎ Ｉｕ又可 算 出 Ｉｕ＋ｆｔ Ｕ）１就 ０ １ （） １ ￣ｌ ｌ
Ｕ ， 此 下 去可 算 出Ｕ 所 有 结 点 上 的值 ， 般 递推 公 式 为ｕ １ ，如 在 一 ＝ 由文 献 ［ ］ ， 解 该 类 方 程分 为单 步法 和多 步 法 。只 要 １知 求 Ｎ ， Ｎ ｕ 即可 算 出ｕ 的算 法 称 为单 步 方 法 ，单 步 法 包 括 Ｊ Ｊ ｔ ｌ Ｎｈ
… ， ｍ＋ｋ－Ｉ才 能 求 Ｕ
出ｕ＋（＞ ） ｋ １ 的算 法 称 为 多 步 方 法 ， 性 多 步 法 包 括 数值 积 分 线 法 和待 定 系 数法 ， 一 般 形 式 为 其 ｕ。 ｈ Ｂ ， 里 ， ： 这
－
（ ０ １Ｋ， ） 常数 ， ｊ ，， ｋ是 ＝ 它需 要 有ｋ 初 值 ｕ， … ，ｋ才 能 得 到 个 。Ｕ ， Ｕ １
０６０ ） ６ ０ ４
用 低 阶微 商 构 造 高 精 度 的 方 法对 实 际计 算 是 很 有 意 义 的 。 ｕ 对 （ １＝ （ ）． ｕ （）ｔ 用ｕ （） 点ｔ １ 处 的 函数 值 及 其 各 ｔ ） ｕ ｔ ＋ ｔｄ利 ｒ ｔ在 ， ｔ
阶 导数 值 构 造 Ｈ ｒｉ 型 插 值 多项 式 ， 代 积 分 。 去 误 差 项 即 ｅｎｔ ｅ 再 舍
【（）ｕ ｕｃ ＝。 ０
㈩
由 于Ｆ 是 已 知 的 ， 是 （ ） ｍ 于 ２ 式可 用 迭 代 法 求 解 。 迭 代 法 都 需 要有 ｕ ， ｎ的 迭 代初 值 。 取 一 个好 的ｕ 即ｕ＋ 选
其 中伪 ｔ 的 已知 函 数 ， 给 定 的初 值 。我 们假 设 函数 和Ｕ ｕ为

计算量大
多一个初值, 可能影响精度
改进欧拉法 /* modified Euler’s method */
§1 Euler’s Method
Step 1: 先用显式欧拉公式作预测，算出 yi+1 = yi + h f ( xi , yi )
Step 2: 再将 yi+1 代入隐式梯形公式的右边作校正，得到
y
K1 = f (xi , yi )
y(xi )
0 1 1
0 2 1
c1 + c2 = 1 o
xi + 1h
xi
xi + h2
x
xi+1
y2 = yi + 2h[c1 f ( xi , yi ) + c2 f ( xi + 1h, y( xi + 1hi )]
y2 = yi + 2h[c1 K1 + c2 f ( xi + 1h, y1 )]
考察改进的欧拉法，可以将其改写为：

yi
+1

K1
= =
1 yi + h[ 2 f (xi , yi
K1 )
+
1 2
K2]

K2 =
f ( xi + h, yi + hK 1 )
斜率 一定取K1 K2 的平均值吗？
步长一定是一个h 吗？
y1 = yi + 1hf ( xi , yi ) = yi + 1hK1
m i=1
yn+1 = yn + h ci f ( xn + i h, y( xn + i h))