大二下-数值分析

数值分析计算机科学与技术系所EDAcaiy

:

62785564东主搂8区407第九章

常微分方程初值问题数值解法§9.1

引言一、背景常微分方程有着广泛的应用，比如：火箭、的飞行控制；机器人控制；某些化学反应过程的描述与控制。又比如：对一个复杂的集成电路，可根据定律列出一组关于节点电压及支路电流的常微分方程。求解这组常微初值问题即可详细分析电路的瞬态特性。但实际当中，只有极少数常微分方程具有解析解，大多需借助数值方法。常微分方程分两大类：初值问题和边值问题，本章仅 前者。二、常微初值问题在初值问题里，集中 由单个的一阶方程构成的初值问题：(1.1)(1.2)

y

'(x)(

f

(x,

y(x))

f

(x,

y)

y(x

)

y

o

o定义：存在常数L>0，使x

[a,

b],

y1,

y2

R,有：

f

(x,

y1)

f

(x,

y2

)

L

y1

y2

,

称

f

满足对y的Lipschitz（）条件，L称Lipschitz常数。f

满足条件还可用更强的条件保证。如

f

连续有界f对y满

y足Lipshitz条件。这概念在第七章非线性方程求解中提到过。明确一下本章对象：a

x

b,定理1：D

x,

y

y

R连续；（1）

f

(

x,

y)（2）f对y满足Lipschitz条件，则初值问题：

y

f

(x,

y)

o

0

o

o

y(

x

)

y

,

x

[a,

b]，y

R存在唯一解。定理2：设f在区域D（如定理1所定义）上连续，且关于y满足条件，设初值问题：y’(x)=f(x,y)，y(x0)=s的解为y(x,

s)，则：

f

(x,

s1)

f

(x,

s2

)

eL

x

x0

s1

s2该定理表明解对初值的敏感性，当f

的L比较小时，不敏感；L大时为

问题。…

…在带形区域D内存在唯一解，或称积分曲线y(x)，满足微分方程及初值条件。三、数值解法基本思想数值方法的思想是将问题的求解要求降低为只求一系列离散点的函数值。具体地说，是将积分区间离散化，并对离散点函数值的近似值：x0

x1

xn

xn1

yn1y0

y1yn的差分方程。所谓差分方程构成形如：yn1

G(yn1,ynynk)是“若干邻近节点函值应满足的关系（方程）”。或者说是

“节点函数值的某种递推关系”。比如从y(x0

)

y0

及其它可,

yk提供的初值

y1,

y2

,

出发，经差分方程步进地求出

yn1

。yn

yn

1

G

(

yn

1

,

ynk

)概括地说是“将连续问题离散化”，差分方程可分类如下：G

0k

0单步法

不含yn

1

显式多步法yn

1

构造差分方程的基本方法有三种：1、基于数值微分的方法（9.2.1）2、基于数值积分的方法（9.2.2，9.3，9.5.2）3、基于Taylor展开的方法

（9.3，9.5）重要介绍方法2、3的全部差分方程，且着重于方法3。因为原则上，方法3可以导出方法1与方法2的全差方，但方法1、2

有直观性，将比较全面地加以介绍。第九章课后练习(20)作业：九3、4、6、7预习：9.4、9.5§9.2简单的数值方法本节介绍四种简单的差分方程构造方法，事实上，它们可以由以上三种构造方法的任一种推导得到。但为了大家对差分方程的构造方法有较全面、完整的了解，第1小节的向前Euler方法与后退Euler方法重点从基于数值微分的方法去推导。9.2.1

（向前）Euler与后退（向后）Euler法一、基于数值微分推导向前Euler公式基于数值微分构造差分方程，有两个重要事实要经常用到：（1）积分曲线上任—点(x,y(x))，满足微分方程，即：y(x)

f

(x,

y(x)),x

[a,

b]（2）节点导数节点数值导数下面推导Euler向前公式：对初值问题，起始点x0的函数值y(x0)已知，那么对一个选定的h：x1

x0

h步长 ，如何构造联系x0与x1两个节点函数值的关系式，即关于节点函数值的差分方程呢？积分曲线在x0满足的微分方程，即以x0代入微方程：y(xo

)

f

(xo

,

y(xo

))左端是x0的节点导数。可选用最简单的数值微分公式：P131，第5行的两点公式近似它：0

10(x

)

1

[

x

x

]

h

h

2只要将其中的符号作相应修改：

f

y

即可。0

10

0

0h

2y(x

)

1

[

y(x

)

y(x

)]

h

y

()

f

(x

,

y(x

))即：1

0

0

0h

21

[

y(x

)

y(x

)]

h

y

(

)

f

(x

,

y(x

))h

y()若设：y(x)C[a,b]（连续）且积分步长h足够小，可将余2略去，则有：y(x1

)

y(x0

)

hf

(x0

,y(x0

))这样，即建立起y(x0)与y(x1

)的递推关系式，由已知初始值y(x0)即可经由差分方程导出

y(x1

),

y(x2

)

，可得一般递推公式：y(xn1

)

y(xn

)

hf

(xn

,

y(xn

))不过仔细检查上面的推导过程会发现，在符号使用不上够严格，这是因为前面已约定：y(x)——积分曲线，是方程的准确解，那么，y(xn

)——

y(x)

在xn

的准确解。但是，上述递推关系式是积分步长足够小，以致数值导数误差可以忽略的条件下得到的。因此，即使

y(x0

)

完全准确，而经过递推公式计算的号使用上的

，y(x1

)

一般就不再准确了。为了消除这种符引入以下记号加以区别：yn——差分方程在xn

的准确值。上述递推 系应 写为：yn1

yn

hf

(xn

,

yn

)——向前Euler公式这是形如yn1

G(yn

)的差分方程。由于yn+1仅依赖于前面一个函数值yn，故称单步法。而且，一旦有了yn，yn+1就可算出，又称显式公式，因此，向前Euler公式为：显式、单步的积分公式。小结：基于数值微分构造求积公式可概括为：“将微分方程离散化，以数值导数近似节点导数，从而构造差分方程。”二、基于数值微分推导向后Euler公式同样是基于数值微分方法推导，与向前Euler构造相似，y(x0)已知，对选定的离散化步长h，有x1

x0

h

，要建立联系x1，x0函数值关系的差分方程。积分方程在x1满足的微分方程：y(x1

)

f

(x1

,

y(x1

))左端的一阶节点导数，可选用P131第6行的两点公式：1

1oh

2f

(x

)

1

[

f

(x

)

f

(x

)]

h

f

()类似地，将符号做相应改变，即

f

y

，有：1

1o

1h

2y(x

)

1

[

y(x

)

y(x

)]

h

y

(

)

f

(x,

y(x

))略去数值导数的误差，又用y1~

y(x1

)，有：y1

yo

hf

(x1

,

y1)得到的非线性代数方程。求出y1后，又可求出y2,y3,……一般地有：yn1

yn

hf

(xn1,

yn1

)——向后Euler公式向后Eluer公式为：隐式、单步法。小结：向前与向后Euler的推导区别仅在于：y

'(xn1

)

f

(xn1

,

y(xn1

))y

'(xn)

f

(xn,

y(xn))——向后Euler——向前Euler即向后Euler是 微分方程在离散点：xn+1的关系。而向前Euler是微分方程在离散点：xn

的关系。三、向后Euler公式的计算向后Euler公式是一个隐式公式，就是说一般从yn计算yn+1时不能像显式公式那样直接计算得到，而需求解一个关于yn+1的非线性代数方程。第七章已经初步介绍求解单变量非线性代数方程的方法，着重在于通过不动点迭代求解。看到，向后Euler公式本身恰好是以某个函数的不动点形式表现的，因为：yn1

yn

hf

(xn1,

yn1

)

(

yn1

)即yn1

(yn1

)的不动点表达形式，并可立即得到相应的不动点迭代格式：(0)n

n,

y(k

)

)(

k

1)n1n

n1

n1

y

y

hf

(xy

y

hf

(x

,

y

)

n1nk

0,1,

(x)

C[a,

b]定理1，

设 满足以下两个条件：复习第七章定理有1、对任意

x

[a,

b]

有

a

(x)

b2、存在正常数L<1，使对任意x,y

[a,b]，(x)

(

y)

L

(x

y)*则

(x)在[a,

b]上存在唯一不动点x*

。定理2，设(x)

C[a,b]满足定理1两个条件，则对任意初值

xo

[a,

b]

，由

xk

1

(xk

)

得到的迭代序列

xk

收敛到

(x)

的Lk(2.5)

。不动点x*，并有误差估计：xk

x

1

L

x1

x0,

y

(

k

)

)y

(

k

1)n

1

yn

hf

(

xn

1设f(x,y)对y满足Lipschitz条件：f(x,y)

f

(x,y)

L

y

y

yn1存在唯一（y(x)存在唯一）利用第七章定理2：（1）

y

R,（2）

Lh

1

y

时，…n

1yn

1

yn

hf

(

xn

1

,

yn

1

)y(

k

1)n1n1n1

n1

y

h f

(x,

y(k

)

)

f

(x

,

y

)n1

n1

hL y(

k

)

yn1n1只要h足够小，使hL<1成立，即可保证迭代法收敛。9.2.2

梯形方法一、基于数值积分的推导：积分：对微分方程y(x)

f

(x,

y(x))，在区间[xn

,xn1

]xxxnxnf

(t,

y(t))dtn1n1y(x)dx

n1ny(xx)

y(

)

xn1f

(t,

y(t))dtxn这里要把积分作近似处理：

ft可利用递形公式计算右端积分即有：2n1

nn

n

n1n1y(x

)

y(x

)

h

[

f

(x

,

y(x

))

f

(x

,

y(x

))]考虑到数值积分不可避免地存在误差，可转化为差分方程形式：2n1

n

n

nn1

,

yn1

)]y

y

h

[

f

(x

,

y

)

f

(x——梯形公式梯形公式为：隐式、单步法。注意：积分区间的取法以及积分 的数值积分处理方法不同，可得到更一般的基于数值积分方法推导的差分格式。比如：仍取积分区间

[xn,

xn1]

，但积分

以左端点的左矩形公式近似，则有：向前Eulery(xn1

)

y(xn

)

hf

(xn

,

y(xn

))

yn1

yn

hf

(xn

,

yn

)f同理，可推出向后Euler公式。yn1

yn

hf

(xn1,

yn1

)t二、梯形公式的计算P283，式（2.8）已构造了一个不动点迭代格式，也做了收敛性分析。n1n

hf

(xn

,

yn

)

y（0）

y(2.8)n

nn1n1

ny

y

h

f

(x

,

y

)

f

(x2

n1

（k+1）,

y（k）)(k

0,1,

2,

)

9.2.3

改进的Euler公式

y梯形公式：n1n

n

n

n1

n1

y

h

[

f(x

, )

(x

,

)]2是一个隐式公式，前面已经，给出初值y(0)后不动点迭代n12nn

n1n1n

n1y(k

1)

y

h

[

f

(x

,

y

)

f

(x

,

y(

k

)

)]给出序列

(

k

)yn1

只要h足够小，有：y(

k

)n1n1k

yyn1

是梯形差分公式精确解。但上述迭代过程一般是个计算量很大的过程，而且从另一角度看，即使准确的

yn1

，也仅是y(xn1

)

的近似值

yn1

y(xn1

)

。因此，耗费大量机时去获得一个不一定很准的函数近似值并不值得，常常采用如下计算过程：n1

y(0)

yn

hf

(xn

,

yn

)2nn

n,y(k

)

)]n1

n1

h

[

f

(x

,

y

)

f

(x(k

0,1,

2,)y(k

1)

y

n1即仅作少量的几步迭代，称为“预报-校正”格式，甚至往往更加简单，只取k=0，即：n1

y(0)nn

n

y

hf

(x

,

y

)2(0)n1

n1h

yn1

yn [

f

(xn

,

yn

)

f

(x

,

y

)]——改进Euler公式（P283

（2.13））也可改写为一个的显式单步公式：2n

n

n

nn

nyn1

yn

h

[

f

(x

,

y

)

f

(x

h

,

y

hf

(x

,

y

))]还可以改写成与R-K公式更相近的表达形式：2

k1

y

n

h

f

(

x

n

,

y

n

)

k

y

n

h

f

(

x

n

1

,

k

1

)y

n

1

1

(

k

k

)2

1

29.2.4

单步法局部截断误差与阶单步法的一般形式：yn1

yn

h(xn

,

yn

,

yn1,

h)(2.11)

h

(xn

,

yn

,

h)显示单步法：yn1

yn与f(x,y)有关，比如：(xn

,

yn

,

h)

f

(xn

,

yn

)向前Euler：梯形公式：2n

n

n1

n

nn1

n1(x

,

y

,

y

,

h)

1

[

f

(x

,

y

)

f

(x

,

y

)]下面引出局部截断误差概念，用于判断一个差分公式的精度，是常微数值解的一

十

基

的 念。Tn1

y(xn1

)

y(xn

)

h(xn

,

y(xn

),

h)定义1，设y(x)是初值问题（1.1），（1.2）的精确解，称为显式单步法：yn1

yn

h(xn

,yn

,h)的局部截断误差。1、向前Euler局部截断误差yn1

yn

hf

(xn

,

yn

)yn1

y(xn

)

hf

(xn

,

y(xn

))向前Euler：若令：则不难看出，二者是有区别的：yn+1表示从y(x0

)出发应用向前Euler差分格式，从yn求出yn+1。一般来：yn

y(xn

)

。而

yn1表示从xn求yn+1

时，使用了积分曲线在xn的准确值y(xn)，并应用向前Euler差分方程得到的xn+1的函数近似值

yn1

。那么，局部截断误差正是

：Tn1

y(xn1

)

yn1

y(xn1

)

y(xn

)

hf

(xn

,

y(xn

))应如何比较右端各量，以便将右端简化呢？基本方法是将各量在xn作Taylor展开。一般都取在xn作Taylor展开，即：y(xn1

)

y(xn

h)

y(x

n)

hy

(x

n)

y

(xn

)

h

22又注意到：Tn1y(xn1

)

yn1

y(xn1

)

y(xn

)

hf

(xn

,

y

xn(xn

, (xn

))

(xn

)y

(xn

)

y(xn

)

hy(xn

)h

2代入上式，有：Tn1

y(xn1

)

yn1

y(xn

)

hy(xn

)

2n

h

2

y(x

)

2

h

2

y(x

)2n——向前Euler局部截断误差首P=1，一阶显示公式定义2，设

y(x)

是（

.

）（.

）的准确解，若存在最大整数P使显式单步法的局部误差满足：Tn1

y(x

h)

y(x)

h(x,

y,

h)

O(hp1

)具有P阶精度。(2.10)称方法yn1

yn

h(xn

,

yn

,

h)n1

(x

p1

p2n

,

y(xn

))h

O(h

)若：T

(xn

,

y(xn

))hp1称：为局部截断误差主首。2、向后Euler局部截断误差向后Euler公式：yn1

yn

hf

(xn1,

yn1

)yn1

y(xn

)

hf

(xn1,

yn1

)记：表示，在xn时，取y(xn)为准确值递推得到的xn+1的差分方程。与推导向前Euler公式局部截断误差相似，应取：数值yn1Tn1

y(xn1

)

yn1

y(xn1

)

y(xn

)

hf

(xn1

,

yn1

)作为向后Euler局部截断误差。一般来说但书上（P285）为了推导局部截断误差主yn1

y(xn1

)比较方便以y(xn1

)得到：代替yn1(xn1

)

(xn

)

h

(xn1, (xn1

))Tn1

从微方程可知：f

(xn1,

y(xn1

))

y(xn1

)

y(xn

h)h

2则有：h

232Tn1

y(xn

)

hy(xn

)

2

y

(xn

)

y(xn

))

O(h

)h[

y(xn

)

hy(xn

)

y

(xn

)]

h

2

y

(x2nP=1，是一阶隐式公式。3、

形公式局部截断误差2,

y

)]n1

nn

nn1

n1y

y

h

[

f

(x

,

y

)

f(x梯形公式：主的方法定义梯形P285仍取便于推导局部截断误差首公式局部截断误差：2,

y(x

))]n1

n1

nn

nn1n1T

y(x

)

y(x

)

h

[

f

(x

,

y(x

))

f

(xh

3y(xn

)首

为：

12P=2

，是二阶隐式公式。§9.3

Runge-Kutta方法本节介绍基于Tayler展开方法推导Runge-Kutta公式为了便于理解，从简单的2阶公式的构造入于，因此，1、2两小节对—

。9.3.2

二阶显式R-K方法一、一般二阶公式约定梯形公式的推导， 是通过与微分方程等价的积分方程：nxf

(t,

y(t))dtxn1y(xn1

)

y(xn

)

还利用积分方程形式去理解如何约定一般给出的。下面的二阶R-K公式。通过等价积分方程构造差分格式，关键是用不同方法对积分做数值积分：ynk1通过两种简单的显式单步公式去理解一般二阶R-K公式。1、向前Euler公式n

nyn1

yn

hf

(x

,

y

)ynk1令：

k1

f(xn

,

yn

)yn1

yn

hk1则：K1是函数f(x,y)在左端点xn的函数近似值。向前Euler可看作用左矩形公式处理积分 得到的最简单的差分公式。2、改进Euler公式2n

nnnn

nn1

ny

y

h

[

f

(x

,

y

)

f

(x

h,

y

hf

(x

,

y

))]若令：k1

f

(xn

,

yn

)

hk1

)k2

f

(xn

h,

ynynk11

1

yn1

yn

h[

k1

k2

]

2

2容易看出，K1仍是f(x,y)在左端点xn的函数近似值，以红点为标识；而K2则可看作f(x,y)在右端点xn+1的函数近似值。因此，改进Euler公式可以看作用梯形积分方法近似处理积分而左、右点函数近似值为K1与K2所给出。将上述数值积分处理积分方程中积分 的思想一般化可在以下两方面加以推广，从而给出一般二阶公式的约定：

h(C1K1

C2

K2

)

yn1

yn1

n

nn2

n

21

1K

2K

f

(x

，y

)

f

(x

h，y

hk

)2n1

n

y

y

h[

1

k

1

k

]12

2k1

f

(xn

,

yn

)12C1

C2

2k

f

(xn

h,

yn

hk1)221

1

h(C1K1

C2

K2

)

yn1

yn

1n

n

2

nK

f

(x

，y

)，K

f

(x

h，y

hk

)2

n

21

1与前面向前Euler与改进Euler公式相比，所作两点扩展表现在：（1）K2的扩展从图中看出：22

h

(xn

,

xn1

],

0

1K2可以取

(xn

,

xn1

]

中任一点：

xn（2）K1与K2的线性组合，（或理解为数值积分公式的扩展，不仅局限于左矩阵或梯形积分模式）。引入线性组合系数C1与C2，在一般二阶R-K公式中含2

,21

,C1,C24个待定常数。二、由Tayler展开定

2

,21

,C1,C2目标是适当选择四个待定参数，使yn+1在xn的Taylor展开，具有与y(xn+1)在xn的直接Taylor展开有尽可能高的相同的阶，也就是使yn+1的阶数尽可能高。为 推导时

写 ，令：yn

y(xn

)，fn

f

(xn

,

y(xn

))

f

(xn

,

yn

)局部截断误差可以写为：nnTn1

(xn1

)

h[C1

(xn

,

n

)

C2(xn

2h,

21h

n

)]f(x,y)作二元函数在局部截断误差表达式中，先看最后一Taylor展开：(xn

2h,

h

n

h

)

'

x

'

xn

h

O

h2n

21

n

n

x

n

n

2

y

21

n三那么，局部截断误差展开式前出现

fn，fx

'

h，fy

'

fnh三

应表达为：C2

f

xn

2h,

yn

21hfnTn1y(xn1

)

yn

h[C1

f

(xn

,

ynh

2h

3Tn1

yn

hyn

2yn

y3!

n

yn

O

h2

C1hfn

C2h fn

fx

'2h

fy

'

21hfny

'n

fn其中：dxy

d

f

(x

,

y(x

))

f'

(x

,

y

)

f

'

(x

,

y

)

y(x

)n

nn

x

n

n

y

n

n

n

f

'

(x

,

y

)

f

'

(x

,

y

)

fx

n

n

y

n

n

n

f

'x

f

'y

fn（P287，（3.8））合并：fn，fx

'h，fy

'fnh

三f

y

'

fn

代入Tn+1表达式，依h

2

h

2Tn1

hfn

2

fx

'

2C1hfnx2

21

y

nh2C2h

n

C22'

h2

C'

212

22

2

x2

21

y

n)

f

'

f

h2

(1

C1

C2

)

fnh

(

C

)

f

'

h

(

1

C

O(h3

)2

2要使Runge-Kutta公式具有2阶精度，应有：1

C1

C2

01

C

0

2（3.9）2

211

C

02三个方程定4个未知数，解不唯一，可以得到一族二阶R-K公式，比如：改进Euler公式是当：2

21

11

22C

C

1—

二

R-K公式，类似可以构造许多其它二阶R-K公式。.

.

显式Runge-Kutta

一回到微分方程的等价积分方程形式：nxf

(t,y(t))dtxn1y(xn1

)

y(xn

)

构造相应差分格式：ri1yn1

yn

h

Ci

Ki

yn

h(C1K1

C

r

Kr)其中：K1

f

xn

,

yni1i

2,

rKi

f

(xn

ih，yn

h

ijk

j

)j

1称为r阶Runge-Kutta法一般形式。从以下的略图理解R-K公式一般形式的构造思想：…通过f(x,y)在r个点上函数近似值做数值积分，而数值积分公式的系数是待定的：C1，……，Cr。括，用数值积分近似积分xnf(t,

y(t))dtxn11

1

2

2

r

rh[C

k

C

K

C

K

]其中，两个因素是待定的：i1j

1yn

h

ijk

jxn

ih，（1）ki

fi即所取作数值积分的函数值为待定。（2）C1，……，Cr，即数值积分公式待定。通过引入大量待定常数，使之达到尽可能高的阶数。类似二阶显式R-K公式的构造与推导，可以构造高阶R-K公式。9.3.3

三阶与四阶R-K公式（了解公式即可）从P288~289看到，在使用三阶、四阶R-K公式时，

一所需计算

f(x,

y)

函数值的次数分别为3次及4次，下面有关于R-K公式阶数与函数值计算的表：R-K阶数≤45计算次数≤4697…….因此，一般取R-K公式≤6阶。第9

章课后练习(21)作业：九8、11、12§

.

单步法的收敛性与稳定性9.4.1

收敛性与相容性一、单步法的收敛性先通过一个简单例子，直观了解一下收敛性的概念。例：用逐次二分积分步长的向前Euler法求解。

y

'

y

0

x

T

y(0)

y0其准确解为：y

y(x)

e

xo相应图形如下：积分区间内任一固定点C：0≤c≤Tc那么，积分曲线在这点的值为：y(c)

y0e下面 使用逐次二分积分步长的向前Euler公式积分上述初值问题时，当积分步长越来越小时，c点处函数近似值的变化情况：（1）将[0,c]分作n等分，则积分步长n

nh

c

0

c（2）取—定n值，从x0出发，用向前Euler一步一步积到c，可得C点近似的函数值yn。yn1

yn

hf

(xn

,

yn

)如，n=1n=2n=4y1y2y4……可得y(c)的近似值序列

yn,n

2

，k为整数。kn

on

k

lim

y

y(c)

y

ec从图上看出题收敛。称向前Euler式积分上述初值问强调一点：在收敛性概念里，c是积分区间任一相对固定的：

x0

c

T

，它不随积分步数或积分步长的变化而变化。这点十分重要，下面给出收敛性定义。n

,

h)定义3，若一种数值方法（如单步法

yn1

yn

h(xn

,

）,

当

k

时，yn

y(c)

称该方对于固定的xn

xo

nh

c法是收敛的。这里对书上定义作两点修改和说明，以避

理

上。

nh（1）xn

x0改为xn

xo

nh

c

，强调是固定点c。（2）ynn

y(x

)h0

n改为

n

h0

，会更确切些。y

n

y(c)书上的表述不很严格，h~n是存在关系的，当h0时，应

有

n

。但一方面

n

，另一方面极限结果仍含n，不够严格。下面显式单步法的收敛性。有P阶精度，且增定理3，设显式单步法yn1

yn

hxn

,

yn

,

h量函数，关于y满足Lipshitz条件：(x,

y,

h)

(x,

y,

h)

L

y

y

y(x0

)又设初值y0是准确的，即

yo

，由其整体截断误差为：py(xn

)

yn

O(h

)一下。在证明前，将一些相关概念及定理结论、条件1、局部截断误差与整体截断误差有几个相关符号再明确一下：y(xn1

)——积分曲线

y(x)

在

xn1

准确值。yn1——

yn1

y(xn

)

h(xn

,

y(xn

),

h)，即设：yn

y(xn

)准确，显式单步法积分一步至的准确值。xn1，从x0出发，显式单

h

xn

,

yn

,

h——

yn1

ynyn1步公式逐步积分至xn+1

的准确值。这里注意：yn1

与yn1

的不同含义，因此，——局部截断误差，前进一步误差。y(xn1

)

yn1en1

y(xn1

)

yn1——整体截断误差，各步误差的积累。可以从图上理解这两种截断误差的含义以及它们之间的关系：y=Tn+1y(xn)ynn12、证

思再把定理条件与结论强调一下：条件：（a）

y(xn1

)

yn1

O(hp1

)，P阶精度。，即增量函数关于y（b）(x,

y,

h)

(x,

y,

h)L

y

y满足Lipschitz条件。（c）

yo

y(x0

)

，即y0准确。即整体截断误差不超过h的P阶小y(x

)

y

O(hp

)n1

n1量。结论：概括地说，定理给出局部截断误差的条件以及其它有关条件，而结论是要求整体截断误差应当满足的关系。因此，应当将整体截断误差与局部截断误差联系起来：y(xn1

)

yn1

y(xn1

)

yn1

yn1

yn1

y(xn1

)

yn1

yn1

yn1其中，第一 为局部截断误差，其误差界是清楚的，问题在第二 的估计。3、证明分别对一、二两 估计：ny(xn1

)

yn1

Chp1yn1

yn1(xn

),

h

n

,

h

y(xn

)

h(xn

,

h

xn

,

y(xn

)

yn

h

(xn

,

y(xn

),

h)

(xn

,

yn

,

h)y(xn

)

yn

L

y(xn

)

yn

(1

hL

)

y(xn

)

yn使用前面引入的整体截断误差记号，有n1e

(1

hL

)

en

chp1是整全截断误差的递推公式，反复递推：n

n

1e

1

hL

e

chp1

(1

hL

)[(1

hL

)+

en

2

chp1

]

chp1=（1

hL）2

en

2

chp1[1+(1+hL

)]nnchp

(1

hL

)ne

chp1[1

(1

hL

)++

（

1

+

hLn-1）]o

(1

hL

)

eo [(1

hL

)

1]L请看P292推导：(1

hL

)n

(ehL

)n

eTL其中，T为积分区间总长度。因此，整体截断误差上界为：n0Le

e eTL

ch

p

(eTL1)nnnLe

y(x

)

y

ch

p

(eTL1)

O(h

p

)如何应用定理1

一个显式单步法积分一个具体问题的收敛性，必须对方法的阶、以及增量函数是否对y满足Lipschiz条件加以分析。P292以改进Euler

法为例对其改敛性作了分析，希望大家掌握它，再给一个思考题：思考：

y

'

f

(x,

y)

y(x

)

y

0设f

(x,y)对y满足Lipschiz条件，y=y(xo

)准确,证明：用P289（3.

）公式做数值积分时收敛。二、相容性定义4，若单步法yn1

yn

h(xn

,yn

,h)的增量函数满足：(x,y,0)

f

(x,y)

，则称单步法与初值问题：

y

'(x)

f

(x,y)

相容

y(x

)

y

o

oP291从显式单步法局部截断误差关于积分步长h的Taylor展开引出相容性概念。

也可以从另一个角度理解相容性含义。设显式单步法：

h(xn

,

yn

,

h)yn1

yn与初值问题：

y

'(x)

f

(x,

y)

y(x

)

y

o

o相容

Tn1

O(hp1

)，

p

1若把单步法改写为：

yn1

ynn

n

(x

,

y

,

h)h两边取极限

h

0

有：

y

'(xn

)

(xn

,

yn

,

0)

f

(xn

,

yn

)这表明：h

0

时，显式单步法格式趋于微分方程本身。这是相容性的直观意义。（1.1）和（1.2）相容的充分必要条件是P1。定理4，p

阶方法

yn1

yn

h

xn

,

yn

,

h

（4.1）与初值问题9.4.2

绝对稳定性与绝对稳定域稳定性是差分方程求解常微分方程另一个重要理论问题，它研究差分方程对误差的性质。它的含义是在积分步长相对固定的条件下，若在某一步引入误差，在继续递推过程中，这个引入的误差会如何？它的性质与什么因素有关？如何防止和控制误差的恶性增长？这些都是稳定性所关心的问题。绝对稳定性是稳定性概念中比较常用、比较基本和经典的。一、例

及通过一个例子，了解稳定性概念和它的基本含义。例4

初值问题：y

100

yy(0)

0为了使问题更直接、更简单，这里把例中初值1改为

o

，即初值本身就表现为引入的某种误差。 分别后Euler

求解这初值问题时表现的不同现象。向前与向积分曲线，即准确解是：y(x)

oe

100

x（1）向前Euleryn1

yn

h(100

yn

)

1.5yn因此，当

yo

o

，h=0.025时y1

1.5oy2

2.25o己放大到P294

表9-4列出向前Euler积分至0.1时，初值的

05.0625

，即放大了5倍多。0h

0.0254-4（2）向后Euleryn1

yn

h(100

yn1

)

yn

2.5

yn113.5nyn1y

看P394：表9-4第3例数据，给出了向后Euler以0.025步长积分同一问题的数值结果，在图上可以看出，初值时引入的误差

0

，在随后各步递推计算过程中，逐渐而这于0，说明初始误差受到控制。这个例子表现了向前Euler

与向后Euler不同的误差

性质，下面给出定义。法是稳定的。用符号表示为定义5，

若一种数值方法在节点值yn上有大小为

的扰动，于以后各节点值ym

(m

n)

上产生的偏差均不超过

,则称该方m

n（m

n）有很大，通常要简化对象，稳定性问题的一般性并称之为模型问题。二、模型方程y

y,

C,称

Re()

0

为模型方程。在模型方程中，假设

Re()

0

是为了使模型问题y

yy(xo

)

yo稳定准 为：

(x)

(

x

)0e

x0模型问题的初值受到扰动时，对模型问题的解产生影响：如果初值扰动产生的影响能逐渐是稳定的，反之是不稳定的。下面说明Re

()

0

保证模型问题稳定。设，有一初值受到扰动的模型问题：y

yy(xo

)

yo

o是与前面给出的模型问题完全相同的系统，区别仅在于后，其准确解为：者在初值带有误差0y

(x)

(

yo

(

x

x

)

o

)e

o那么：e

(

x

xo

)(x)

(x)

oeRe(

)(

x

xo

)

oy(x)

y(x)

当：Re()

oRe()

o稳定不稳定y(x)

y(x)

三、显式单步法公式稳定性1

、向前Euler公式模型方程的向前Euler公式为：yn1

yn

h

yn

(1

h)

yn记第n步，带有误差的解为：yn

*

ynnnnn1递推一步有：y

(x

h)

y

(1

h)(

y

)n

yn1

(1

h)n

(1

h)

n即递推一步的误差为：n1n1n1

y

yn1

1

h

n即：这是向前Euller

公式对误差或扰动的递推关系式，而且可以看出，误差 的关系式与原来的差分格式完全相同。若要求误差不增长，即1

h

1

n1

yn1

n

yn称

1

h

1

为向前Euler的绝对稳定域。=h平面定义6，单步法

yn1

yn

h(xn

,

yn

,

h)

用于的模型方程y

y

(Re()

o)

1

，则称单步法是绝对稳定称为绝对稳定域，它与实轴的满足

E

hE(h)

1若得到的解yn1

E(h)yn的。在h

平面上，区域交称绝对稳定区间。所以，对向Euler公式，绝对稳定区间是：2

h

Re()

0

0

h

2

o

R回到例4

100,

2

100h

0应有

0

h

2

0.02100是其绝对稳定区间。因此，当取0.025时是不稳定的。2、二阶R-K公式（

.

）关键是推导模型方程下的二阶R-K公式：K1

yn

h2

h

y

)

y2K2

(

yn2nnyn（h）22

hk

h2nnn那么：n1y

yn)

ynh

2

(1

h

2h

2

E(h)

1

h

2因此，绝对稳定区域为：2h

21

h

1三、四阶R-K方法自己推导一下相应的表达式四、隐式单步法稳定性1．向后Euler公式模型方程的向后Euler公式为：yn1

yn

h

yn1因此，1nyn1y

1

h=h11

h即，

E(h)

1绝对稳定域为：1

h

1绝对稳定区间：

h

Re()

00

h

,

0

R即对任意积分步长，向后Euler公式稳定（无条件稳定）2．

形

（P296，

己

）同样地，

0

R

其绝对稳定区间为：

o

h

,

0

Rh

Re(h)

0定义7

，如果数值方法的绝对稳定域包含了那么称此方法是A-稳定的。由定义可知，A-隐定的方法对步长h

没有限制。小结：一般而言，隐式公式的稳定性优于显式公式，对单步法或多步法都是成立的。这是选择数值积分公式时应当考虑的因一。§

.

线性与步法前面己比较详细地介绍了多种单步法：显式、隐式。他们都有广泛的应用，主要的优缺点是：优点：1、自启动，或称自开始。只要给出初值y0按差分公式即可一步一步递推地计算。逐次地求出y1,y22、易于改变步长缺点：1、低阶公式精度低，而高阶公式计算量大。上节课曾 ，以计算

f

(x,

y)

函数值次数作为计算一个衡量标准时，5阶R-K

有6次计算量，6阶R-K公式有7次计算量。为了克服单步法存在的上述缺点，希望能构造出精度高，计算量小的误差公式。它的基本思想是：“求yn+1时，充分利用如下信息：y

,

y

,ynn n

1

kf

n

,

f

n

1

,以构造线性多步积分f

n

公k

式，其中，k>0，即为k步法。”线性多步公式的构造方法同样有三种：基于数值微分、数值积分以及Taylor

展开的方法。本书仅介绍基于Taylor

展开方法，下面着重介绍基于Taylor展开法。对基于数值积分构造方法略做介绍。9.5.1

线性多步法的一般公式概括地说，基于Taylor展开构造线性多步公式是一种“待定常数法”。与R-K公式的构造原理、方法完全相同。可分两步：第一步是设定带若干待定常数的线性多步公式的一般形式，第二步可由Taylor展开定出特定常数，具体实现如下。一、一般线性多步法公式的约定时，利用计算ynk,

ynfnfnk

,ynk

1,fnk

1,ynk

2

,fnk

2,其中，yni

y(xni),fni

f

(xni

,ni

)xni

xn

ih做线性组合可构造出如下带待定常数的一般线性多步公式：(P360,

5.1)k

1

kyn

k

i

yn

i

h

i

fn

ii

0

i

0k

0

显式

k步法

0

隐式k步法二、由Taylor展开定义首先介绍线性多步法局部截断误差及阶的定义。y

f

(x,

y)定义8，y(x)设是初值问题的准确解，线性多步上的局部截断误差为：y(x0

)

y0公式（5.1）在xnkTnkk