第八章常微分方程數(shù)值解法

第8章

常微分方程數(shù)值解法§1.1§1

引言為什么要研究數(shù)值解法(8.1)

y(a)=a一階常微分方程初值問題的一般形式為

y￠=?(x,y)

,a￡x￡b其中?(x,y)是已知函數(shù),a為給定的初值.如果函數(shù)?(x,y)在區(qū)域{a￡x￡b,-￥

<y<￥}上連續(xù)且關(guān)于y滿足Lipschitz條件f

(x,

y)

-

f

(x,

y)

￡

L

y

-

y

,"

x,

y其中L>0為L(zhǎng)ipschitz常數(shù),則初值問題(8.1)有唯一解.所謂數(shù)值解法,就是設(shè)法將常微分方程離散化,建立差分方程,給出解在一些離散點(diǎn)上的近似值.§1.2

構(gòu)造數(shù)值解法的基本思想假設(shè)初值問題(8.1)的解y=y(x)唯一存在且足夠光滑.對(duì)求解區(qū)域[a,b]做剖分a=x0<x1<x2<…<xn<…<xN=b其中剖分節(jié)點(diǎn)xn=a+nh,n=0,1,…,N,h稱為剖分步長(zhǎng).數(shù)值解法就是求精確解y(x)在剖分節(jié)點(diǎn)xn上的近似值yn?y(xn),n=1,2,…,N.我們采用數(shù)值積分方法來建立差分公式.

在區(qū)間[xn,xn+1]上對(duì)方程(8.1)做積分,則有對(duì)右邊的積分應(yīng)用左矩形公式，則有(8.2)n

+1xxnf

(x,

y(x))dxy(xn+1

)

-

y(xn

)

=梯形公式xyo

ab左矩形公式y(tǒng)=?(x)bab

-

a+

f

(b)][

f

(a)2f

(x)dx

?baf

(x)dx

?

(b

-

a)

f

(a)右矩形公式baf

(x)dx

?

(b

-

a)

f

(b)中矩形公式baa

+

b)2f

(x)dx

?

(b

-

a)

f

((8.2)n

+1xxnn+1

nf

(x,

y(x))dxy(x

)

-

y(x

)

=稱為梯形公式.對(duì)右邊的積分應(yīng)用左矩形公式，則有y(xn+1

)

?

y(xn

)

+

hf

(xn

,

y(xn

))因此，建立節(jié)點(diǎn)處近似值yn滿足的差分公式+

hf

(xn

,

yn)

yn+1

=

yn0

y=a

,

n

=

0,1,2,

N

-1稱之為Euler公式.若對(duì)(8.2)式右邊的積分應(yīng)用梯形求積公式，則可導(dǎo)出差分公式

y0=a

,

n

=

0,1,2,

N

-12+

h

[

f

(x

,

y

)

+

f

(x

,

y

)]n

n

n+1

n+1

yn+1

=

ynn

+1若在區(qū)間[xn-1,xn+1]上對(duì)方程(8.1)做積分,則有xxn

-1n+1

n-1f

(x,

y(x))dxy(x

)

-

y(x

)

=對(duì)右邊的積分應(yīng)用中矩形求積公式，則得差分公式y(tǒng)n+1

=

yn-1

+

2hf

(xn

,

yn

)

0

y=a

,

n

=

0,1,2,

N

-1稱為Euler中點(diǎn)公式或稱雙步Euler公式.例1

利用Euler方法求初值問題12-

2

y

2

,0

￡

x

￡

21

+

xy￠=

y(0)

=

0的數(shù)值解.此問題的精確解是y(x)=x/(1+x2).解

此時(shí)的Euler公式為分別取步長(zhǎng)h=0.2

,0.1

,0.05,計(jì)算結(jié)果如下1nnn-

2

y

2

)1

+

x

2n+1y

=

y

+

h(=

0 ,

n

=

0,1,2

y0hxnyny(xn)y(xn)-ynh=0.20.000.000000.000000.000000.400.376310.34483-0.031480.800.542280.48780-0.054481.200.527090.49180-0.035291.600.466320.44944-0.016892.000.406820.40000-0.00682h=0.10.000.000000.000000.000000.400.360850.34483-0.016030.800.513710.48780-0.025901.200.509610.49180-0.017811.600.458720.44944-0.009282.000.404190.40000-0.00419h=0.050.000.000000.000000.000000.400.352870.34483-0.008040.800.500490.48780-0.012681.200.500730.49180-0.008921.600.454250.44944-0.004812.000.402270.40000-0.00227yn-1在Euler公式和梯形公式中,為求得yn+1,只需用到前一步的值yn,這種差分方法稱為單步法,這是一種自開始方法.Euler中點(diǎn)公式則不然,

計(jì)算yn+1時(shí)需用到前兩步的值yn

,,稱其為兩步方法,兩步以上的方法統(tǒng)稱為多步法.

在Euler公式和Euler中點(diǎn)公式中,需要計(jì)算的yn+1已被顯式表示出來,稱這類差分公式為顯式公式,而梯形公式中,需要計(jì)算的yn+1隱含在等式兩側(cè),稱其為隱式公式.隱式公式中,每次計(jì)算yn+1都需解方程,要比顯式公式需要更多的計(jì)算量,但其計(jì)算穩(wěn)定性較好.計(jì)算數(shù)值解的精度要比Euler公式好,但它屬于隱式公式,不便于計(jì)算.實(shí)際上,常將Euler公式與梯形公式結(jié)合使用:§2

改進(jìn)的Euler方法和Taylor展開方法§2.1

改進(jìn)的Euler方法從數(shù)值積分的角度來看,梯形公式2+

h

[

f

(x

,

y

)

+

f

(x

,

y

)]n

n

n+1

n+1

yn+1

=

yn

y0=a

,

n

=

0,1,2,

N

-1=

yn

+

hf

(xn

,

yn)n+1)]h2y[

k

]n+1

n+1[

k

+1]n+1[

f

(xn

,

yn

)

+

f

(x

,

y=

yn

+y[0]

y0=a

,

n

=

0,1,2,

N

-1由迭代法收斂的角度看，當(dāng)(e是給定的精度要求)時(shí),取就可以|<

en+1

n+1[

k

+1]

-

y[

k

]|

y.[

k

+1]n+1n+1y

=

y￡

L

<1,2

?yh

?f而且,只要+

hf

(xn

,

yn

)+

[

f

(xn

,

yn

)

+

f

(xn+1

,

yn+1

)]2hyn+1

=

ynyn+1

=

yn

y0=a

,

n

=

0,1,2,

N

-1保證迭代公式收斂,

而當(dāng)h很小時(shí),收斂是很快的.通常采用只迭代一次的算法:稱之為改進(jìn)的Euler方法.

這是一種單步顯式方法.改進(jìn)的Euler方法也可以寫成2+

h

(K

+

K

)1

2K1

=

f

(xn

,

yn

)yn+1

=

yn=a

,

n

=

0,1,2,

N

-1K

2

=

f

(xn

+

h,

yn

+

hK1

)

y0的數(shù)值解,取步長(zhǎng)h=0.1.[精確解為y(x)=(1+2x)1/2.]y(0)=1例2

求初值問題y￠=y-2x/y ,

0￡x￡1解

(1)

利用Euler方法yn+1

=1.1yn

-

0.2xn

/

yn

0=1

,

n

=

0,1,2,9

y+

0.05(K1

+

K

2

)nn1

n=

y

-

2x

/

y

K

yn+1

=

yn

y0計(jì)算結(jié)果如下:2

Kn

1=1 ,

n

=

0,1,2,9y

+

0.1K2(x

+

0.1)-

n=

yn

+

0.1K1(2)利用改進(jìn)Euler方法nxnEuler方法yn改進(jìn)Euler法yn精確解y(xn)0011110.11.11.0959091.09544520.21.1918181.1840961.18321630.31.2774381.2662011.26499140.41.3582131.3433601.34164150.51.4351331.4164021.41421460.61.5089661.4859561.48324070.71.5803381.5525151.54919380.81.6497831.6164761.61245290.91.7177791.6781681.6733201011.7847701.7378691.732051§2.2

差分公式的誤差分析在節(jié)點(diǎn)xn+1的誤差y(xn+1)-yn+1,不僅與yn+1這一步計(jì)算有關(guān),而且與前n步計(jì)算值yn,yn-1,…,y1都有關(guān).為了簡(jiǎn)化誤差的分析,著重研究進(jìn)行一步計(jì)算時(shí)產(chǎn)生的誤差.即假設(shè)yn=y(xn),求誤差y(xn+1)-yn+1,這時(shí)的誤差稱為局部截?cái)嗾`差,它可以反映出差分公式的精度.如果單步差分公式的局部截?cái)嗾`差為O(hp+1),則稱該公式為p階方法.這里p為非負(fù)整數(shù).顯然,階數(shù)越高,方法的精度越高.研究差分公式階的重要手段是Taylor展開式,一元函數(shù)和二元函數(shù)的Taylor展開式為:23!2!y

(x

)y(x

)

=

y(xh

+

y

(xn

)

h3

+n+

h)

=

y(xn

)

+

y￠(xn

)h

+nn+1++2

2k?2

f

(x

,

y

)?2

f

(x

,

y

)h

+

21

?2

f

(x

,

y

)k?yhk

+h

+?x?f

(x

,

y

)

n

n

?y2

n

n

?x?y

n

n

2!

?x2?f

(xn

,

yn

)n

nf

(x

+

h,

y

+

k)

=

f

(xn

,

yn

)

+nn另外,在yn=y(xn)的條件下,考慮到y(tǒng)￠(x)=?(x,y(x)),則有y￠(xn)=?(xn,y(xn))=?(xn,yn)=?nny

￠(x

)=nnnn

nf

?f+

fdx

?x

?yd

?[

f

(x

,

y(x

))]

=2

2?x

?y

?y?y2?x2

?x?y?2

f

?2

f

?2

f

?f

?f

?f=

n

+

2

n

fn

+

n

fn

+

n

n

+(

n

)

fnnnnnf

?f[

+

f

]dx

?x

?yd

?y￠(x )

=對(duì)Euler方法，有

yn+1=yn+h?(xn,yn)22!h

+y

(x

)y(x

)

=

y(x

+

h)

=

y(x

)

+

y￠(x

)h

+nn+1

n

n

n從而有:=yn+?(xn,yn)h+O(h2)y(xn+1)-yn+1=O(h2)所以Euler方法是一階方法.再看改進(jìn)Euler方法,因?yàn)?hKn?yh

+

?f

n?x+

?f

nK

2

=

f

(xn

+

h,

yn

+

hK1

)

=

f312122h

Kh

K

n

+

O(h

)2

n

?y2?2

f+

n

?x?yh

+

22

?x21

?2

f

?2

f+可得而n

nnf

?f

n

+

?n

y

=

y

+

f h

+2

?x

?yh

2n+14fn

n

+

O(h

)n

n

?y2?2

ff

+

n

?x?y+

24

?x2h3

?2

f

?2

f+3!22y

(x

)n+

O(h

4

)h

+

y

(xn

)

h3y(x

)

=

y(xn

)

+

y￠(xn

)h

+n+1n

f

?f

n

+

?f

n

=

yn

+

f

n

h

+2

?x

?yh

2422nn

n

+

O(h

)n

+(

)

f

n

n

n

?x

?y

?y?f

?f

?ff

+

n

?y2?2

ff

+

n

?x?y?2

f+

26

?x2h3

?2

f+從而有: y(xn+1)-yn+1=O(h3)所以,改進(jìn)的Euler方法是二階方法.§2.3

Taylor展開方法設(shè)y(x)是初值問題(8.1)的精確解,利用Taylor展開式可得………………n+1

n

nh

p

+1(

p

+1)!2!

P!y

(

p

)

(x

)y

(x

)y

(

p

+1)

(x)y(x

)

=

y(x

)

+

y￠(x

)h+

n

h

2

+

+

n

h

p

+p

+1(

p

-1)n

n(xn

,

y(xn

))

+

O(h

)2!

P!h

2

h

p(1)=

y(xn

)

+

hf

(xn

,

y(xn

))

+

f

(x

,

y(x

))

+

+

fn

n

n

n2!

P!因此，可建立節(jié)點(diǎn)處近似值yn滿足的差分公式h

2

h

pn+1

n

n

ny

=

y

+

hf(x

,

y

)

+

f

(1)

(x

,

y

)

+

+

f

(

p

-1)

(x

,

y

)0

y=a

,

n

=

0,1,2,

N

-1ff

2)

2?y+

(?x

?y?f

?f

?f+?y

2?2

f+

2

f

+?x

2

?x?y?2

f

?2

ff

(

2)

(x,

y)

=稱之為p階Taylor展開方法.其中f

(1)

(x,

y)

=

?f

(x,

y)

+

?f

(x,

y)

f

(x,

y)?x

?y可見,公式的局部截?cái)嗾`差為:y(xn+1)-yn+1=O(hp+1).所以,此差分公式是p階方法.由于Taylor展開方法涉及很多復(fù)合函數(shù)?(x,y(x))的導(dǎo)數(shù)的計(jì)算,比較繁瑣,因而很少直接使用,經(jīng)常用它為多步方法提供初始值.然而,Taylor展開方法給出了一種構(gòu)造單步顯式高階方法的途徑.§3

Runge-Kutta方法§3.1 Runge-Kutta方法的構(gòu)造

Euler方法可寫為+

hK

yn+1

=

ynn

nK

=

f

(x

,

y

)構(gòu)造差分公式……改進(jìn)的Euler方法可寫為2+

h

(K

+

K

)1

2K1

=

f

(xn

,

yn)

yn+1

=

ynK

2=

f

(xn

+

h,

yn

+

hK1

)+

h(l1

K1

+

l2K

2

+

+

lp

K

p

)21

122

nn……=

f

(x

+

a

h,

y

+

hb

K

)K

K1

=

f

(xn

,

yn

)

yn+1

=

ynp

-1i

=1KP=

f

(xn

+

a

P

h,

yn

+

h

b

pi

Ki

)其中{li,ai,bij}為待定參數(shù).若此公式的局部截?cái)嗾`差為O(hp+1),稱公式為p階Runge-kutta方法,簡(jiǎn)稱p階R-K方法.對(duì)于p=2的情形,應(yīng)有+

h(l1

K1

+

l2

K

2

)1

K

yn+1

=

yn(8.3)K

2=

f

(xn

,

yn

)=

f

(xn

+

ah,

yn

+

bhK1

)由于yn+1=yn+hl1?n+hl2(?n+ah?xn+bh?n?yn)+O(h3)=yn+h(l1+l2)?n+h2l2(a?xn+b?n

?yn)+O(h3)32hxn+

O(h

)f

+

f

n

f

yn

)y(xn+1

)

=

yn

+

hf

n

+

2

(所以,只要令l1+l2=1,

al2=1/2,

bl2=1/2(8.4)稱之為中點(diǎn)公式,或可寫為若取a=1,則得l1=l2=1/2,b=1,此時(shí)公式(8.3)就是改進(jìn)的

Euler公式;若取l1=0,則得l2=1,a=b=1/2,公式(8.3)為+

hK

21=

f

(x

,

y

)

K

yn+1

=

yn12121+

hK

)K

2nn

n=

f

(xn

+

h,

yn+

1

h,

y

+

1

hf

(x

,

y

))2

n

2

n

n+

hf

(xyn+1

=

yn一般地,參數(shù)由(8.4)確定的一族差分公式(8.3)統(tǒng)稱為二階R-K方法.高階R-K公式可類似推導(dǎo).下面列出常用的三階、四階R-K公式.6+

h

(K

+

4K

+

K

)1

2

31=

f

(xn

,

yn

)

K三階R-K公式

yn+1

=

ynK3

=

f

(xn

+

h,

yn

-

hK1

+

2hK

2

)四階標(biāo)準(zhǔn)R-K公式6+

h

(K

+

2K

+

2K

+

K

)1

2

3

41=

f

(xn

,

yn

)2

K=

f

(xn

+

1

h,

y

+

1

hK

)2

n

2

1K

yn+1

=

ynK

4=

f

(xn

+

h,

yn

+

hk3

)2

1121n+

hK

)

K

2

=

f

(xn

+

h,

y3

K=

f

(xn

+

1

h,

y

+

1

hK

)2

n

2

2解

四階標(biāo)準(zhǔn)R-K公式為21+

2K

+

2K

3

+

K

4

)y

=

y

+

h(Kn+1

n

6

1-

2xn

/

yny(0)=1例3

用四階標(biāo)準(zhǔn)R-K方法求初值問題y￠=y-2x/y ,

0￡x￡1的數(shù)值解,取步長(zhǎng)h=0.2.212212hK

)-

(2xn

+

h)

/(

yn

+112112hK

)n-

(2x

+

h)

/(

yn

+=

yn

+

hKK

3

=

yn

+

hKK1

=

ynK

2334=

y

+

hK

-

2(x

+

h)

/(

y

+

hK

)Knnn計(jì)算結(jié)果如下:nxnyny(xn)nxnyny(xn)00.01.001.0030.61.48331.483210.21.18321.183240.81.61251.612520.41.34171.341651.01.73211.7321也可以構(gòu)造隱式R-K方法,其一般形式為pr

=1+

h

lr

Kr

yn+1

=

ynKpr

n=

f

(xi

=1,

r

=1,2,,

plri

Ki

)+

a

r

h,

yn

+

h稱之為p級(jí)隱式R-K方法,同顯式R-K方法一樣確定參數(shù).如yn

=

y

+

1

h(K

+

K

)+1

n

2

1

2K

1221121+

hK

)K

2=

f

(xn

,

yn

)=

f

(xn

+

h,

yn

+

hK是二級(jí)二階隱式R-K方法,也就是梯形公式.但是p級(jí)隱式R-K方法的階可以大于p,例如,一級(jí)隱式中點(diǎn)公式為+

hK1yn+1

=

yn2

1+

1

h,

y

+

1

hK

)2

1=

f

(xKnn或?qū)憺閚+1+1

n

n+

1

h,

1

(

y

+

y

))2

2

n+

hf

(xyn

=

y它是二階方法.§3.2

變步長(zhǎng)Runge-Kutta方法以p階R-K方法為例討論.設(shè)從xn以步長(zhǎng)h計(jì)算y(xn+1)的n+1近似值為

y

(

h

)

,局部截?cái)嗾`差為n+1

n+1)

-

y

(

h)

=

Ch

p

+1y(x其中,C是與h無關(guān)的常數(shù).n+1y(x

)的近似值記為,其局部截?cái)嗾`差為于是有從而,得到事后誤差估計(jì)(

)如果將步長(zhǎng)減半,取h/2為步長(zhǎng),從xn經(jīng)兩步計(jì)算得到hy2n+1(

)2

p2n+1n+1p

+1

=

1

Ch

p

+1h?

2C(

)2y(x

)

-

yh可見,當(dāng)h2

p1y(x(

)2?)

-

y

(

h)y(x

)

-

yn+1n+1

n+1

n+11(

)(

)

2

p(

yhhn+12n+12n+1n+1-

y

(

h

)

)-1?y(x

)

-

y-

|￡(

)2n+1

n+1y

(

h)|

yh(

)he

成立時(shí),可取y(x2n+1n+1)

?

y.否則,應(yīng)將步長(zhǎng)再次減半進(jìn)行計(jì)算.(8.5)對(duì)于Euler方法,有§4

單步方法的收斂性和穩(wěn)定性的單步顯式方法可一統(tǒng)一寫為如下形式y(tǒng)n+1=yn+hF

(xn,yn,h)

y(a)=a§4.1

單步方法的收斂性求解初值問題

y￠=?(x,y)

,a￡x￡b其中F

(x,y,h)稱為增量函數(shù).F

(x,y,h)=?(x,y)對(duì)于改進(jìn)的Euler方法,有F

(x,y,h)=1/2[?(x,y)+?(x+h,y+h?(x,y))]設(shè)y(x)是初值問題(8.1)的解,yn是單步法y(xn),則稱單步法(8.5)是收斂的.可見,若方法(8.5)是收斂的,則當(dāng)hfi0時(shí),整體截?cái)嗾`差en=y(xn)-yn將趨于零.定理8.1

設(shè)單步方法(8.5)是p?1階方法,

增量函數(shù)F

(x,y,h)在區(qū)域{a￡x￡b,-￥

<y<+￥,0￡h￡h0}上連續(xù),且關(guān)于y滿足Lipschitz條件,初始近似y0=y(a)=a,則方法(8.5)是收斂的,且存在與h無關(guān)的常數(shù)C,使|y(xn)-yn|￡Chp證明

因?yàn)閱尾椒椒?8.5)是p階方法,則y(x)滿足定義8.1=hfi

0(8.5)產(chǎn)生的近似解.如果對(duì)任意固定的點(diǎn)xn,均有l(wèi)im

yn注意到y(tǒng)(xn+1)=y(xn)+hF

(xn,y(xn),h)+Rn(h)其中,局部截?cái)嗾`差|Rn(h)|￡C1hp+1,

記en=y(xn)-yn

,則有en+1=en+h[F

(xn,y(xn),h)-F

(xn,yn,h)]+Rn(h)利用Lipschitz條件得|en+1|￡(1+hL)|en|+C1hP+1遞推得到n-1ip

+1(1

+

hL)e

￡

(1

+

hL)

n

e

+

C

hn

0

10hLi

=0C

h

p

+1￡

(1

+

hL)

n

e

+

1

[(1

+

hL)

n -1]1+hL￡ehL,

(1+hL)n￡enhL￡eL(b-a)則有|en|￡|e0|eL(b-a)+C1hp/L(eL(b-a)-1)由于e0=y(a)-y0=0,所以有|en|￡C1hp/L(eL(b-a)-1)=Chp設(shè)?(x,y)連續(xù)且關(guān)于y滿足Lipschitz條件,對(duì)于Euler方法,由于F

(x,y,h)=?(x,y),故Euler方法是收斂的.對(duì)于改進(jìn)的Euler方法,由?(x,y)的Lipschitz條件有|F

(x,y,h)-F

(x,y*,h)|￡1/2|?(x,y)-?(x,y*)|+1/2|?(x+h,y+h?(x,y))-?(x+h,y*+h?(x,y*))|￡1/2L(2+hL)|y-y*|則當(dāng)h￡h0時(shí),F

關(guān)于y滿足常數(shù)為1/2L(1+h0L)的Lipschitz條件,因此改進(jìn)的Euler方法是收斂的.可類似驗(yàn)證各階R-K方法是收斂的.§4.2定義8.2單步方法的穩(wěn)定性對(duì)于初值問題(8.1),取定步長(zhǎng)h,用某個(gè)差分方法進(jìn)行計(jì)算時(shí),假設(shè)只在一個(gè)節(jié)點(diǎn)值yn上產(chǎn)生計(jì)算誤差d,即計(jì)算值‘yn=yn+d,如果這個(gè)誤差引起以后各節(jié)點(diǎn)值ym(m>n)的變化均不超過d,則稱此差分方法是絕對(duì)穩(wěn)定的.討論數(shù)值方法的穩(wěn)定性,通常僅限于典型的試驗(yàn)方程y￠=ly其中l(wèi)是復(fù)數(shù)且Re(l)<0.在復(fù)平面上,當(dāng)方法穩(wěn)定時(shí)要求變量lh的取值范圍稱為方法的絕對(duì)穩(wěn)定域,它與實(shí)軸的交集稱為絕對(duì)穩(wěn)定區(qū)間.將Euler方法應(yīng)用于方程y￠=ly,得到

yn+1=(1+lh)yn設(shè)在計(jì)算yn時(shí)產(chǎn)生誤差dn,計(jì)算值‘yn=yn+dn,則dn將對(duì)以后各節(jié)點(diǎn)值計(jì)算產(chǎn)生影響.記‘ym=ym+dm,m?n,由上式可知誤差dm滿足方程dm=(1+lh)dm-1=…=(1+lh)m-ndn

,

m?n可見,若要|dm|<|dn|,必須且只須|1+lh|<1,因此Euler法的絕對(duì)穩(wěn)定域?yàn)閨1+lh|<1,絕對(duì)穩(wěn)定區(qū)間是-2<Re(l)h<0.對(duì)隱式單步方法也可類似討論.如將梯形公式用于方程y￠=ly,則有yn+1=yn+h/2

l(yn+yn+1)解出yn+1得y1

-

1

lh1

+

1

lh=

2

yn+1

n2類似前面分析,可知絕對(duì)穩(wěn)定區(qū)域?yàn)?

2

<11

-

1

lh1

+

1

lh由于Re(l)<0,所以此不等式對(duì)任意步長(zhǎng)h恒成立,這是隱式公式的優(yōu)點(diǎn).一些常用方法的絕對(duì)穩(wěn)定區(qū)間為方

法方法的階數(shù)穩(wěn)定區(qū)間Euler方法1(-2

,

0)梯形方法2(-￥

,

0)改進(jìn)Euler方法2(-2

,

0)二階R-K方法2(-2

,

0)三階R-K方法3(-2.51

,

0)四階R-K方法4(-2.78

,

0)例4

考慮初值問題y(0)=1y￠=-30y ,

0￡x￡1取步長(zhǎng)h=0.1

,利用Euler方法計(jì)算y10?y(1). [y(x)=e-30x

]解

因y0=1,計(jì)算得y10=1024,而y(1)=9.357623·10-14.這是因?yàn)閘h=-3不屬于Euler方法的絕對(duì)穩(wěn)定區(qū)間.若取h=0.01,計(jì)算得y100=3.234477·10-16.若取h=0.001,計(jì)算得y1000=5.911998·10-14.若取h=0.0001,計(jì)算得y10000=8.945057·10-14.若取h=0.00001,計(jì)算得y100000=9.3156·10-14.單步顯式方法的穩(wěn)定性與步長(zhǎng)密切相關(guān),在一種步長(zhǎng)下是穩(wěn)定的差分公式,取大一點(diǎn)步長(zhǎng)就可能是不穩(wěn)定的.收斂性是反映差分公式本身的截?cái)嗾`差對(duì)數(shù)值解的影響;穩(wěn)定性是反映計(jì)算過程中舍入誤差對(duì)數(shù)值解的影響.只有即收斂又穩(wěn)定的差分公式才有實(shí)用價(jià)值.§5

線性多步方法由于在計(jì)算yn+1時(shí),已經(jīng)知道yn

,yn-1

,…,及?(xn,yn),?(xn-1,yn-1),…,利用這些值構(gòu)造出精度高、計(jì)算量小的差分公式就是線性多步法.§5.1

利用待定參數(shù)法構(gòu)造線性多步方法

r+1步線性多步方法的一般形式為當(dāng)b-1?0時(shí),公式為隱式公式,反之為顯式公式.參數(shù){ai,bi}的選擇原則是使方法的局部截?cái)嗾`差為y(xn+1)-yn+1=O(h)r+2這里,局部截?cái)嗾`差是指,在yn-i=y(xn-i),i=0,1,…,r的前提下,誤差y(xn+1)-yn+1.選取參數(shù)a,b0,b1,b2,使三步方法

yn+1=ayn+h(b0?n+b1?n-1+b2?n-2)為三階方法.rri

=0

i

=-1+

h

bi

f

n-iyn+1

=

ai

yn-i例5解

設(shè)yn=y(xn),yn-1=y(xn-1),yn-2=y(xn-2),則有?n=?(xn,y(xn))=y￠(xn)?n-1=?(xn-1,y(xn-1))=y￠(xn-1)=y￠(xn-h)=y￠(xn)-hy

￠(xn)+1/2h2y

￠(xn)-1/6h3y(4)(xn)+O(h4)?n-2=y￠(xn)-2hy

￠(xn)+2h2y

￠(xn)-4/3h3y(4)(xn)+O(h4)于是有yn+1=ay(xn)+h(b0+b1+b2)y￠(xn)-h2(b1+2b2)y

￠(xn)+h3(1/2b1+2b2)y

￠(xn)-h4/6(b1+8b2)y(4)(xn)+O(h5)y(xn+1)=y(xn)+hy￠(xn)+1/2h2y

￠(xn)+1/6h3y

￠(xn)+1/24h4y(4)(xn)+O(h5)若使: y(xn+1)-yn+1=O(h4)

,只要a,b0,b1,b2滿足:a=1,

b0+b1+b2=1,

b1+2b2=-1/2

,

b1+4b2=1/3解之得:12

3

121

20=

23

,

=

5b

=

-

4

,

ba

=1,

bn

+1于是有三步三階顯式差分公式y(tǒng)n+1=yn+h/12(23?n-16?n-1+5?n-2)§5.2

利用數(shù)值積分構(gòu)造線性多步方法因?yàn)閤xnn+1

nf

(x,

y(x))dxy(x

)

=

y(x

)

+n

+1設(shè)pr(x)是函數(shù)?(x,y(x))的某個(gè)r次插值多項(xiàng)式,則有xxnr

nnn+1p

(x)dx

+

Ry(x

)

=

y(x

)

+其中由此,可建立差分公式n

+1xxnrp

(x)dxyn+1

=

yn

+n

+1xxnrn(

f

(x,

y(x))

-

p

(x))dxR

=r選取不同的插值多項(xiàng)式pr(x),就可導(dǎo)出不同的差分公式.下面介紹常用的Adams公式.1.Adams顯式公式設(shè)已求得精確解y(x)在步長(zhǎng)為h的等距節(jié)點(diǎn)xn-r,…,xn上的近似值yn-r

,…,yn

,記?k=?(xk,yk),利用r+1個(gè)數(shù)據(jù)

(xn-r,?n-r),…,(xn,?n)構(gòu)造r次Lagrange插值多項(xiàng)式pr

(x)

=

ln-

j

(x)

f

n-

jj

=0其中由此,可建立差分公式由于rj

=

0,1,rl

(x)

=k

?

jk

=0

(xn-

j

-

xn-k

)(x

-

xn-k

)n-

jhbxxrnnn-

jn-

jy

=

y

+n

+1l

(x)dx

fj

=0n+1nxrxxnxn+

th)l

(x)dx

=n

+1n

+1k

?

jk

=0

(x

-

x

)n-

j

n-k(x