9常微分方程數(shù)值解法2021/7/171§9.0引言§9.1

歐拉（Euler）方法§9.2

改進(jìn)的歐拉（Euler）方法§9.3

龍格-庫塔（Runge-Kutta）法§9.4

線性多步法§9.5

相容性、收斂性與穩(wěn)定性9.0引言2021/7/1729.0.1基本概念9.0.2微分方程數(shù)值解法的一般過程9.0.3數(shù)值解法的三種類型9.0引言9.0.1基本概念常微分方程可分為初值問題和邊值問題。一階常微分方程初值問題的一般形式為dy

=

f

(x,

y)2021/7/173dx

y(x0

)

=

y0如果f(x，y)在[a，b]上連續(xù)，且滿足Lipschitz條件：$L>0(常數(shù))，"y1，y2，使得|

f(x，y1)

-

f(x，y2)|≤L

|y1

-y2|則微分方程的解存在唯一。注意：這只是個充分條件9.0引言2021/7/1749.0.1基本概念對于常微分方程，兩點(diǎn)邊值問題是要尋找方程在區(qū)間[a，b]邊界上滿足已知條件的解。區(qū)間兩端給出不同的邊界條件構(gòu)成不同類型的邊值問題這里只討論初值問題。9.0引言2021/7/1759.0.2

微分方程數(shù)值解法的一般過程對一些典型的微分方程(可分離變量方程，一階線性方程等等)，有可能找出它們的一般解表達(dá)式，然后用初始條件確定表達(dá)式中的任意常數(shù)，這樣解即能確定。這種解常稱為常微分方程的一般解或解析解。9.0引言2021/7/1769.0.2

微分方程數(shù)值解法的一般過程微分方程的數(shù)值解法：不求y=y(x)的精確表達(dá)式，而求區(qū)間[a，b]中離散點(diǎn)x0，x1，…，xN處的函數(shù)值。取[a，b]內(nèi)的一系列節(jié)點(diǎn)x0，x1，…，xN，用數(shù)值方法求得y(x)在每個節(jié)點(diǎn)xn

的值y(xn)的近似值yn，即yn

?y(xn)。而y0，y1，…，yN稱為微分方程的數(shù)值解。9.0引言2021/7/177xx0

x1

x2

…

xn

…yy0

y1

y2

…

yn

…對于常微分初值問題一般解

y

=

y(x)數(shù)值解符號規(guī)定y(xn)：一般解y=y(x)在x=xn處的精確值.yn：一般解y=y(x)在x=xn處的近似值.9.0引言2021/7/178求y(x)—>求y0，y1，…，yN方法：采用步進(jìn)式和遞推法將x離散化，x0<x1<…<xn

<xn+1…，一般采用等距步長xn=xn-1+

h，或

xn=x0+nh

n=1,2，…，Ny2fi

…fi

yn

fiyn-mfi

yn-m+1…fi

yn+1n-1n-m

y0

=

y(x0

)

y=

g(h,

xn

,

yn

,

y

n+1,,

y

)建立遞推公式通常利用Taylor公式或數(shù)值積分法。計(jì)算公式計(jì)算過程y0fi

y1fi9.0引言9.0.3微分方程數(shù)值解法的三種類型hhfi

0y

'(x)

=

lim

y(x

+

h)

-

y(x)hy

'(x)

?

y(x

+

h)

-

y(x)1.用差商代替導(dǎo)數(shù)由導(dǎo)數(shù)定義若h的值較小，則可用向前差商近似導(dǎo)數(shù)代入y’=

f

(x,y)可得h2021/7/179

y(x

+

h)

-

y(x)

?

f

(x,

y)即y(x+h)≈y(x)+hf(x,y)2021/7/17109.0引言0故原初值問題可離散化為

yn+1

=

yn

+

hf

(xn

,

yn

)

0

y

=

y(x

)

n

=

0,1,

2于是由初始值y(x0)=y0出發(fā)，可依次地計(jì)算出y1=y0+hf(x0,y0)y2=y1+hf(x0+h,y1)……yn=yn-1+hf(x0+(n-1)h,yn-1)……9.0引言2.使用數(shù)值積分的方法對于微分方程y’=f(x,y)兩邊求x0到x的定積分x0

dx

x0

x

dy

dx

=

x

f

(

x,

y(

x))dx即x0y(

x)

-

y(

x0

)

=

x

f

(t

,

y(t

))dt或?qū)憺?2021/7/1711xy(x)

=

y0

+

x

f

(t,

y(t))dt這就是與初值問題0

0

y(

x

)

=

y

y'=

f

(

x,

y)等價的積分方程。只要用某種數(shù)值積分方法便可建立起近似公式。9.0引言例：對積分部分應(yīng)用左矩形公式，則有xn+1

f

(t,

y(t))dtny(xn+1

)

-

y(xn

)

=

x?

(xn+1

-

xn

)

f

(xn

,

y(xn

))02021/7/1712故其對應(yīng)的數(shù)值公式為

yn+1

=

yn

+

hf

(xn

,

yn

)

0

y

=

y(x

)

n

=

0,1,

2,單步法：計(jì)算yn+1只需要知道yn即上一步計(jì)算值.9.0引言例：對積分部分應(yīng)用梯形公式，則有故其數(shù)值公式為2,

y

)]n+1

nn

n

n+1

n+1y

=

y

+

h

[

f

(x

,

y

)

+

f

(x[2021/7/1713n

nf

(x

,

y(x

))n+1n+1xn+1

-

xn2?+

f

(x

,

y(x

))]nxn+1

f

(t,

y(t))dty(xn+1

)

-

y(xn

)

=

x9.0引言公式22021/7/1714n

nnn+1n+1

n+1y

=

y

+

h

[

f

(x

,

y

)

+

f

(x

,

y

)]特點(diǎn)：①未知量yn+1不但出現(xiàn)在左端，而且出現(xiàn)在右端隱式公式：未知數(shù)yk+1同時出現(xiàn)在等式的兩邊。隱式公式一般需要通過迭代方法求解yk+1②計(jì)算yn+1需要知道yn，yn+1，即n步和n+1步的計(jì)算值.多步法：計(jì)算yn+1除需要知道yn外，還需要知道其他步的計(jì)算值.9.0引言3.使用泰勒公式在微分方程y′=f(x,y)中，y′是x及y(x)的函數(shù).由于精確值y(x+h)在h=0處的泰勒展式為(4)2021/7/171526

24h2

h3

h4y

"(x)

+

y

￠(x)

+

y(x)

+y(x

+

h)

=

y(x)

+

hy

'(x)

+取前兩項(xiàng)，有y(x+h)?y(x)+hy’(x)0根據(jù)y

′=f(x,y)，得公式

yn+1

=

yn

+

hf

(xn

,

yn

)

0

y

=

y(x

)

n

=

0,1,

2,9.1

Euler方法2021/7/17169.1.19.1.2Euler方法

Euler方法的誤差分析9.1

Euler方法dy

=

f

(x,

y)2021/7/1717dx

y(x0

)

=

y0求解一階微分方程初值問題9.1.1 Euler方法即用差分方程初值問題

yn+1

=

yn

+

hf

(xn

,

yn

)

0

0

y

=

y(x

)

n

=

0,1,

2,9.1

Euler方法9.1.1 Euler方法yn+1

=yn

+h

f

(xn,

yn)

(n=0,...,N-1)亦稱為歐拉折線法x0

x1其幾何解釋見右圖，即在區(qū)間[x0，x1]上取左端點(diǎn)的導(dǎo)數(shù)值計(jì)算右端點(diǎn)函數(shù)值。由于公式中采用向前差分近似導(dǎo)數(shù)，上式稱為向前歐拉公式。由于公式右端變量均為已知，可直接求出yn+1，又稱為顯式歐拉公式。2021/7/17189.1

Euler方法9.1.1 Euler方法如果采用向后差商近似導(dǎo)數(shù)1hy(x1

)

-

y(x0

)y￠(x

)

?x0x1y(

x1

)

?

y

+

h f

(

x

,

y(

x

))0

1

1yn+1

=

yn

+

h

f

(xn+1,

yn+1

)

(n

=

0,

...

,

N

-1)這是一個隱式公式，稱為隱式歐拉公式，又稱為向后歐拉公式。2021/7/179.1

Euler方法9.1.1 Euler方法例1.

y(0)

=

1y

y￠=

y

-

2x

0

￡

x

￡

1解:)nnnn2xyyn+1

=

yn

+hf

(xn

,

yn

)

=

y

+

h(

y

-n

=

1,

2,,10h

=

0.1用Euler

方法求解初值問題：y2021/7/17f

(x,

y)

=

y

-

2xx0

=

0,

N

=10,

b

=1,

y0

=1顯然由Euler

公式,有9.1

Euler方法9.1.1 Euler方法得)00001y2

xy

=

y

+

h(

y

-1=

1

+

0.1(1

-

2

·

0

)

=

1.1)2021/7/17112y2

xy

=

y1

+

h(

y1

-1.1=1.1+

0.1(1.1-

2

·0.1)

=1.1918依此類推。9.1

Euler方法9.1.1 Euler方法nyn+1n+1n+1)

=

y

+

h(

y-

2xn+1

)n

=1,

2,,1012021/7/17n+1y

=(

y

+

y2

-8h(1-

h)x

)n

n

n+12(1-

h)由向后Euler

公式,有yn+1

=

yn

+

hf

(xn+1

,

y解得nxny(xn)ynqynh10.11.095445115010331.100000000000001.0907375368352120.21.183215956619921.191818181818181.1740757612934830.31.264911064067351.277437833714721.2512485067965140.41.341640786499871.358212599560291.3230934977520950.51.414213562373101.435132918657801.3901780746271960.61.483239697419131.508966253566331.4528699233246470.71.549193338482971.580338237655221.5113768371646580.81.612451549659711.649783431047711.5657672354519890.91.673320053068151.717779347860091.615977254482521011.732050807568881.784770832497981.661807042621092021/7/170.10.20.30.40.50.60.70.80.91101.11.21.31.41.51.61.71.8向向Euler法準(zhǔn)準(zhǔn)準(zhǔn)向向Euler法圖9.22021/7/179.1

Euler方法2021/7/179.1.1 Euler方法Euler方法的幾何意義是用一條自點(diǎn)(x0,y0)出發(fā)的折線去逼近積分曲線y=f(x)。因此，這種方法又稱為折線法。顯然，歐拉法簡單地取折線的端點(diǎn)作為數(shù)值解，精度不高。提高精度的一種方法是取較小的h，例1中若將h減少為原來的1/10后的計(jì)算結(jié)果見圖9.3。當(dāng)然，h減少必然導(dǎo)致計(jì)算量增加，需要在精度和計(jì)算量之間取得平衡。圖9.30.20.30.40.50.60.70.80.9110

0.12021/7/171.11.21.31.41.51.61.71.89.1

Euler方法2021/7/179.1.2 Euler方法的誤差分析單步法計(jì)算公式一般可寫成yn+1

=

yn

+

hj

(xn，yn，xn+1，

yn+1，h)如果j

中不含yn+1

，則方法為顯式的，否則為隱式的一步計(jì)算產(chǎn)生的誤差：設(shè)第n步得到是精確解y(xn)，據(jù)此按單步公式計(jì)算第n+1步得到y(tǒng)n+1與精確解y(xn+1)之間的誤差。9.1

Euler方法2021/7/179.1.2 Euler方法的誤差分析定義9.1

設(shè)y(x)是方程9.1的精確解，則Rn+1=y(xn+1)-

y(xn)-

hj

(xn，

y(xn)，xn+1，

yn+1，h)稱為單步法的局部截?cái)嗾`差，簡稱截?cái)嗾`差。可由Taylor展開得到截?cái)嗾`差。9.1

Euler方法29.1.2 Euler方法的誤差分析

Euler法的截?cái)嗾`差Rn+1

=

y(xn+1

)

-

y(xn

)

-

hf

(xn

,

y(xn

))=

y(xn

+

h)

-

y(xn

)

-

hy￠(xn

)=

1

h2

y

￠(x)也可表示為22021/7/1712nh

y

(x

)n+1￠R

=+

O(h3

)

=

O(h2

)該公式只考慮了誤差的一個來源。9.1

Euler方法2021/7/179.1.2 Euler方法的誤差分析如果某種數(shù)值方法的局部截?cái)嗾`差為O(hP+1)，稱該方法是p階的。Euler方法是1階的。定義9.2在不計(jì)舍入誤差的條件下，方程(9-1)的精確解y(xn+1)與數(shù)值解yn+1之差稱為整體截?cái)嗾`差，記為en+1=

y(xn+1)-

yn+19.1

Euler方法22021/7/172n+1|

R

|=|

1

h

y

(2￠x)

|￡

1

h2

Mn

=1,

2,9.1.2 Euler方法的誤差分析

Euler方法的整體截?cái)嗾`差：函數(shù)f(x，y)關(guān)于y滿足Lipschitz條件|f(x，y1)-f(x，y2)|≤

L|y1-y2|L為Lipschitz常數(shù)設(shè)函數(shù)f(x，y)充分光滑，y(x)在[a，b]上二階連續(xù)可微，則局部截?cái)嗾`差有界，即存在M>0，對任意

x?[a，b]，都有9.1

Euler方法9.1.2 Euler方法的誤差分析記yn+1

=

y(xn

)

+

hf

(xn

,

y(xn

))h22021/7/17則有|

en+1

|=|

y(xn+1

)-yn+1

|￡|

y(xn+1

)-yn+1

|

+|

yn+1

-yn+1

|=|

Rn+1

|

+

|

y(xn

)

+

hf

(xn

,

y(xn

))

-

yn

-

hf

(xn

,

yn

)

|￡|

Rn+1

|

+

|

y(xn

)

-

yn

|

+h

|

f

(xn

,

y(xn

))

-

f

(xn

,

yn

)

|￡|

Rn+1

|

+

|

y(xn

)

-

yn

|

+hL

|

y(xn

)

-

yn

|=|

Rn+1

|

+(1+

hL)

|

y(xn

)

-

yn

|￡

M

+

(1+

hL)

|

en

|29.1

Euler方法令22021/7/179.1.2 Euler方法的誤差分析h2|

en+1

|￡

2

M

+

(1+

hL)

|

en

|h2T

=

M

,S

=1+

hL依次下去，有0n|

e

|n-1|

en+1

|￡

T

+

S

|

en

|￡

T

+

S

(T

+

S

|

en-1|)

=

T

+

ST

+

S

2

|

e

|￡

T

+

ST

+

S

2

(T

+

S

|

e

|)

=

T

+

ST

+

S

2T

+

S

3

|

e

|n-2

n-2S

k

+

Sn+1k

=0￡

￡

T9.1

Euler方法9.1.2 Euler方法的誤差分析但

e0=y(x0)-

y0=0因此nkS

n+1n+1n+1k

=0h2

M

(1+

hL)n+1

-1

hM=

T

==

[(1+

hL)

-1]S

-1

2 (1+

hL)

-1

2L|

e

|￡

T

S?根據(jù)公式2!xx2(x

?0)e

=1+

x

+ +

?

0有|

e2021/7/172Ln+1|￡

hM

[ehL(n+1)

-1]9.1

Euler方法9.1.2 Euler方法的誤差分析因?yàn)?/p>

(n+1)h≤Nh≤b-a得因此，整體誤差en+1與步長h同階當(dāng)hfi

0時，eNfi

0。hM2L|

en+1

|￡[eL(b-a

)

-1]

=

O(h)hM2021/7/172LhL(n+1)n+1|

e

|￡

[e-1]9.2

改進(jìn)的(歐拉)Euler方法2,

y

)]n+1

nn

n

n+1

n+1y

=

y

+

h

[

f

(x

,

y

)

+

f

(x9.2.1

梯形公式在9.0節(jié)里提到了梯形公式42021/7/1712n梯形公式的局部截?cái)嗾`差為h3n+1R

=

-

y

￠(x

)

+

O(h

)梯形公式是二階方法9.2

改進(jìn)的(歐拉)Euler方法9.2.1

梯形公式梯形公式需迭代求解，由簡單迭代法，迭代公式為n

n

nnh2,

y(k

)

)n+1(k

+1)n+1

n+1

y(0)=

y

+

hf

(x

,

y

)y

=

y

+[

f

(xn

,

yn

)

+

f

(x

n+1(k

=

0,1,)f(x,y)滿足Lipschitz條件，則兩次迭代值的差為h

hL2

22021/7/17|

y(k

+1)

(k

)(k

-1)(k

)

(k

-1)n+1

n+1n+1

n+1

n+1

n+1

n+1

n+1-

y

|=|

f

(x

,

y(k

)

)

-

f

(x

,

y)

|￡

|

y

-

y

|可以限定每步的迭代次數(shù)。由此導(dǎo)出改進(jìn)Euler方法9.2

改進(jìn)的(歐拉)Euler方法9.2.1

改進(jìn)Euler法思想：在梯形公式迭代時，每步只迭代一次。先用Euler公式求yn+1的一個初步近似值，稱為預(yù)測值；再用梯形公式進(jìn)行校正，得到近似值yn+1。22021/7/17n

nnn+1n+1

n+1

yn+1

=

yn

+

hf

(xn

,

yn

)

y=

y

+

h

[

f

(x

,

y

)

+

f

(x

,

y

)]預(yù)測校正9.2

改進(jìn)的(歐拉)Euler方法9.2.1

梯形公式可以改寫成qn12n

p

yp

=

yn

+

hf

(xn

,

yn

)

y=

y

+

hf

(x

+

h,

y

)

yn+1

=(

yp

+

yq

)22021/7/17n

nnn+1n+1

n+1

yn+1

=

yn

+

hf

(xn

,

yn

)

y=

y

+

h

[

f

(x

,

y

)

+

f

(x

,

y

)]9.2

改進(jìn)的(歐拉)Euler方法2021/7/179.2.1

梯形公式算法9.11)輸入a，b、f(x,y)，計(jì)算總步數(shù)N，初值y00輸出(x，y)2)置

h

b

-

a

n=0，x=a，y=

yN3)計(jì)算yp=y+h

f(x,y)yq=y+h

f(x+h,

yp)置(yp

+

yq)/2y，x+h

=

x

輸出（x，y）4)如果n<N-1，置n+1n，轉(zhuǎn)3)；否則，停機(jī)11.11.21.31.41.51.61.71.810

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9用預(yù)測-修正方法計(jì)算y(1)的近似值比將h減少為原來的1/10的向前和向后Euler方法的計(jì)算結(jié)果還要好點(diǎn)。2021/7/172021/7/17解：由題意改進(jìn)歐拉法的具體形式為要求步長h=0.2

，計(jì)算y(1.2)，

y(1.4)

的近似值

y(1)

=

19.2

改進(jìn)的(歐拉)Euler方法例：用歐拉預(yù)報-校正公式求解初值問題

y'+y

+

y

2

sin

x

=

0p

pcnn

nnp(y

p

+

yc

)2n1-y

-y

sin

x

y

=

y

+

h(-

y

-

y

sin

x

)

y

=

y

+

h(

)

yn+1

=n+122f

(x,

y)=

-y

-

y

2

sin

x9.2

改進(jìn)的(歐拉)Euler方法由y(1)=y0=1知y0=1，取步長h=0.2，計(jì)算如下n=0時(

))200

0

021212pc

0

p

p

1

y

=

y

+

h

-y

-

y

sin

x=1+

0.2

(-1-12

sin1)?

0.631706

y

=

y

+

h

(-y

-y

sin

x=1+

0.2

(-0.631706

-

0.6317062

sin1.2)?

0.799272

y

=1(yp

+

yc

)=(0.631706

+

0.799272)?

0.7154892021/7/172021/7/179.2

改進(jìn)的(歐拉)Euler方法n=1時(

))21

1

121212p

1c

1

p

p

2

y

=

y

+

h

-y

-

y

sin

x=

0.715489

+

0.2

(-0.715489

-

0.7154892

sin1.2)?

0.476964

y

=

y

+

h

(-y

-y

sin

x=

0.715489

+

0.2

(-0.476964

-

0.4769642

sin1.4)?

0.575259

y

=2(yp

+

yc

)=(0.476964

+

0.575259)?

0.526112y(1.2)

?

y1

?

0.715489?

0.526112y(1.4)

?

y29.3龍格-庫塔(Runge-Kutta)法9.3.1

Runge-Kutta法的基本思想設(shè)初值問題)(

p

)kk有解y(x),由Tayler公式得:h

2

h

p2!

p

!p

+1y

(

xk

+1

)

=

y

(

xk

)

+

h

y￠(

xk

)

+(

x

)

+

O

(hy

￠(

x

)

+

+

y令p

階Taylor方法.p=1時即為Euler公式.稱之為Taylor級數(shù)法.其中(

p

)k2021/7/1745h

2

h

pp

!y

k

+1=

y

k

+

h

y

k￠

+

2!

y

k￠+

+

ykky

(i

)?

y

(i

)(

x

),

i

=

0,1,

2,,

pk

k當(dāng)y

(i

)?y

(i

)(x

)

時,有y(

xk

+1

)

-

yk

+1=O(hp

+1

)

.此時①為y

=

f

(

x,

y),

y(

x0

)

=

y0

,

a

￡

x

￡

b9.3龍格-庫塔(Runge-Kutta)法例:取步長h=0.1,用一階、二階和四階Taylor方法求解下列初值問題12021/7/1746

y￠=

y

2, 0

￡

x

￡

.2

y

(0)

=

1解: (1)

一階Taylor法nn

+1y

=

y

+

hy￠=

y

+

0.1

y

2n

n

nn01234yn+11.11.2211.370081.557791.800469.3龍格-庫塔(Runge-Kutta)法22021/7/174730.122!2

ynyn

+1(2)二階Taylor法y

￠=

(

y

2

)￠=

2

yy

￠=

2

y

3=

yn

+

0.1

yn

+k01234yk+11.111.246891.421751.652631.970889.3龍格-庫塔(Runge-Kutta)法22021/7/17483

450.12

0.130.142!

3!4!24

ynn+1

n

n(3)

四階Taylor法y

￠=

(2

y3

)￠=

6

y2

y￠=

6

y4y(4)

=

24

y3

y￠=

24

y5

y

=

y

+

0.1y

+2

yn

+

6

yn

+k01234yk+11.11111.249961.428481.666441.999429.3龍格-庫塔(Runge-Kutta)法2021/7/1749yk+11.11111.249961.428481.666441.99942k01234y(xk+1)1.11111.250001.428571.666672.00000yn+11.11.2211.370081.557791.80046yk+11.111.246891.421751.652631.970889.3龍格-庫塔(Runge-Kutta)法2021/7/179.3.1

Runge-Kutta法的基本思想例中的函數(shù)比較簡單，如果函數(shù)表達(dá)式復(fù)雜，相應(yīng)的導(dǎo)數(shù)計(jì)算非常繁瑣，如y￠(x)

=

f

(x,

y),y

￠(x)

=

fx

(x,

y)

+

f

y

(x,

y)

f

(x,

y),如何快速估計(jì)y的各階導(dǎo)數(shù)影響是簡化計(jì)算公式，提高計(jì)算效率的關(guān)鍵。9.3龍格-庫塔(Runge-Kutta)法2021/7/1751這里稱為[xn

,xn+1]上的平均斜率。y(xn+1)

=

y(xn)

+

hf(x，y(x))

=

y(xn)

+

hK*K

*

=

f

(x,

y

(x))9.3.1Runge-Kutta法的基本思想根據(jù)Lagrange中值定理y(xn+1)

=

y(xn)

+hy

￠(x)

x?

(xn+1，xn)由方程可知y￠(x)

=

f(x，y(x))因此9.3龍格-庫塔(Runge-Kutta)法2021/7/1752即，取K*=K1，其精度為一階1

n

n

K

=

f

(x

,

y

)9.3.1Runge-Kutta法的基本思想從已知的兩個公式看顯式Euler公式y(tǒng)n+1

=

yn+hf(xn，yn)，改寫成

yn+1

=

yn

+

hK19.3龍格-庫塔(Runge-Kutta)法可寫成122

21

n

nK

Kn+1

n2

n+1

n

1

y=

y

+

h(+

)K

=

f

(x

,

y

)K=

f

(x

,

y

+

hK

)即取K*=K1/2+K2/2，其精度為二階9.3.1Runge-Kutta法的基本思想改進(jìn)的Euler公式22021/7/1753n

nnn+1n+1

n+1

yn+1

=

yn

+

hf

(xn

,

yn

)

y=

y

+

h

[

f

(x

,

y

)

+

f

(x

,

y

)]9.3龍格-庫塔(Runge-Kutta)法2021/7/17549.3.1Runge-Kutta法的基本思想共同點(diǎn)：利用f(x,y)的某些函數(shù)值的線性組合計(jì)算y(xn+1)的近似值yn+1

。Euler公式每步算一次f(x,y)的值K1作為K*的近似值，是y(xn+1)在xn處的一階Taylor多項(xiàng)式，為一階方法。改進(jìn)的Euler公式每步算兩次f(x,y)的值K1和K2，用它們的加權(quán)平均值作為K*近似值，是y(xn+1)在xn處的二階Taylor多項(xiàng)式，為二階方法。p2021/7/17559.3龍格-庫塔(Runge-Kutta)法9.3.2

RK方法的構(gòu)造在每一步[xn，xn+1]內(nèi)多取幾個點(diǎn)(如p個)的Ki值，用其加權(quán)平均值作為K*的近似值yn+1

=

yn

+

h

ci

Kii=1i-1j

=1

Ki

=

f

xn

+

aih,

yn

+

hbij

K

j

其中ci

、ai

、bij

為待定參數(shù),

且

a1=0

.稱為p級Runge-Kutta方法計(jì)算公式.9.3龍格-庫塔(Runge-Kutta)法2021/7/17569.3.2

RK方法的構(gòu)造其中K1

=

f(xn，yn)K2

=

f(xn+a2h，yn+b21hK1)K3

=

f(xn+a3h，yn+h(b31K1+

b32K2))………9.3龍格-庫塔(Runge-Kutta)法22021/7/17573231

12!

3!2!(

p

)nnr

hh2h

pp

!n

+1

n

1p

+1

)n

+1y

=

y

+

r

h

+

r

h

+y(

x

)

=

y(

xn

)

+

h

y￠(

xn

)

+y

(

x

)

+

O(h+y

￠(

x

)

+

+9.3.2RK方法的構(gòu)造式中待定參數(shù)的確定:先將公式右端在(xn

,yn

)處展成h的冪級數(shù)(即將yn+1

展成h的冪級數(shù));再將y(xn+1)作Taylor

級數(shù)展開;最后比較兩式中hn(n=0，1，2，…)的系數(shù),以確定出所有待定參數(shù).9.3龍格-庫塔(Runge-Kutta)法2021/7/17589.3.2RK方法的構(gòu)造如要求:r

=

y

￠(

x

),

r

=

y

￠(

x

),

,

r

=

y

(

p

)

(

x

)1

n

2

n

p

n即可得p個方程,從而確定出待定參數(shù).代入表達(dá)式即可得到計(jì)算公式.如果要求兩個表達(dá)式的前p+1項(xiàng)完全重合,即局部截?cái)嗾`差達(dá)到

O(h

p

+1

)

,

則稱上式為p

階r級

的Runge-Kutta方法.常用的是r=2,3,4級的R-K方法.且適當(dāng)選取參數(shù)使得p

=r9.3龍格-庫塔(Runge-Kutta)法2021/7/17591

K=

f

(

xn

,

yn

)

K=

f

(

xn

+

a2

h,

yn

+

b21hK1

)

2yn

+1

=

yn

+

hc1

f

(

xn

,

yn

)

+

hc2

f

(

xn

+

a2

h,

yn

+

B21hK1

)記f

=

f

(xn

,

yn

),

fx

=

fx

(xn

,

yn

),

f

y

=

f

y

(xn

,

yn

)9.3.2

RK方法的構(gòu)造Runge-Kutta方法的推導(dǎo)(以r

=2為例):當(dāng)r

=2

時

yn

+1

=

yn

+

h(c1

K1

+

c2

K

2

)9.3龍格-庫塔(Runge-Kutta)法2021/7/1760又=

y

+

(c

+

c

)

f

h

+

c

(a

f

+

b

K

f

)h2

+

O(h3

)n

1

2

2

2

x

21

1

yy￠=

f

(

xn

,yn

)

=

f

,

y

￠=

f

x

+

f

y

f9.3.2

RK方法的構(gòu)造則K1

=

f

(

xn

,yn

)K2

=

f

(

xn

+

a2

h,

yn

+

b21hK1

)=

f

(

x

,

y

)

+

h(a

f

+

b

K

f

)

+

O(h2

)

,n

n

2

x

21

1

yyn+1

=

yn

+

(c1

K1

+

c2

K2

)h9.3龍格-庫塔(Runge-Kutta)法32021/7/1761322nnx

yh

2n

+1

n

ny

(

x

)

=

y

(

x

)

+

h

y￠(

x

)

+

y

￠(

x)

+

O

(h

)=

y

+

h

f

+(

f

+

f f

)

+

O

(h

)為使局部截?cái)嗾`差為O

(h3

),比較上述兩式右端同次冪系數(shù),應(yīng)取9.3.2

RK方法的構(gòu)造h

29.3龍格-庫塔(Runge-Kutta)法2021/7/1762這是一個四個參數(shù)三個方程的非線性方程組.它有一個自由度，因此有無窮多個解.稱滿足上述方程組的一組公式為二級二階Runge-Kutta方法.2

2c

a

=

1/

2c

b

=

1/

2

2

219.3.2

RK方法的構(gòu)造c1

+

c2

=

12021/7/17639.3龍格-庫塔(Runge-Kutta)法(1)預(yù)估-校正公式c1

=

c211nn=

1/

2,

a2

=

b21

=

1=

yn

+

h

(

K

1

+

K

2

)

/

2

yn

+1

K=

f

(

xn

,

yn

)

K=

f

(

x

+

h

,

y

+

h

K

)

2(2)122c

=

0,

c=

1,

a2

=

b21=

112n=

yn

+

h

K

2

yn

+1

K=

f

(

xn

,

yn

)

K=

f

(

xn

+

h

,

y

+

h

K

)2

2

1中點(diǎn)公式9.3.2

RK方法的構(gòu)造常用的二階Runge-Kutta方法:9.3龍格-庫塔(Runge-Kutta)法2021/7/1764注:

二級Runge-Kutta方法的精度最高是二階的,不可能達(dá)到三階.要提高計(jì)算方法的階,就必須增加點(diǎn)數(shù).Heun(休恩)公式(3)c1

=

1/

4,

c2

=

3

/

4,

a2

=

b21

=

2

/

312=

yn

+

h

(

K

1

+

3

K

2

)

/

4

yn

+1

K=

f

(

xn

,

yn

)

K=

f

(

xn

+

2

h

,

y3

n+

2

h

K

)3

19.3.1

Runge-Kutta法的基本思想9.3龍格-庫塔(Runge-Kutta)法2021/7/1765(1)12=

yn

+

h

(

K

1

+

4

K

2

+

K

3

)

/

6

yn

+1

K=

f

(

xn

,

yn

)

K=

f

(

xn

+

1

h

,

y2

n+

1

h

K

)2

1=

f

(

xn

+

h

,

yn

-

h

K

1

+

2

h

K

2

)

K

3三階Kutta方法三階Heun方法(2)233=

y

k

+

h

(

k1

+

3

k

3

)

/

4

y

k

+1

k=

f

(

xk

,

y

k

)

1

k=

f

(

xk

+

1

h

,

y

+

1

hk

)3

k

3

1=

f

(

xk

+

2

h

,

y

k+

2

hk

2

)

k

39.3.1

Runge-Kutta法的基本思想常用的三階Runge-Kutta方法(r

=3):9.3龍格-庫塔(Runge-Kutta)法2021/7/176613nn=

yn

+

h

(

K

1

+

2

K

2

+

2

K

3

+

K

4

)

/

6標(biāo)準(zhǔn)(經(jīng)典)四階Runge-Kutta方法

yn

+1

K=

f

(

x

,

y

)n

K

2

=

f

(

xn

+

1

h

,

yn

+

1

h

K

1

)2

2

K=

f

(

x=

f

(

xn+

1

h

,

y

+

1

h

K

)2

n

2

2+

h

,

yn

+

h

K

3

)

K

4(1)9.3.1

Runge-Kutta法的基本思想常用的四階Runge-Kutta方法(r

=4):9.3龍格-庫塔(Runge-Kutta)法2021/7/179.3.1Runge-Kutta法的基本思想算法9.21)輸入a，b、f(x,y)，計(jì)算總步數(shù)N，初值y002)置

h

b

-

a

n=0，x=a，y=

y3)計(jì)算Nx=

x0+hK1=

f(x0，y0)K2=

f(x0+h/2，y0+

K1h/2)K3=

f(x0+h/2，y0+

K2h/2)9.2

改進(jìn)的(歐拉)Euler方法2021/7/179.2.1

梯形公式算法9.2K4=

f(x0+h，y0+

K3h)y=y0+h

(K1+2

K2+2

K3+

K4)/6輸出（x，y）4)如果n<N-1，置

n+1n，x

x0，y

y0，轉(zhuǎn)3)；否則，停機(jī)9.3龍格-庫塔(Runge-Kutta)法32021/7/17694221

2422

3nn

+1 2

-122-

2

y=

yn

+

h[

K1

+

(2

-

2

)

K

2

+

(2

+

2

)

K+

K

]

/

6

K1

=

f

(

xn

,

yn

)

K=

f

(

xn

+

1

h,

y

+

1

hK

)2

2

1

K=

f

(

xn

+

1

h,

y

+2

nhK

+

hK

)

3

K=

f

(

x

+

h,

y

-

2

hKn

n

2+

2+

2

hK

)9.3.1

Runge-Kutta法的基本思想(2)

Gill(吉爾)公式9.3龍格-庫塔(Runge-Kutta)法2021/7/1770預(yù)報點(diǎn)的個數(shù)r123456789r

≥10精度的階數(shù)123445667≤

r

-2一般情況下,四階Runge-Kutta方法已可滿足精度要求.9.3.1

Runge-Kutta法的基本思想注:

從理論上講,可以構(gòu)造任意高階的計(jì)算方法.但事實(shí)上,

精度的階數(shù)與預(yù)報點(diǎn)的個數(shù)之間并非等量關(guān)系.9.3龍格-庫塔(Runge-Kutta)法2021/7/1771例3:

用經(jīng)典Runge-Kutta方法求解下列初值問題(取

h

=

0.1)

y

￠=

2

x

+

y

, 0

￡

x

￡

1

y

(0)

=

1解:h

=

0.1,

x0

=

0,

y0

=

1,

xk

+1

=

xk

+

0.1

(k

=

0,1,

,

9)標(biāo)準(zhǔn)Runge-Kutta公式為:13n

nnn

2=

yn

+

0

.1 (

K

1

+

2

K

2

+

2

K

3+

K

4

)

/

6

yn

+1

K=

2

x

+

y

K

2

=

2(

xn

+

0

.05

)

+

yn

+

0

.05

K

1

K=

2(

x

+

0

.05

)

+

y

+

0

.05

K

K

4

=

2(

xn

+

0

.1)

+

yn

+

0

.1K

39.3龍格-庫塔(Runge-Kutta)法2021/7/1772計(jì)算結(jié)果見下表.為比較在相同計(jì)算量條件下近似解的精度,表中列出了Euler法(h=0.025)和改進(jìn)的Euler法(h=0.05)在相應(yīng)節(jié)點(diǎn)上的計(jì)算結(jié)果.注:用表中每種方法計(jì)算yi

都需要計(jì)算四次f

的值,即它們的計(jì)算量基本相等.xiEuler法h=0.025改進(jìn)Euler法h=0.05經(jīng)典R-K法h=0.1準(zhǔn)確解0.11.1114391.1153801.1155121.1155130.21.2552091.2639141.2642081.2642080.31.4346671.4490891.4495761.4495760.41.6535171.6747561.6754731.6754740.51.9158491.9451711.9461621.9461640.62.2261782.2650402.2663542.2663560.72.5894852.6395612.6412552.6412580.83.0112713.0744793.0766193.0766230.93.4976063.5761443.5788043.5788091.04.0551924.1515734.1548394.1548452021/7/17739.3龍格-庫塔(Runge-Kutta)法2021/7/17749.3.4

變步長的PK方法微分方程解函數(shù)y(x)的變化不均勻，等步長求數(shù)值解，可能造成某些點(diǎn)精度高，有些點(diǎn)精度低，為保證精度，必須取較小的步長。但隨著步長的縮小，在一定求解范圍內(nèi)所要完成的步數(shù)就增加了.步數(shù)的增加不但引起計(jì)算量的增大，而且可能導(dǎo)致舍入誤差的嚴(yán)重積累.因此同積分的數(shù)值計(jì)算一樣，微分方程的數(shù)值解法也有個選擇步長的問題.9.3龍格-庫塔(Runge-Kutta)法2021/7/1775變步長的PK方法在選擇步長時，需要考慮兩個問題：怎樣衡量和檢驗(yàn)計(jì)算結(jié)果的精度？如何依據(jù)所獲得的精度處理步長？2021/7/17769.3龍格-庫塔(Runge-Kutta)法O(h5)，故9.3.4

變步長的PK方法考察四階R-K公式，從節(jié)點(diǎn)xn出發(fā)，先以h為步長求出一個近似值，記為y(h

)，由于公式的局部截?cái)嗾`差為n+1y(xn+1

n+1)

-

y(

h

)?

ch5再求得一個近似值，每跨一步的局部截?cái)嗾`差然后將步長折半，即取為步長

2

，從xn跨兩步到xn+1，hhy(

)2n+1)5

，因此有2h是c(5

h

2

+1

n+1

h

y(xn

)

-

y

?

2c

2

9.3龍格-庫塔(Runge-Kutta)法即有9.3.4

變步長的PK方法兩式比可看到，步長折半后，誤差大約減少到1/16，由此易得下列事后估計(jì)式116

h

2

+1

n+1n+1n+1y(xn

)

-

y?

.y(x

)

-

y(

h)12021/7/1777].15

h

h

2

2

+1

n+1n+1(h)n+1y(xn

)

-

y

?[

y

-

y9.3龍格-庫塔(Runge-Kutta)法偏差9.3.4

變步長的PK方法因此，可以通過檢查步長，折半前后兩次計(jì)算結(jié)果的

h

2

n+1(h)n+1D

=

y

-

y直至Δ<ε為止,這時取最終得到的作為結(jié)果；來判定所選的步長是否合適，具體地說，將區(qū)分以下兩種情況處理：1.對于給定的精度ε,如果Δ>ε,我們反復(fù)將步長折半計(jì)算,(

)2021/7/17782hyn+19.3龍格-庫塔(Runge-Kutta)法2021/7/17799.3.4

變步長的PK方法2.如果Δ<ε，反復(fù)將步長作加倍計(jì)算，直至Δ>ε為止，此時前面一次的結(jié)果為yn+1，步長即為合適步長。這種通過加倍或折半處理步長的方法稱為變步長方法.表面上看，為了選擇步長，每一步的計(jì)算量增加了，但總體考慮往往是合算的.9.4

線性多步法2021/7/1780單步法優(yōu)點(diǎn)：只要給定初值，就可進(jìn)行計(jì)算。單步法缺點(diǎn)：計(jì)算yn+1只用到前一步的信息yn。提高精度需通過增加非節(jié)點(diǎn)處函數(shù)值的計(jì)算，如4階R-K公式需計(jì)算4個函數(shù)值。線性多步法基本思想：利用前面已知的信息yn，yn-1，…，yn-r構(gòu)造算法，提高精度。常用線性多步法，即計(jì)算公式中只出現(xiàn)yn+1，yn-1，…，yn-r及f(xn+1，yn+1)，f(xn，yn)，…，f(xn-r，yn-r)的一次項(xiàng)。9.4

線性多步法2021/7/1781r

ri=0

i=-1一般的線性多步法公式可表示為yn+1

=

ai

yn-i

+

h

bi

fn-i其中，ai，bi，均為常數(shù)，fk=f(xk，yk)（k=n+1，n，…，n-r)如果β-1=0，則上式稱為顯式r步法，這時yn+1可直接算出；如果β-1≠0,則上式稱為隱式r步法,求解時要用迭代法方可算出yn+1.公式中系數(shù)αi及βi可根據(jù)方法的局部截?cái)嗾`差及階確定。9.4

線性多步法2021/7/1782當(dāng)r=1時即為單步法此時若取a0=1，β-1

=0，b0=1，得到顯式Euler公式；0

-1

0若取a

=1，β

=

b

=

1

，得到梯形公式。2可以看出，線性多步法每次只需計(jì)算一次f(x，y)的值。構(gòu)造線性多步公式常用Taylor展開和數(shù)值積分的方法。9.4

線性多步法2021/7/1783將公式9.4.1

線性多步法公式的導(dǎo)出r

ryn+1

=

ai

yn-i

+

h

bi