第九章§1y(x

f(x,)y

ax 方

xy'-2y=4x

y'=2y4

(x,y)=2yx

且給出初

y(1)=-就得到一階常微分方程的初值問題dy

f(x,

y)

2y y(1)只要函數(shù)

(x,

y)適當(dāng)光滑連續(xù)，且 李于于f(x,

y)

f(x,y)

y[a,b]，令a=x0<x1<…<xn=b,其中hk=xk+1-xk,如是等距節(jié)點(diǎn)h=(b-a)/n,h稱為步長(zhǎng)。即yk≈y(xk)，這樣y0,y1,...,yn稱為微分方程的數(shù)值解。微分方程離散化常用方A用差商代替微

y'dy

f(x,y)

y(xn1)

y(xn)

f(x

,y(x(xn,yn

h

xn

yn

y(xn

yn1

y(xn1

yn1

f(x,y yn1

yn

(xn

yn

n0,1,B.數(shù)值積 用數(shù)值積分方法離散化

y'

f(x,y)

dydx

f

用yn1

yn代替y(xn1

y(xn))n f(x,n

y)dx

h

(xn

ynyn1

(xn

yn

(n

C 在 附近y(x)的

展開h2y(xnh)y(xn)

hy'(xn)2y'(xn)h2y(xn

(xn

y(xn))

y'(xn)

f(x,

ax bby(x) h的線性部分

yn

y(xn

y(xn1

的近似值yn1

yn

(xn

yn

n §2（Eular一、計(jì)算公式

f(x,

a

b

yn1 y(x

)y(x)

1. x0

y0

yyx

fx,ydx(x0,y0

y0

y0

y

x

y(x

)f(x,y

(x0,y0

h

hy1

y0

(x0

y0

y

xn1點(diǎn)的yn1

yn

(xn

yn

n用分段的折線近函數(shù),此為 2、EulerTn1

y(xn1)

yn1其中yn1是當(dāng)

y(xn

y(xn1)在xn

h2y(xn1)y(xnh)

y(xn)hf(xn,y(xn))+

y'()xny(xn1)

y

h

y(xn)hf(xn,y(xn

hfx,y

yn1

y(xn)hfxn,y(xn)h2Tn1

y(xn1)

yn1

2

xn

令M

a

y'(x)

y(x)

h2

M

O23、總體方法誤差 充分光滑，或滿Lipschitz條件fxn

y(xn)

fxn

yn

y(xn)第

ny(xn1)

y(xn1)

以下估計(jì)

yn1

其中yn1

y(xn)h

(xn, 總體方法誤差yn1

y(xn)

xn

y(xn)

xn

yny(xn)

fxn

y(xn)

fxn

yn由Lipschitz條件

hL

y(xn)

hL) en1

y(xn1)yn1

y(xn1)yn1

yn1

yn1

1

hL

N

e0

y(x0)

1 1

hL

eN

hL

hL)N1由局部

O(h2

hL

hl)N1N k

O(h2

N11k(1hL)N

)Oh1hL1lim(1

lim(1

xNh

eL(xNx0

無關(guān)常數(shù)

Oh4、微分方程數(shù)值解的穩(wěn)定

n1

yn1

yy 定義：一種數(shù)值

y'=

其中

是復(fù)常數(shù)定的

h0若計(jì)算誤差

在計(jì)

ynkk

不產(chǎn)生增大的誤差，即nk

n，稱

與這種方法對(duì)

y'

yn1

Euler法的絕對(duì)穩(wěn)定區(qū)

y*h

- 從而 n1h當(dāng)1h

是絕對(duì)穩(wěn)定區(qū)如果整個(gè)左半平面是絕對(duì)穩(wěn)定區(qū)域稱A穩(wěn)定的。二、向后（后退的）Euler用向后差商：y'(xn1)

y(xn1)h

y(xn

y'dy

f(x,y)為隱式算法

yn1

yn

yn1y(x) 為避免解非線性方程，與Euler法結(jié)合y(0)

hfx,y y(k1)

hf(x

,y(k) k0,1,

向后Euler方法收斂條件與截?cái)嗾`y(0)

hfx,y y(k1)

hf(x

,y(k) k0,1,

y(ky

(kyy

(ky,y,

收斂條

O(h2

Oh向后Euler對(duì)y

則

1

[(1h)(1[1

)h2]0.5

只要Re(h0

因此向

法是A穩(wěn)定的。但收斂要0

仍受三、梯形公由積分途徑：y(xn1

y(xn)

f

y

dy

f(x,y)積分用梯形公式，且令：yn1

y(xn1),

y(xn則得 yh

f(x,y)f(x , )

同樣與Euler法結(jié)合，形成迭代算法，對(duì)

y(0)

hf(x,y y(k1)

h

f(x,

)f(x

,y(k)),

k

由迭代算法，對(duì)

y(0)

hf(x,y y(k1)

h

f(x,

)f(x

,y(k)),

k

y(ky

(k y y

fxn1

(kyy

fxn1

(kyy2

(kyy

(ky,y,

0hL2

0(h2

h

f(x,

)f(x , )

對(duì)yy

yh(y

h

h

1

1

Re(h)

2n 12

2當(dāng)Re(h)0 梯形公式是四、改進(jìn)的尤拉公 yn1

yn

(xn

yn

h[f(x,

)f(x ,

h[f(x,

)f(x h, hf(x,

yp

hf(xn

yn hf(x , y )/y例：用尤拉公式和改進(jìn)的尤拉公式解初值問y'

y2x

(0

x

hf(x,y

yp

yn1

ynhf(xn,yn

)/解：取步長(zhǎng)

0.， 尤拉公式為

h( 2xnyy

n h( 2xny y 改進(jìn)的尤拉公式為

h(yp

2xn1yp計(jì)算結(jié)果略

尤拉兩步公 y(xn1

y(xn1)替代yx

),nyn1n

(xn

yn稱其為尤拉兩步公式（首先要確定yn-1,yn兩個(gè)值)。yn1

(xn

yn校

yh[f(x,

)f(x ,

考慮上面給 －校正系統(tǒng)的局部截?cái)嗾` 公式：假

y(xn

yn1

y(xn1 其局部截?cái)嗾`差

y(xn1)

yn1

3y(xny(xn1)

y(xn

h)

y(xn)

(xn

y(xn))2!

'(xnhh

h4 (4)23!y'(xn)4!2

xn

2hf(x,y)y(x

)2hy'(x)y(xh)2hy'(x

'y(xn)''y(xn)'

h(xn)hh2h2(xn)2

h(xn)hh3h(xn)

h(xn)hh4h4!(xn)4!

y(4)(x) nny(4)(x)nn

(xn對(duì)校正公式：假定yn1是準(zhǔn)確的，即yn1y(xn1 其局部截?cái)嗾`差為

y(xn1)

yn112

(xn h[f(x,

)f(x ,

y(x)h[y'(x

h2 y]]y(x)hy(x)'2y(x)4y(x)nnnn

y(xn)

2[y(xn)

y(xn)

(xn)

y(xn

y(xn1)y(xn1)

yn1

14y(xn1)

yn1

555

(yn1

yn1y(xn1)

1(

n1

yn1利用估計(jì)出的誤差作為計(jì)算結(jié)果的一種補(bǔ)償,有能使精度得到改善,此時(shí)的計(jì)算公式

2hy'

yn1

yn1

(xn,yn

h[f(x,

)f(x ,

改 mn1

pn1

5(

cn

m'm

f

,mn1

y(xn1)

yn1

4(y

n1

yn1校

yh(y'm'

y(xn1)

1(

yn1

改

1(

y'y

f(xn1,

yn1 —yn1

y(xn1)

y(xn)其中yx)f(x,y),yx)fx,yf(x 但由于公式中各階偏導(dǎo)數(shù)計(jì)算復(fù)雜，不實(shí)用y'

f(0)

(0)

f(0)

f

f

(2)

(j

f

j2)

f

j2)

(j

j

如果將Euler公式與改進(jìn)Euler

yn1

yn f(x,y y y

yn

h(2

K1

2K2

f(x,y

hf(x,y

K2 f(xnh,

hK1

校 yh[f(x,

)f(x ,

以上兩組公式都使用函數(shù)

(x,

y)線性組合來計(jì)算y(xn1)的近似值Euler公式：每步計(jì)算一次

(x,

y)改進(jìn)Euler公式：需計(jì)算兩次

(x,

y)

pKK1

f(xn,yn

,p),p)Ki

f(xnaih,

hbijj

j

2,于是可考慮用函數(shù)

(x,

y)線性組合來構(gòu)造近似公式，構(gòu)造是要求近似公式(xn,yn)處的Taylor展開式與解y(x)在xn處的Taylor展二、二 －庫塔方一般地，RK方法設(shè)近似公式K1

f(xn

yn

Ki

f

aih,

hbijKj j

2, ,p)在(xn,yn)處的Taylor展開式與y(x)在xn處的Taylor展開的前面項(xiàng)盡可能多地重合yn1

yn

c2K2當(dāng)p時(shí)，近似公式為Kf(xy Kf (xna2h,Kf

hb21K1上yn1

ynyn

(xn

yn f'(x,y)f(x,

)]}

c2)

(xn

ync[af'(x,y) f'(x,y)f(x,

O(h y(xn1)在xn處的Taylor展開式 h2 y(xn1)

y(xn)

(xn)2y(xn)O(hyn

f(xn

ynh2 2[fx(xn,yn)fy(xn,yn)f(xn,yn)]O(hc1c2

c

1/

c

1/

(h3 取c1

近似公式

K1

ff(x

h, 取

1

yn1

ynhK2

f(x,y f(x , hK2) 三、三 －庫塔方

可導(dǎo)出三階RK常用的三階RK yh(K4KK Kf(x,y K K

(

,

2K1Kf(xh,yhK2hK 對(duì)p

四、四 －庫塔方常用的四階RK yh(K2K2KK Kf(x,y K2

f

2,

2K1KfKf

h,

hK K4f

h,

hK3例：設(shè)取步長(zhǎng)h=0.2,從x=0直到x=1用四 —庫方法求解初值問方法求解初值問題

yy

(0

x

解：由經(jīng)典的四 —庫塔公式

y

h

Kyy hyy

K

yy2xny

nKf(x,y h h

hK

2xn KfKf

h,

2K1

yn

2KfKf

h,

hK

Ky Ky

h

2xn yhy

f(xh,yhK

yn2K2

2(xnh)n3 n3

兩點(diǎn)說明RK方法的導(dǎo)出基于Taylor展開，故要求所求題的解具有較高的光滑度如果解的光滑性差，則用四階RK方法解的效不如改進(jìn)Euler五、變步長(zhǎng) —庫塔方 步長(zhǎng)求一近似值y(h)y(xn1)

(h)yy

ch5將步長(zhǎng)折半，即

h為步長(zhǎng)從x跨兩步到 ，求一近 值y(2)每跨一步的截?cái)嗾`差是c

h2

,y(xn1)

(2)yy

2c

h2y(x

)y(2) 由上兩

yy(xn1)y

(h) y(xn1)

(2) y y

(2)[[

(h) K1f(xn K3

f

h,y

hyn

h(yn h(yhK)yh1(h)2 n

h(y

)

h

1(h)314yn1

h6h

2K2

n

K4y1h1(h)21(h)31(h)4n 1h

1(h)21(h)31(h)4 n

絕對(duì)穩(wěn)定區(qū)域1h

1(h)21(h)31(h)4 21-3- -0--單步

§4.需重新計(jì)算多個(gè)點(diǎn)處如何通過較多地利用前面的已知信構(gòu)造高精度的算法計(jì)算yn1，這就是多步法的基本多步法中最常用的是線性多步法，它的計(jì)算公式中出現(xiàn)yn1，ynynr及

(xn1

yn1

f(xnr

ynr)的一 項(xiàng)，其一般形

ii0

hi1

fni其中i,i均為常數(shù)

f(xk

yk)

n1

0，隱式構(gòu)造線性多步公式常用Taylor展開和數(shù)值積分方開式相比較，要求它們前面的項(xiàng)重合，由此確定參數(shù)i,i以r為例：設(shè)初值問題的解y(x)充分光滑，待定兩步公式

ii0

i1

fni

0

1

0

fn1記y(k

y(k)(x

)

則y(x)在

處的Taylor展開2p 2py(x)

yn

y'(x

yxn)y2

(x

xn

(py (xyp!

xnnnO((xnn

x)p1假設(shè)前n步計(jì)算結(jié)果都是準(zhǔn)確的，即 y(xiy'(x)f(x,

)

n則

y( h)

y'h

h2 n

y(4)

2! y(5) h4

h5

O(h64! f(x ,

)y'(x )y'y''h nh2

y(4)

2!y(5) h3 h4

O(h5 f(x,y)y'

4!

f(x , )y'(x

)y'

y''h h2

y(4)

2!y(5) h3 h4

O(h5 4!將以上各公式代入并整理， ()

( )y'h(1

)y''h2

(

1)

1)y(4)h4 (

1

1)y(5)h5

O(h6 為使上式有p階精度，只須使其與y(xn1)在xnTayloryy (5)yyy(xn1)ynyhn n n

h5O(h62! 的前p+12 21111 a1 6

1i,i，使其滿足 1,0, 1,

yn1

yn

h(2

n1

fn又如：解上面方程組得

0,1, 1, 相應(yīng)的線性二步四階公式（Simpson公式

h(

4

其截?cái)嗾`差

h5y(5)

O(h6n由此可知，線性二步公式至多是四階公式n a1

0 2 211

1 a1

16 6 hphR1 O(顯然，線性（一 取r=3，并令

2

i1 i1i0 (i)kk(i)k

3i0

i1可解得

0

,

,2

,

相應(yīng)的線性多步公

ii0

i1

f

h(55

59

37

9

n2

n3

(5) R251h5y(5)

O(h6

i1

]yn

O(h

2

0=1，1

9,

,

,2 11

h(9

19

5

n2nn

h5y(5)

O(h6y(xn1)

y(xn)

x

f(x,

y(x))dx

x

F(x)dx

xnr或xn+1,xn (x次插值多項(xiàng)式r(x)代替F(x)求積分,即得r

3，過節(jié)點(diǎn)xnxn1xn2xn3F(x)的三次插值多項(xiàng)式3(x)

li(x)F(xnii0其中

(x)

(x

xn)(x

xn1)(x3

xn2)(x

xn3

(i

(x

xni)(xnij0ji

xnjFFh[55F對(duì)上式用yn

yn1代替y(xn

y(xn1用

(xk

yk)F(xk)

f(xk

y(xk

))

n,

2,

3則

h(55

59

37

9

n2

n3這就是四階Adams顯式公式。由于積分區(qū)間在插[xn3,xn]外面，又稱為四階Adams外插公式由插值余項(xiàng)公式可得其局部截?cái)嗾`差

F(4)()

y(5)()Rn1

(x

xnj)dx

(x

xnj

j

jRn1

y(5)()

(x3j3

xnj)dx

h5y(5)

同樣，如果過節(jié)點(diǎn)xn＋1,xn,xn1,xn2的F(x)三次插值2項(xiàng)式

(x)

li(x)F(xnii1其中

(x)

(x

xn1)(x

xn)(x3

xn1)(x

xn2

(x

xni)(xnijj

xnj代替F(x)求積分，即得四階Adams隱式公

h(9

19

5

n2其局部截?cái)嗾`差

h5y(5)(

由于積分區(qū)間在插[xn2,xn1內(nèi)，故Adams隱式稱為MilieMilie公式是（三

(Ham

g)公取r=2，并令

2

0,可得到Ham

g公

(9y )3h(

2f

n2

其局部截?cái)嗾`差

h5y(5)

O(h6nHamming公式是四階三步隱式公式n隱式法與顯式法的比一般地，同階的隱式顯式法精確，而且數(shù)值穩(wěn)定性也好。但在隱式公式中，通常很難解出yn+1,需而是將它們聯(lián)合使用：先用顯式公式求出y(xn+1)的預(yù)測(cè)值，記作yn1,再用隱式公式對(duì) —校正系

h(55

59

37

9

h

n2

y

[9f(x ,

)19

5

n2

4h(2

2 n3

n2 (9y ) h[f(x ,

)2

n2

說明 有時(shí)為提高精度，校正公式可迭代進(jìn)行多次，但代次數(shù)一般不超過3次用局部截?cái)嗾`差進(jìn)一步修 －校正公nny(xn1)

yn1

nn

h5y(5)

O(h6y(xn1)

nn

h5y(5)

O(h6兩式相減

yn1

yn1

h5y(5)

O(h6n h5y(5)n

720(

n1

yn1 y(xn1)

yn1

(yn1

yn1y(xn1)

(

n1

yn1由上面就得到多環(huán)節(jié)的

h(55

59

37

9

n2

n3

(cp 改

h y

[9f(x ,

)19

5

校

9

n2yn1

0(cn1

pn1 改cc

1(9ynn1n

n1輸

a,

f

算N,

ba

,

a,n計(jì)

fn1

f(xn1

yn1K

f

h,

K1

Kf(x h,

K2

f(x

h,

K

y

h

K

xa 輸出(xn

yn

若

3置n1

n返回否則，置n1 0,轉(zhuǎn)6計(jì)f3

f(x3

x

p

4(2f

2f

m

112

p

c1(9yy)3h[f(x,m)2ff y若n

N，置n1

xj1

xj

yj1

yj

fj1f(j

0,1,2),x

x3,y

y3,p

p0,

c0轉(zhuǎn)6；否則停機(jī)§5.常微分方程組與高階方程的數(shù)值一階方程組y'f(x,y,y , y(x)y0

(i1, y

采用向量的記y(y,y )T y0

(y0,y0 y0)T

y'f(x,

y(x0)

y0;f(f,f )T y'f(x,

y(x0)

y0;求解這一初值問題的四 —庫塔公式

h(k2k2kk h其h

k1

f

,yn k2

f

h,

2k1kf(xh,

hk

f(xh,yhk

hyi,n1h

22Ki3

Ki4 其

NKi1

fi(xn

y2n

yNn f(x

h,

hK,

h

,

h i

2n