大二下-數(shù)值分析

數(shù)值分析計(jì)算機(jī)科學(xué)與技術(shù)系所EDAcaiy

:

62785564東主摟8區(qū)407第九章

常微分方程初值問題數(shù)值解法§9.1

引言一、背景常微分方程有著廣泛的應(yīng)用，比如：火箭、的飛行控制；機(jī)器人控制；某些化學(xué)反應(yīng)過程的描述與控制。又比如：對一個(gè)復(fù)雜的集成電路，可根據(jù)定律列出一組關(guān)于節(jié)點(diǎn)電壓及支路電流的常微分方程。求解這組常微初值問題即可詳細(xì)分析電路的瞬態(tài)特性。但實(shí)際當(dāng)中，只有極少數(shù)常微分方程具有解析解，大多需借助數(shù)值方法。常微分方程分兩大類：初值問題和邊值問題，本章僅 前者。二、常微初值問題在初值問題里，集中 由單個(gè)的一階方程構(gòu)成的初值問題：(1.1)(1.2)

y

'(x)(

f

(x,

y(x))

f

(x,

y)

y(x

)

y

o

o定義：存在常數(shù)L>0，使x

[a,

b],

y1,

y2

R,有：

f

(x,

y1)

f

(x,

y2

)

L

y1

y2

,

稱

f

滿足對y的Lipschitz（）條件，L稱Lipschitz常數(shù)。f

滿足條件還可用更強(qiáng)的條件保證。如

f

連續(xù)有界f對y滿

y足Lipshitz條件。這概念在第七章非線性方程求解中提到過。明確一下本章對象：a

x

b,定理1：D

x,

y

y

R連續(xù)；（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：設(shè)f在區(qū)域D（如定理1所定義）上連續(xù)，且關(guān)于y滿足條件，設(shè)初值問題：y’(x)=f(x,y)，y(x0)=s的解為y(x,

s)，則：

f

(x,

s1)

f

(x,

s2

)

eL

x

x0

s1

s2該定理表明解對初值的敏感性，當(dāng)f

的L比較小時(shí)，不敏感；L大時(shí)為

問題?！?/p>

…在帶形區(qū)域D內(nèi)存在唯一解，或稱積分曲線y(x)，滿足微分方程及初值條件。三、數(shù)值解法基本思想數(shù)值方法的思想是將問題的求解要求降低為只求一系列離散點(diǎn)的函數(shù)值。具體地說，是將積分區(qū)間離散化，并對離散點(diǎn)函數(shù)值的近似值：x0

x1

xn

xn1

yn1y0

y1yn的差分方程。所謂差分方程構(gòu)成形如：yn1

G(yn1,ynynk)是“若干鄰近節(jié)點(diǎn)函值應(yīng)滿足的關(guān)系（方程）”?；蛘哒f是

“節(jié)點(diǎn)函數(shù)值的某種遞推關(guān)系”。比如從y(x0

)

y0

及其它可,

yk提供的初值

y1,

y2

,

出發(fā)，經(jīng)差分方程步進(jìn)地求出

yn1

。yn

yn

1

G

(

yn

1

,

ynk

)概括地說是“將連續(xù)問題離散化”，差分方程可分類如下：G

0k

0單步法

不含yn

1

顯式多步法yn

1

構(gòu)造差分方程的基本方法有三種：1、基于數(shù)值微分的方法（9.2.1）2、基于數(shù)值積分的方法（9.2.2，9.3，9.5.2）3、基于Taylor展開的方法

（9.3，9.5）重要介紹方法2、3的全部差分方程，且著重于方法3。因?yàn)樵瓌t上，方法3可以導(dǎo)出方法1與方法2的全差方，但方法1、2

有直觀性，將比較全面地加以介紹。第九章課后練習(xí)(20)作業(yè)：九3、4、6、7預(yù)習(xí)：9.4、9.5§9.2簡單的數(shù)值方法本節(jié)介紹四種簡單的差分方程構(gòu)造方法，事實(shí)上，它們可以由以上三種構(gòu)造方法的任一種推導(dǎo)得到。但為了大家對差分方程的構(gòu)造方法有較全面、完整的了解，第1小節(jié)的向前Euler方法與后退Euler方法重點(diǎn)從基于數(shù)值微分的方法去推導(dǎo)。9.2.1

（向前）Euler與后退（向后）Euler法一、基于數(shù)值微分推導(dǎo)向前Euler公式基于數(shù)值微分構(gòu)造差分方程，有兩個(gè)重要事實(shí)要經(jīng)常用到：（1）積分曲線上任—點(diǎn)(x,y(x))，滿足微分方程，即：y(x)

f

(x,

y(x)),x

[a,

b]（2）節(jié)點(diǎn)導(dǎo)數(shù)節(jié)點(diǎn)數(shù)值導(dǎo)數(shù)下面推導(dǎo)Euler向前公式：對初值問題，起始點(diǎn)x0的函數(shù)值y(x0)已知，那么對一個(gè)選定的h：x1

x0

h步長 ，如何構(gòu)造聯(lián)系x0與x1兩個(gè)節(jié)點(diǎn)函數(shù)值的關(guān)系式，即關(guān)于節(jié)點(diǎn)函數(shù)值的差分方程呢？積分曲線在x0滿足的微分方程，即以x0代入微方程：y(xo

)

f

(xo

,

y(xo

))左端是x0的節(jié)點(diǎn)導(dǎo)數(shù)?？蛇x用最簡單的數(shù)值微分公式：P131，第5行的兩點(diǎn)公式近似它：0

10(x

)

1

[

x

x

]

h

h

2只要將其中的符號作相應(yīng)修改：

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()若設(shè)：y(x)C[a,b]（連續(xù)）且積分步長h足夠小，可將余2略去，則有：y(x1

)

y(x0

)

hf

(x0

,y(x0

))這樣，即建立起y(x0)與y(x1

)的遞推關(guān)系式，由已知初始值y(x0)即可經(jīng)由差分方程導(dǎo)出

y(x1

),

y(x2

)

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

)

y(xn

)

hf

(xn

,

y(xn

))不過仔細(xì)檢查上面的推導(dǎo)過程會(huì)發(fā)現(xiàn)，在符號使用不上夠嚴(yán)格，這是因?yàn)榍懊嬉鸭s定：y(x)——積分曲線，是方程的準(zhǔn)確解，那么，y(xn

)——

y(x)

在xn

的準(zhǔn)確解。但是，上述遞推關(guān)系式是積分步長足夠小，以致數(shù)值導(dǎo)數(shù)誤差可以忽略的條件下得到的。因此，即使

y(x0

)

完全準(zhǔn)確，而經(jīng)過遞推公式計(jì)算的號使用上的

，y(x1

)

一般就不再準(zhǔn)確了。為了消除這種符引入以下記號加以區(qū)別：yn——差分方程在xn

的準(zhǔn)確值。上述遞推 系應(yīng) 寫為：yn1

yn

hf

(xn

,

yn

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

G(yn

)的差分方程。由于yn+1僅依賴于前面一個(gè)函數(shù)值yn，故稱單步法。而且，一旦有了yn，yn+1就可算出，又稱顯式公式，因此，向前Euler公式為：顯式、單步的積分公式。小結(jié)：基于數(shù)值微分構(gòu)造求積公式可概括為：“將微分方程離散化，以數(shù)值導(dǎo)數(shù)近似節(jié)點(diǎn)導(dǎo)數(shù)，從而構(gòu)造差分方程。”二、基于數(shù)值微分推導(dǎo)向后Euler公式同樣是基于數(shù)值微分方法推導(dǎo)，與向前Euler構(gòu)造相似，y(x0)已知，對選定的離散化步長h，有x1

x0

h

，要建立聯(lián)系x1，x0函數(shù)值關(guān)系的差分方程。積分方程在x1滿足的微分方程：y(x1

)

f

(x1

,

y(x1

))左端的一階節(jié)點(diǎn)導(dǎo)數(shù)，可選用P131第6行的兩點(diǎn)公式：1

1oh

2f

(x

)

1

[

f

(x

)

f

(x

)]

h

f

()類似地，將符號做相應(yīng)改變，即

f

y

，有：1

1o

1h

2y(x

)

1

[

y(x

)

y(x

)]

h

y

(

)

f

(x,

y(x

))略去數(shù)值導(dǎo)數(shù)的誤差，又用y1~

y(x1

)，有：y1

yo

hf

(x1

,

y1)得到的非線性代數(shù)方程。求出y1后，又可求出y2,y3,……一般地有：yn1

yn

hf

(xn1,

yn1

)——向后Euler公式向后Eluer公式為：隱式、單步法。小結(jié)：向前與向后Euler的推導(dǎo)區(qū)別僅在于：y

'(xn1

)

f

(xn1

,

y(xn1

))y

'(xn)

f

(xn,

y(xn))——向后Euler——向前Euler即向后Euler是 微分方程在離散點(diǎn)：xn+1的關(guān)系。而向前Euler是微分方程在離散點(diǎn)：xn

的關(guān)系。三、向后Euler公式的計(jì)算向后Euler公式是一個(gè)隱式公式，就是說一般從yn計(jì)算yn+1時(shí)不能像顯式公式那樣直接計(jì)算得到，而需求解一個(gè)關(guān)于yn+1的非線性代數(shù)方程。第七章已經(jīng)初步介紹求解單變量非線性代數(shù)方程的方法，著重在于通過不動(dòng)點(diǎn)迭代求解?？吹?，向后Euler公式本身恰好是以某個(gè)函數(shù)的不動(dòng)點(diǎn)形式表現(xiàn)的，因?yàn)椋簓n1

yn

hf

(xn1,

yn1

)

(

yn1

)即yn1

(yn1

)的不動(dòng)點(diǎn)表達(dá)形式，并可立即得到相應(yīng)的不動(dòng)點(diǎn)迭代格式：(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，

設(shè) 滿足以下兩個(gè)條件：復(fù)習(xí)第七章定理有1、對任意

x

[a,

b]

有

a

(x)

b2、存在正常數(shù)L<1，使對任意x,y

[a,b]，(x)

(

y)

L

(x

y)*則

(x)在[a,

b]上存在唯一不動(dòng)點(diǎn)x*

。定理2，設(shè)(x)

C[a,b]滿足定理1兩個(gè)條件，則對任意初值

xo

[a,

b]

，由

xk

1

(xk

)

得到的迭代序列

xk

收斂到

(x)

的Lk(2.5)

。不動(dòng)點(diǎn)x*，并有誤差估計(jì)：xk

x

1

L

x1

x0,

y

(

k

)

)y

(

k

1)n

1

yn

hf

(

xn

1設(shè)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

時(shí)，…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

梯形方法一、基于數(shù)值積分的推導(dǎo)：積分：對微分方程y(x)

f

(x,

y(x))，在區(qū)間[xn

,xn1

]xxxnxnf

(t,

y(t))dtn1n1y(x)dx

n1ny(xx)

y(

)

xn1f

(t,

y(t))dtxn這里要把積分作近似處理：

ft可利用遞形公式計(jì)算右端積分即有：2n1

nn

n

n1n1y(x

)

y(x

)

h

[

f

(x

,

y(x

))

f

(x

,

y(x

))]考慮到數(shù)值積分不可避免地存在誤差，可轉(zhuǎn)化為差分方程形式：2n1

n

n

nn1

,

yn1

)]y

y

h

[

f

(x

,

y

)

f

(x——梯形公式梯形公式為：隱式、單步法。注意：積分區(qū)間的取法以及積分 的數(shù)值積分處理方法不同，可得到更一般的基于數(shù)值積分方法推導(dǎo)的差分格式。比如：仍取積分區(qū)間

[xn,

xn1]

，但積分

以左端點(diǎn)的左矩形公式近似，則有：向前Eulery(xn1

)

y(xn

)

hf

(xn

,

y(xn

))

yn1

yn

hf

(xn

,

yn

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

yn

hf

(xn1,

yn1

)t二、梯形公式的計(jì)算P283，式（2.8）已構(gòu)造了一個(gè)不動(dòng)點(diǎn)迭代格式，也做了收斂性分析。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

改進(jìn)的Euler公式

y梯形公式：n1n

n

n

n1

n1

y

h

[

f(x

, )

(x

,

)]2是一個(gè)隱式公式，前面已經(jīng)，給出初值y(0)后不動(dòng)點(diǎn)迭代n12nn

n1n1n

n1y(k

1)

y

h

[

f

(x

,

y

)

f

(x

,

y(

k

)

)]給出序列

(

k

)yn1

只要h足夠小，有：y(

k

)n1n1k

yyn1

是梯形差分公式精確解。但上述迭代過程一般是個(gè)計(jì)算量很大的過程，而且從另一角度看，即使準(zhǔn)確的

yn1

，也僅是y(xn1

)

的近似值

yn1

y(xn1

)

。因此，耗費(fèi)大量機(jī)時(shí)去獲得一個(gè)不一定很準(zhǔn)的函數(shù)近似值并不值得，常常采用如下計(jì)算過程：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即僅作少量的幾步迭代，稱為“預(yù)報(bào)-校正”格式，甚至往往更加簡單，只取k=0，即：n1

y(0)nn

n

y

hf

(x

,

y

)2(0)n1

n1h

yn1

yn [

f

(xn

,

yn

)

f

(x

,

y

)]——改進(jìn)Euler公式（P283

（2.13））也可改寫為一個(gè)的顯式單步公式：2n

n

n

nn

nyn1

yn

h

[

f

(x

,

y

)

f

(x

h

,

y

hf

(x

,

y

))]還可以改寫成與R-K公式更相近的表達(dá)形式：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

單步法局部截?cái)嗾`差與階單步法的一般形式：yn1

yn

h(xn

,

yn

,

yn1,

h)(2.11)

h

(xn

,

yn

,

h)顯示單步法：yn1

yn與f(x,y)有關(guān)，比如：(xn

,

yn

,

h)

f

(xn

,

yn

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

n

n1

n

nn1

n1(x

,

y

,

y

,

h)

1

[

f

(x

,

y

)

f

(x

,

y

)]下面引出局部截?cái)嗾`差概念，用于判斷一個(gè)差分公式的精度，是常微數(shù)值解的一

十

基

的 念。Tn1

y(xn1

)

y(xn

)

h(xn

,

y(xn

),

h)定義1，設(shè)y(x)是初值問題（1.1），（1.2）的精確解，稱為顯式單步法：yn1

yn

h(xn

,yn

,h)的局部截?cái)嗾`差。1、向前Euler局部截?cái)嗾`差yn1

yn

hf

(xn

,

yn

)yn1

y(xn

)

hf

(xn

,

y(xn

))向前Euler：若令：則不難看出，二者是有區(qū)別的：yn+1表示從y(x0

)出發(fā)應(yīng)用向前Euler差分格式，從yn求出yn+1。一般來：yn

y(xn

)

。而

yn1表示從xn求yn+1

時(shí)，使用了積分曲線在xn的準(zhǔn)確值y(xn)，并應(yīng)用向前Euler差分方程得到的xn+1的函數(shù)近似值

yn1

。那么，局部截?cái)嗾`差正是

：Tn1

y(xn1

)

yn1

y(xn1

)

y(xn

)

hf

(xn

,

y(xn

))應(yīng)如何比較右端各量，以便將右端簡化呢？基本方法是將各量在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局部截?cái)嗾`差首P=1，一階顯示公式定義2，設(shè)

y(x)

是（

.

）（.

）的準(zhǔn)確解，若存在最大整數(shù)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稱：為局部截?cái)嗾`差主首。2、向后Euler局部截?cái)嗾`差向后Euler公式：yn1

yn

hf

(xn1,

yn1

)yn1

y(xn

)

hf

(xn1,

yn1

)記：表示，在xn時(shí)，取y(xn)為準(zhǔn)確值遞推得到的xn+1的差分方程。與推導(dǎo)向前Euler公式局部截?cái)嗾`差相似，應(yīng)?。簲?shù)值yn1Tn1

y(xn1

)

yn1

y(xn1

)

y(xn

)

hf

(xn1

,

yn1

)作為向后Euler局部截?cái)嗾`差。一般來說但書上（P285）為了推導(dǎo)局部截?cái)嗾`差主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、

形公式局部截?cái)嗾`差2,

y

)]n1

nn

nn1

n1y

y

h

[

f

(x

,

y

)

f(x梯形公式：主的方法定義梯形P285仍取便于推導(dǎo)局部截?cái)嗾`差首公式局部截?cái)嗾`差：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方法本節(jié)介紹基于Tayler展開方法推導(dǎo)Runge-Kutta公式為了便于理解，從簡單的2階公式的構(gòu)造入于，因此，1、2兩小節(jié)對—

。9.3.2

二階顯式R-K方法一、一般二階公式約定梯形公式的推導(dǎo)， 是通過與微分方程等價(jià)的積分方程：nxf

(t,

y(t))dtxn1y(xn1

)

y(xn

)

還利用積分方程形式去理解如何約定一般給出的。下面的二階R-K公式。通過等價(jià)積分方程構(gòu)造差分格式，關(guān)鍵是用不同方法對積分做數(shù)值積分：ynk1通過兩種簡單的顯式單步公式去理解一般二階R-K公式。1、向前Euler公式n

nyn1

yn

hf

(x

,

y

)ynk1令：

k1

f(xn

,

yn

)yn1

yn

hk1則：K1是函數(shù)f(x,y)在左端點(diǎn)xn的函數(shù)近似值。向前Euler可看作用左矩形公式處理積分 得到的最簡單的差分公式。2、改進(jìn)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)在左端點(diǎn)xn的函數(shù)近似值，以紅點(diǎn)為標(biāo)識；而K2則可看作f(x,y)在右端點(diǎn)xn+1的函數(shù)近似值。因此，改進(jìn)Euler公式可以看作用梯形積分方法近似處理積分而左、右點(diǎn)函數(shù)近似值為K1與K2所給出。將上述數(shù)值積分處理積分方程中積分 的思想一般化可在以下兩方面加以推廣，從而給出一般二階公式的約定：

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與改進(jìn)Euler公式相比，所作兩點(diǎn)擴(kuò)展表現(xiàn)在：（1）K2的擴(kuò)展從圖中看出：22

h

(xn

,

xn1

],

0

1K2可以取

(xn

,

xn1

]

中任一點(diǎn)：

xn（2）K1與K2的線性組合，（或理解為數(shù)值積分公式的擴(kuò)展，不僅局限于左矩陣或梯形積分模式）。引入線性組合系數(shù)C1與C2，在一般二階R-K公式中含2

,21

,C1,C24個(gè)待定常數(shù)。二、由Tayler展開定

2

,21

,C1,C2目標(biāo)是適當(dāng)選擇四個(gè)待定參數(shù)，使yn+1在xn的Taylor展開，具有與y(xn+1)在xn的直接Taylor展開有盡可能高的相同的階，也就是使yn+1的階數(shù)盡可能高。為 推導(dǎo)時(shí)

寫 ，令：yn

y(xn

)，fn

f

(xn

,

y(xn

))

f

(xn

,

yn

)局部截?cái)嗾`差可以寫為：nnTn1

(xn1

)

h[C1

(xn

,

n

)

C2(xn

2h,

21h

n

)]f(x,y)作二元函數(shù)在局部截?cái)嗾`差表達(dá)式中，先看最后一Taylor展開：(xn

2h,

h

n

h

)

'

x

'

xn

h

O

h2n

21

n

n

x

n

n

2

y

21

n三那么，局部截?cái)嗾`差展開式前出現(xiàn)

fn，fx

'

h，fy

'

fnh三

應(yīng)表達(dá)為：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表達(dá)式，依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階精度，應(yīng)有：1

C1

C2

01

C

0

2（3.9）2

211

C

02三個(gè)方程定4個(gè)未知數(shù)，解不唯一，可以得到一族二階R-K公式，比如：改進(jìn)Euler公式是當(dāng)：2

21

11

22C

C

1—

二

R-K公式，類似可以構(gòu)造許多其它二階R-K公式。.

.

顯式Runge-Kutta

一回到微分方程的等價(jià)積分方程形式：nxf

(t,y(t))dtxn1y(xn1

)

y(xn

)

構(gòu)造相應(yīng)差分格式：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公式一般形式的構(gòu)造思想：…通過f(x,y)在r個(gè)點(diǎn)上函數(shù)近似值做數(shù)值積分，而數(shù)值積分公式的系數(shù)是待定的：C1，……，Cr。括，用數(shù)值積分近似積分xnf(t,

y(t))dtxn11

1

2

2

r

rh[C

k

C

K

C

K

]其中，兩個(gè)因素是待定的：i1j

1yn

h

ijk

jxn

ih，（1）ki

fi即所取作數(shù)值積分的函數(shù)值為待定。（2）C1，……，Cr，即數(shù)值積分公式待定。通過引入大量待定常數(shù)，使之達(dá)到盡可能高的階數(shù)。類似二階顯式R-K公式的構(gòu)造與推導(dǎo)，可以構(gòu)造高階R-K公式。9.3.3

三階與四階R-K公式（了解公式即可）從P288~289看到，在使用三階、四階R-K公式時(shí)，

一所需計(jì)算

f(x,

y)

函數(shù)值的次數(shù)分別為3次及4次，下面有關(guān)于R-K公式階數(shù)與函數(shù)值計(jì)算的表：R-K階數(shù)≤45計(jì)算次數(shù)≤4697…….因此，一般取R-K公式≤6階。第9

章課后練習(xí)(21)作業(yè)：九8、11、12§

.

單步法的收斂性與穩(wěn)定性9.4.1

收斂性與相容性一、單步法的收斂性先通過一個(gè)簡單例子，直觀了解一下收斂性的概念。例：用逐次二分積分步長的向前Euler法求解。

y

'

y

0

x

T

y(0)

y0其準(zhǔn)確解為：y

y(x)

e

xo相應(yīng)圖形如下：積分區(qū)間內(nèi)任一固定點(diǎn)C：0≤c≤Tc那么，積分曲線在這點(diǎn)的值為：y(c)

y0e下面 使用逐次二分積分步長的向前Euler公式積分上述初值問題時(shí)，當(dāng)積分步長越來越小時(shí)，c點(diǎn)處函數(shù)近似值的變化情況：（1）將[0,c]分作n等分，則積分步長n

nh

c

0

c（2）取—定n值，從x0出發(fā)，用向前Euler一步一步積到c，可得C點(diǎn)近似的函數(shù)值yn。yn1

yn

hf

(xn

,

yn

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

yn,n

2

，k為整數(shù)。kn

on

k

lim

y

y(c)

y

ec從圖上看出題收斂。稱向前Euler式積分上述初值問強(qiáng)調(diào)一點(diǎn)：在收斂性概念里，c是積分區(qū)間任一相對固定的：

x0

c

T

，它不隨積分步數(shù)或積分步長的變化而變化。這點(diǎn)十分重要，下面給出收斂性定義。n

,

h)定義3，若一種數(shù)值方法（如單步法

yn1

yn

h(xn

,

）,

當(dāng)

k

時(shí)，yn

y(c)

稱該方對于固定的xn

xo

nh

c法是收斂的。這里對書上定義作兩點(diǎn)修改和說明，以避

理

上。

nh（1）xn

x0改為xn

xo

nh

c

，強(qiáng)調(diào)是固定點(diǎn)c。（2）ynn

y(x

)h0

n改為

n

h0

，會(huì)更確切些。y

n

y(c)書上的表述不很嚴(yán)格，h~n是存在關(guān)系的，當(dāng)h0時(shí)，應(yīng)

有

n

。但一方面

n

，另一方面極限結(jié)果仍含n，不夠嚴(yán)格。下面顯式單步法的收斂性。有P階精度，且增定理3，設(shè)顯式單步法yn1

yn

hxn

,

yn

,

h量函數(shù)，關(guān)于y滿足Lipshitz條件：(x,

y,

h)

(x,

y,

h)

L

y

y

y(x0

)又設(shè)初值y0是準(zhǔn)確的，即

yo

，由其整體截?cái)嗾`差為：py(xn

)

yn

O(h

)一下。在證明前，將一些相關(guān)概念及定理結(jié)論、條件1、局部截?cái)嗾`差與整體截?cái)嗾`差有幾個(gè)相關(guān)符號再明確一下：y(xn1

)——積分曲線

y(x)

在

xn1

準(zhǔn)確值。yn1——

yn1

y(xn

)

h(xn

,

y(xn

),

h)，即設(shè)：yn

y(xn

)準(zhǔn)確，顯式單步法積分一步至的準(zhǔn)確值。xn1，從x0出發(fā)，顯式單

h

xn

,

yn

,

h——

yn1

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

的準(zhǔn)確值。這里注意：yn1

與yn1

的不同含義，因此，——局部截?cái)嗾`差，前進(jìn)一步誤差。y(xn1

)

yn1en1

y(xn1

)

yn1——整體截?cái)嗾`差，各步誤差的積累?？梢詮膱D上理解這兩種截?cái)嗾`差的含義以及它們之間的關(guān)系：y=Tn+1y(xn)ynn12、證

思再把定理?xiàng)l件與結(jié)論強(qiáng)調(diào)一下：條件：（a）

y(xn1

)

yn1

O(hp1

)，P階精度。，即增量函數(shù)關(guān)于y（b）(x,

y,

h)

(x,

y,

h)L

y

y滿足Lipschitz條件。（c）

yo

y(x0

)

，即y0準(zhǔn)確。即整體截?cái)嗾`差不超過h的P階小y(x

)

y

O(hp

)n1

n1量。結(jié)論：概括地說，定理給出局部截?cái)嗾`差的條件以及其它有關(guān)條件，而結(jié)論是要求整體截?cái)嗾`差應(yīng)當(dāng)滿足的關(guān)系。因此，應(yīng)當(dāng)將整體截?cái)嗾`差與局部截?cái)嗾`差聯(lián)系起來：y(xn1

)

yn1

y(xn1

)

yn1

yn1

yn1

y(xn1

)

yn1

yn1

yn1其中，第一 為局部截?cái)嗾`差，其誤差界是清楚的，問題在第二 的估計(jì)。3、證明分別對一、二兩 估計(jì)：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使用前面引入的整體截?cái)嗾`差記號，有n1e

(1

hL

)

en

chp1是整全截?cái)嗾`差的遞推公式，反復(fù)遞推：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推導(dǎo)：(1

hL

)n

(ehL

)n

eTL其中，T為積分區(qū)間總長度。因此，整體截?cái)嗾`差上界為：n0Le

e eTL

ch

p

(eTL1)nnnLe

y(x

)

y

ch

p

(eTL1)

O(h

p

)如何應(yīng)用定理1

一個(gè)顯式單步法積分一個(gè)具體問題的收斂性，必須對方法的階、以及增量函數(shù)是否對y滿足Lipschiz條件加以分析。P292以改進(jìn)Euler

法為例對其改斂性作了分析，希望大家掌握它，再給一個(gè)思考題：思考：

y

'

f

(x,

y)

y(x

)

y

0設(shè)f

(x,y)對y滿足Lipschiz條件，y=y(xo

)準(zhǔn)確,證明：用P289（3.

）公式做數(shù)值積分時(shí)收斂。二、相容性定義4，若單步法yn1

yn

h(xn

,yn

,h)的增量函數(shù)滿足：(x,y,0)

f

(x,y)

，則稱單步法與初值問題：

y

'(x)

f

(x,y)

相容

y(x

)

y

o

oP291從顯式單步法局部截?cái)嗾`差關(guān)于積分步長h的Taylor展開引出相容性概念。

也可以從另一個(gè)角度理解相容性含義。設(shè)顯式單步法：

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

時(shí)，顯式單步法格式趨于微分方程本身。這是相容性的直觀意義。（1.1）和（1.2）相容的充分必要條件是P1。定理4，p

階方法

yn1

yn

h

xn

,

yn

,

h

（4.1）與初值問題9.4.2

絕對穩(wěn)定性與絕對穩(wěn)定域穩(wěn)定性是差分方程求解常微分方程另一個(gè)重要理論問題，它研究差分方程對誤差的性質(zhì)。它的含義是在積分步長相對固定的條件下，若在某一步引入誤差，在繼續(xù)遞推過程中，這個(gè)引入的誤差會(huì)如何？它的性質(zhì)與什么因素有關(guān)？如何防止和控制誤差的惡性增長？這些都是穩(wěn)定性所關(guān)心的問題。絕對穩(wěn)定性是穩(wěn)定性概念中比較常用、比較基本和經(jīng)典的。一、例

及通過一個(gè)例子，了解穩(wěn)定性概念和它的基本含義。例4

初值問題：y

100

yy(0)

0為了使問題更直接、更簡單，這里把例中初值1改為

o

，即初值本身就表現(xiàn)為引入的某種誤差。 分別后Euler

求解這初值問題時(shí)表現(xiàn)的不同現(xiàn)象。向前與向積分曲線，即準(zhǔn)確解是：y(x)

oe

100

x（1）向前Euleryn1

yn

h(100

yn

)

1.5yn因此，當(dāng)

yo

o

，h=0.025時(shí)y1

1.5oy2

2.25o己放大到P294

表9-4列出向前Euler積分至0.1時(shí)，初值的

05.0625

，即放大了5倍多。0h

0.0254-4（2）向后Euleryn1

yn

h(100

yn1

)

yn

2.5

yn113.5nyn1y

看P394：表9-4第3例數(shù)據(jù)，給出了向后Euler以0.025步長積分同一問題的數(shù)值結(jié)果，在圖上可以看出，初值時(shí)引入的誤差

0

，在隨后各步遞推計(jì)算過程中，逐漸而這于0，說明初始誤差受到控制。這個(gè)例子表現(xiàn)了向前Euler

與向后Euler不同的誤差

性質(zhì)，下面給出定義。法是穩(wěn)定的。用符號表示為定義5，

若一種數(shù)值方法在節(jié)點(diǎn)值yn上有大小為

的擾動(dòng)，于以后各節(jié)點(diǎn)值ym

(m

n)

上產(chǎn)生的偏差均不超過

,則稱該方m

n（m

n）有很大，通常要簡化對象，穩(wěn)定性問題的一般性并稱之為模型問題。二、模型方程y

y,

C,稱

Re()

0

為模型方程。在模型方程中，假設(shè)

Re()

0

是為了使模型問題y

yy(xo

)

yo穩(wěn)定準(zhǔn) 為：

(x)

(

x

)0e

x0模型問題的初值受到擾動(dòng)時(shí)，對模型問題的解產(chǎn)生影響：如果初值擾動(dòng)產(chǎn)生的影響能逐漸是穩(wěn)定的，反之是不穩(wěn)定的。下面說明Re

()

0

保證模型問題穩(wěn)定。設(shè)，有一初值受到擾動(dòng)的模型問題：y

yy(xo

)

yo

o是與前面給出的模型問題完全相同的系統(tǒng)，區(qū)別僅在于后，其準(zhǔn)確解為：者在初值帶有誤差0y

(x)

(

yo

(

x

x

)

o

)e

o那么：e

(

x

xo

)(x)

(x)

oeRe(

)(

x

xo

)

oy(x)

y(x)

當(dāng)：Re()

oRe()

o穩(wěn)定不穩(wěn)定y(x)

y(x)

三、顯式單步法公式穩(wěn)定性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

公式對誤差或擾動(dòng)的遞推關(guān)系式，而且可以看出，誤差 的關(guān)系式與原來的差分格式完全相同。若要求誤差不增長，即1

h

1

n1

yn1

n

yn稱

1

h

1

為向前Euler的絕對穩(wěn)定域。=h平面定義6，單步法

yn1

yn

h(xn

,

yn

,

h)

用于的模型方程y

y

(Re()

o)

1

，則稱單步法是絕對穩(wěn)定稱為絕對穩(wěn)定域，它與實(shí)軸的滿足

E

hE(h)

1若得到的解yn1

E(h)yn的。在h

平面上，區(qū)域交稱絕對穩(wěn)定區(qū)間。所以，對向Euler公式，絕對穩(wěn)定區(qū)間是：2

h

Re()

0

0

h

2

o

R回到例4

100,

2

100h

0應(yīng)有

0

h

2

0.02100是其絕對穩(wěn)定區(qū)間。因此，當(dāng)取0.025時(shí)是不穩(wěn)定的。2、二階R-K公式（

.

）關(guān)鍵是推導(dǎo)模型方程下的二階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因此，絕對穩(wěn)定區(qū)域?yàn)椋?h

21

h

1三、四階R-K方法自己推導(dǎo)一下相應(yīng)的表達(dá)式四、隱式單步法穩(wěn)定性1．向后Euler公式模型方程的向后Euler公式為：yn1

yn

h

yn1因此，1nyn1y

1

h=h11

h即，

E(h)

1絕對穩(wěn)定域?yàn)椋?

h

1絕對穩(wěn)定區(qū)間：

h

Re()

00

h

,

0

R即對任意積分步長，向后Euler公式穩(wěn)定（無條件穩(wěn)定）2．

形

（P296，

己

）同樣地，

0

R

其絕對穩(wěn)定區(qū)間為：

o

h

,

0

Rh

Re(h)

0定義7

，如果數(shù)值方法的絕對穩(wěn)定域包含了那么稱此方法是A-穩(wěn)定的。由定義可知，A-隱定的方法對步長h

沒有限制。小結(jié)：一般而言，隱式公式的穩(wěn)定性優(yōu)于顯式公式，對單步法或多步法都是成立的。這是選擇數(shù)值積分公式時(shí)應(yīng)當(dāng)考慮的因一?！?/p>

.

線性與步法前面己比較詳細(xì)地介紹了多種單步法：顯式、隱式。他們都有廣泛的應(yīng)用，主要的優(yōu)缺點(diǎn)是：優(yōu)點(diǎn)：1、自啟動(dòng)，或稱自開始。只要給出初值y0按差分公式即可一步一步遞推地計(jì)算。逐次地求出y1,y22、易于改變步長缺點(diǎn)：1、低階公式精度低，而高階公式計(jì)算量大。上節(jié)課曾 ，以計(jì)算

f

(x,

y)

函數(shù)值次數(shù)作為計(jì)算一個(gè)衡量標(biāo)準(zhǔn)時(shí)，5階R-K

有6次計(jì)算量，6階R-K公式有7次計(jì)算量。為了克服單步法存在的上述缺點(diǎn)，希望能構(gòu)造出精度高，計(jì)算量小的誤差公式。它的基本思想是：“求yn+1時(shí)，充分利用如下信息：y

,

y

,ynn n

1

kf

n

,

f

n

1

,以構(gòu)造線性多步積分f

n

公k

式，其中，k>0，即為k步法?！本€性多步公式的構(gòu)造方法同樣有三種：基于數(shù)值微分、數(shù)值積分以及Taylor

展開的方法。本書僅介紹基于Taylor

展開方法，下面著重介紹基于Taylor展開法。對基于數(shù)值積分構(gòu)造方法略做介紹。9.5.1

線性多步法的一般公式概括地說，基于Taylor展開構(gòu)造線性多步公式是一種“待定常數(shù)法”。與R-K公式的構(gòu)造原理、方法完全相同?？煞謨刹剑旱谝徊绞窃O(shè)定帶若干待定常數(shù)的線性多步公式的一般形式，第二步可由Taylor展開定出特定常數(shù)，具體實(shí)現(xiàn)如下。一、一般線性多步法公式的約定時(shí)，利用計(jì)算ynk,

ynfnfnk

,ynk

1,fnk

1,ynk

2

,fnk

2,其中，yni

y(xni),fni

f

(xni

,ni

)xni

xn

ih做線性組合可構(gòu)造出如下帶待定常數(shù)的一般線性多步公式：(P360,

5.1)k

1

kyn

k

i

yn

i

h

i

fn

ii

0

i

0k

0

顯式

k步法

0

隱式k步法二、由Taylor展開定義首先介紹線性多步法局部截?cái)嗾`差及階的定義。y

f

(x,

y)定義8，y(x)設(shè)是初值問題的準(zhǔn)確解，線性多步上的局部截?cái)嗾`差為：y(x0

)

y0公式（5.1）在xnkTnkk