計算方法4常微分方程數(shù)值解法課件

常微分方程的數(shù)值解法NumericalSolutionstoOrdinaryDifferentialEquations

常微分方程的數(shù)值解法對象一階常微分方程初值問題:一階常微分方程組初值問題:高階常微分方程初值問題:(4.1)對象一階常微分方程初值問題:一階常微分方程組初值問題:高階常一階常微分方程初值問題:實際工程技術(shù)、生產(chǎn)、科研上會出現(xiàn)大量的微分方程問題很難得到其解析解，有的甚至無法用解析表達(dá)式來表示，因此只能依賴于數(shù)值方法去獲得微分方程的數(shù)值解。一階常微分方程初值問題:實際工程技術(shù)、生產(chǎn)、科研上會出現(xiàn)大量ab

x0x1x2...xn-1xn用數(shù)值方法，求得y(x)在每個節(jié)點xk的值y(xk)的近似值，用yk表示，即yk

≈y(xk)，這樣y0,y1,...,yn稱為微分方程的數(shù)值解求y(x)————>求y0,y1,...,yn

？微分方程的數(shù)值解法：不求y=y(x)的精確表達(dá)式,而求離散點x0,x1,…xn處的函數(shù)值設(shè)(4.1)

的解y(x)的存在區(qū)間是[a,b]，初始點x0=a，取[a,b]內(nèi)的一系列節(jié)點x0,x1,...,xn。a=x0<x1<…<xn=b，一般采用等距步長abx0x1x2...思路計算過程：方法：采用步進(jìn)式和遞推法將[a,b]n等分,a=x0<x1<…<xn=b,步長h=(b-a)/n,xk=a+kh

怎樣建立遞推公式？Taylor公式數(shù)值積分法思路計算過程：方法：采用步進(jìn)式和遞推法將[a,b]n等分,4.1Euler

公式思想:用向前差商近似代替微商.(4.2)歐拉公式(EulerScheme）4.1Euler公式思想:用向前差商近似代替微商.(幾何意義 y(x)過點P0(x0,y0)且在任意點(x,y)的切線斜率為f(x,y)y(x)在點P0(x0,y0)的切線方程為：

y=y0+f(x0,y0)(x-x0)在切線上取點P1(x1,y1)

y1=y0+f(x0,y0)(x1-x0)y1正是Euler公式所求4.類似2，過P1以f(x1,y1)為斜率作y(x)的切線，在其上取點

P2(x2,y2)，依此類推…5.折線P0P1P2…Pn…作為曲線y(x)的近似

——歐拉折線法xp0p1p2p3p4x0x1x2x3x4yy(x)幾何意義 y(x)過點P0(x0,y0)且在任意點(x,y)思想:用向后差商近似代替微商.歐拉法（續(xù)）用隱式歐拉法，每一步都需解方程（或先解出yn+1的顯式表達(dá)式），但其穩(wěn)定性好。隱式歐拉公式(4.3)思想:用向后差商近似代替微商.歐拉法（續(xù)）用隱式歐拉法，計算方法4_常微分方程數(shù)值解法課件整體誤差ek=y(xk)-yk，下面對其加以分析y1=y0+hf(x0,y0)=1+0.1×(1-0/1)=1.1y2=y1+hf(x1,y1)=1.1+0.1×(1.1-2×0.1/1.1)=1.191818y3=y2+hf(x2,y2)=1.277438…其精確解為整體誤差ek=y(xk)-yk，下面對其加以分析y1=y0+計算方法4_常微分方程數(shù)值解法課件歐拉法（續(xù)）思想:用中心差商近似代替微商.注：計算時，先用歐拉法求出y1，以后再用二步歐拉法計算。二步歐拉法(4.4)歐拉法（續(xù)）思想:用中心差商近似代替微商.注：計算時，先歐拉法（續(xù)）公式單步否顯式否單步顯式單步隱式二步顯式截斷誤差y(xn+1)-yn+1

歐拉法（續(xù)）公式單步否顯式否單步顯式單步隱式二步顯式截斷誤差截斷誤差Def4.1設(shè)y(xn)

是(4.1)式的精確解，yn是(4.2)式歐拉法得到的近似解，稱y(xn)-yn為歐拉法的整體截斷誤差.

Def4.3若某算法的局部截斷誤差為O(hp+1),稱該算法有p階精度.Def4.2假設(shè)yn=y(xn)

,即第n步計算是精確的前提下,稱Rn+1=y(xn+1)-yn+1為歐拉法的局部截斷誤差.

截斷誤差Def4.1設(shè)y(xn)是(4.1)式的精確解，分析：證明其局部截斷誤差為O(h2),可通過Taylor展開式分析。證明：Euler

公式為令yn=y(xn)，下證：

y(xn+1)-yn+1=

O(h2)由y’(x)

=f(x,y(x))

定理4.4

歐拉法的精度是一階。分析：證明其局部截斷誤差為O(h2),可通過Taylor展開二步歐拉法的局部截斷誤差Def4.5假設(shè)yn=y(xn),yn－1=y(xn－1),稱Rn+1=y(xn+1)-yn+1為二步歐拉法的局部截斷誤差.

定理4.6

隱式歐拉法的精度是一階,二步歐拉法的精度是二階。證明：對二步歐拉法進(jìn)行證明，考慮其局部截斷誤差，令yn=y(xn),yn－1=y(xn－1),將上兩式左右兩端同時相減：∴二步歐拉法的局部截斷誤差為O(h3)，其精度是二階。二步歐拉法的局部截斷誤差Def4.5假設(shè)yn=y(xn)數(shù)值積分法對右端的定積分用數(shù)值積分公式求近似值：（1）用左矩形數(shù)值積分公式：（2）用梯形公式：——梯形公式梯形公式:將顯示歐拉公式,隱式歐拉公式平均可得梯形公式是隱式、單步公式,其精度為二階數(shù)值積分法對右端的定積分用數(shù)值積分公式求近似值：（1）用左矩證：令yn=y(xn),由Talor公式有分析：梯形公式是隱式公式，證明其局部截斷誤差為O(h3)

要用到二元函數(shù)的Taylor公式。f(xn+1,yn+1)=f(xn+1,y(xn+1)+(yn+1-y(xn+1))=f(xn+1,y(xn+1))+fy(xn+1,η)(yn+1-y(xn+1))，η∈(xn，xn+1)

=y’(xn+1)+fy(xn+1,η)(yn+1-y(xn+1))=y’(xn)+hy”(xn)+O(h2)+fy(xn+1,η)(yn+1-y(xn+1))=f(xn,yn)+hy”(xn)+fy(xn+1,η)(yn+1-y(xn+1))+O(h2)又y(xn+1)=y(xn+h)=y(xn)+hy’(xn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)/2+h[f(xn,yn)+hy”(xn)]/2+O(h3)定理4.7

梯形公式的精度是二階。證：令yn=y(xn),由Talor公式有分析：梯形公式是隱從而y(xn+1)=yn+1+hfy(xn+1,η)[y(xn+1)-yn+1]/2+O(h3)∴y(xn+1)-yn+1=

hfy(xn+1,η)[y(xn+1)-yn+1]/2+O(h3)∴y(xn+1)-yn+1=O(h3)/[1-hfy(xn+1,η)/2]=O(h3)∴梯形公式的截斷誤差為O(h3)，其精度是2階。f(xn+1,yn+1)=f(xn,yn)+hy”(xn)+fy(xn+1,η)(yn+1-y(xn+1))+O(h2)y(xn+1)=yn+hf(xn,yn)/2+h[f(xn,yn)+hy”(xn)]/2+O(h3)=yn+hf(xn,yn)/2+h[f(xn+1,yn+1)-fy(xn+1,η)(yn+1-y(xn+1))+O(h2)]/2+O(h3)=yn+h[f(xn,yn)+f(xn+1,yn+1)]/2+hfy(xn+1,η)(y(xn+1)-yn+1)/2+O(h3)從而y(xn+1)=yn+1+hfy(xn+1,η)[y(解:

取h=0.01x0=0y0=y(0)=1則

y(0.01)≈y1

f(x,y)=y，由梯形公式,基于穩(wěn)定性考慮解析解解:取h=0.01x0=0y0=y(0)=歐拉公式的比較歐拉法最簡單,精度低隱式歐拉法穩(wěn)定性好二步歐拉法顯式,但需要兩步初值,且第2個初值只能由其他方法給出,可能對后面的遞推精度有影響梯形公式法精度有所提高,但為隱式,需要迭代求解,計算量大歐拉公式的比較歐拉法最簡單,精度低隱式歐拉法穩(wěn)定性好二步歐拉4.2

改進(jìn)的Euler法

Euler公式計算量小，精度低梯形公式計算量大，精度高綜合兩個公式，提出預(yù)報—校正公式：

預(yù)報

校正——改進(jìn)的Euler法寫成嵌套公式：平均化形式：4.2改進(jìn)的Euler法Euler公式計算量小[例4.4]用改進(jìn)的Euler法解初值問題在區(qū)間[0,0.4]上，步長h=0.1的解，并比較與精確解(y=1/(1-x))的差異。解：Euler法的具體形式為：yn+1=yn+hyn2，改進(jìn)的Euler法的具體形式為：∵x0=0,h=0.1,則

x1=0.1,x2=0.2,x3=0.3,x4=0.4

計算y1:

yp=y0+0.1y02=1+0.1·12=1.1yc=1+0.1X1.12=1.121y1=(1.1+1.121)/2≈1.1118同樣可求y2,、y3、,y4[例4.4]用改進(jìn)的Euler法解初值問題在區(qū)間[0,0.注：（1）令y(xn)=yn，可推導(dǎo)改進(jìn)的Euler法的局部截斷誤差

為O(h3),具有二階精度。（2）

改進(jìn)的Euler法也可寫成如下平均化形式注：（2）改進(jìn)的Euler法也可寫成如下平均化形式4.3龍格—庫塔方法如何構(gòu)造高階的方法？(4.5)為了構(gòu)造函數(shù)φ使得(4.5)式成為高階方法，Taylor對于一般的顯式單步法若討論其精度（階數(shù)），即考慮4.3龍格—庫塔方法如何構(gòu)造高階的方法？(4.5)為了構(gòu)造令則有考慮到上述方法求高階微分，實際上不實用令則有考慮到上述方法求高階微分，實際上不實用改進(jìn)的Euler公式

:Euler公式：yn+1=yn+hf(xn,yn)寫成精度：一階

精度：二階

由Lagrange中值定理，稱為y(x)在[xn,xn+1]上的平均斜率取—Euler公式

取—改進(jìn)Euler公式

Euler公式用一個點的值k1作為k*的近似值，而改進(jìn)的Euler公式用二個點的值k1和k2的平均值作為k*近似值，改進(jìn)的Euler法比Euler法精度高。改進(jìn)的Euler公式:Euler公式：yn+1=yn+hf(4.6)Runge-Kutta法的思想：在［xn,xn+1］內(nèi)多預(yù)報幾個點的ki值并用其加權(quán)平均值作為k*近似值。而構(gòu)造出具有更高精度的計算公式。多個斜率加權(quán)平均可提高精度Runge-Kutta法一般形式：(4.6)Runge-Kutta法的思想：在［xn,xn+1以k1與k2的加權(quán)平均來近似k*，設(shè)其中c1,c2,α2,b21為待定參數(shù)。適當(dāng)選取參數(shù)，使(*)式的精度為二階，即使其局部截斷誤差為O(h3)令y(xn)=yn，由泰勒公式：二階龍格—庫塔方法(*)以k1與k2的加權(quán)平均來近似k*，設(shè)其中c1,c2,α2,b由多元函數(shù)的泰勒公式和(*)式則有

比較(a)與(b)要使y(xn+1)-yn+1=O(h3)(a)(b)由多元函數(shù)的泰勒公式和(*)式則有比較(a)與(b)要使上述方程組有四個未知量，只有三個方程，有無窮多組解。取任意一組解便得一種二階龍－庫公式。

當(dāng)c1=c2=1/2,a2=b21=1時二階Runge-Kutta公式為yn+1=yn+k1/2+k2/2k1=hf(xn,yn)k2=hf(xn+h,yn+k1)此即改進(jìn)Euler法

取c2=0，c2=1，a2=1/2，b21=1/2yn+1=yn+k2k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)此為中點法或變形的

Euler公式上述方程組有四個未知量，只有三個方程，有無窮多組解。當(dāng)c1三階龍格－庫塔法是用三個值k1,k2,k3的加權(quán)平均來近似k*，即有yn+1=yn+c1k1+c2k2+c3k3k1=hf(xn,yn)k2=hf(xn+a2h,yn+b21k1)k3=hf(xn+a3h,yn+b31k1+b32k2)要使其具有三階精度，必須使局部截斷誤差為O(h4)類似二階龍格－庫塔法的推導(dǎo)，c1,c2,c3,a2,a3,b21,b31,b32應(yīng)滿足c1+c2+c3=1a2=b21a3=b31+b32c2a2+c3a3=1/2c2a22+c3a32=1/3c3a32=1/6由該方程組任意解可得三階龍格-庫塔公式例：Kutta公式kn+1=yn+(k1+4k2+k3)/6k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)k3=hf(xn+h,yn-k1+2k2)

三階龍格－庫塔法是用三個值k1,k2,k3的加權(quán)平均來近似k類似可推出四階龍格-庫塔公式，常用的有例：經(jīng)典Runge-Kutta法

yn+1=yn+(k1+2k2+2k3+k4)/6k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)k3=hf(xn+h/2,yn+k2/2)k4=hf(xn+h,yn+k3)局部截斷誤差O(h5)

還有:Gill公式及m(m>4)階龍格－庫塔法。

m>4時：計算量太大，精確度不一定提高，有時會降低。類似可推出四階龍格-庫塔公式，常用的有yn+1=yn+(k1Gill公式節(jié)省存儲單元控制舍入誤差Gill公式節(jié)省存儲單元對于經(jīng)典的四階Runge-Kutta法給出如下算法：[算法4.2]求解：

dy/dx=f(x,y)a≤x≤by(a)=y0

Step

1:

輸入a,b,y0

及NStep

2:

(b-a)/N=>h,a=>x,y0=>yStep

3:

輸出

(x,y)Step4:ForI=1T0Nhf(x,y)=>k1hf(x+h/2,y+k1/2)=>k2hf(x+h/2,y+k2/2)=>k3hf(x+h,y＋k3)=>k4y+(k1+2k2+2k3+k4)/6=>yx+h=>x輸出(x,y)Step5:停止對于經(jīng)典的四階Runge-Kutta法給出如下算法：[算法4[例4.3]用四階經(jīng)典Runge－Kutta方法解初值問題：

(1)求，

[例4.3]用四階經(jīng)典Runge－Kutta方法解初值問題：(2)求，

(2)求，自適應(yīng)龍格－庫塔法用戶提出問題I：問題：①：如何判斷|y(xn)-yn|<ε∵精度值y(xn)未知。

②：如何取h=?解①：如用p階龍格－庫塔法計算，局部截斷誤差為O(hp+1)

y’=f(x,y)y(x0)=y0

要求誤差<ε＝10-8求數(shù)值解｛｝xnh/2xn+1h如

xn-----------------xn+1

令yn=y(xn)yn+1(h)

則y(xn+1)-yn+1(h)≈chp+1步長折半xnxn+h/2xn+1分兩步計算y(xn+1)的近似值yn+1(h/2)。則y(xn+1)-yn+1(h/2)≈2c(h/2)p+1

自適應(yīng)龍格－庫塔法用戶提出問題I：問題：①：如何判斷|y(定理：對于問題I若用P階龍格－庫塔法計算y(xn+1)在步長折半前后的近似值分別為yn+1(h),yn+1(h/2)則有誤差公式

注：10

誤差的事后估計法

20

停機(jī)準(zhǔn)則：△<ε(可保證|y(xn+1)-yn+1(h/2)|<ε)

解②：⑴h取大，局部截斷誤差chp+1大，不精確⑵h取小，運算量大（步多），舍入誤差積累大解決策略：變步長龍格－庫塔法

If（△>ε）將步長折半反復(fù)計算，直至△<ε為止，取h為最后一次的步長，yn+1為最后一次計算的結(jié)果。

Elseif(△<ε)將步長加倍反復(fù)計算，直至△>ε為止，最后一次運算的前一次計算結(jié)果即為所需。定理：對于問題I若用P階龍格－庫塔法計算y(xn+1)在步長4.5收斂與穩(wěn)定性對于常微分方程初值問題(4.1)求y=y(x)--------求y(xn)------------

求yn

xn=x0+nh離散化某種公式例：Euler法，改進(jìn)的Euler法，龍格－庫塔法都是單步法。數(shù)值解法：單步法：計算yn+1時只用到前一步的結(jié)果yn。顯式單步法：

yn+1=yn+hφ(xn,yn,h)

φ(x,y,h)為增量函數(shù)，它依賴于f,僅是xn,yn,h的函數(shù)Def:若某數(shù)值方法對于任意固定的xn=x0+nh,當(dāng)

h->0(n->∞)時，

yn->y(xn)，則稱該法是收斂的。即(xn=x0+nh為固定值)4.5收斂與穩(wěn)定性對于常微分方程初值問題(4.1)求y=改進(jìn)Euler法的收斂性

判定改進(jìn)Euler法的收斂性：

yn+1=yn+h[f(xn,yn)+f(xn+h,yn+hf(xn,yn))]/2

則φ(x,y,h)=[f(x,y)+f(x+h,y+hf(x,y))]/2若：①yo=y(x0)

②f(x,y)關(guān)于y滿足李--條件:則：

限定h≤h0，則即Φ(x,y,h)

滿足李普希茲條件,由Th改進(jìn)的Euler法收斂同樣可驗證，若f(x,y)關(guān)于y滿足李普希茲條件,且y0=y(x0)時，Euler法，標(biāo)準(zhǔn)四階龍格－庫法也收斂。改進(jìn)Euler法的收斂性判定改進(jìn)Euler法的收斂性：若：4.5.2

穩(wěn)定性在討論收斂性時，一般認(rèn)可：數(shù)值方法本身計算過程是準(zhǔn)確的。實際上，并非如此：①

初始值y0有誤差δ0＝y(tǒng)0-y(x0).②

計算的每一步有舍入誤差。這些初始誤差在計算過程的傳播中，是逐步衰減的，還是惡性增長，這就是穩(wěn)定性問題。4.5.2穩(wěn)定性在討論收斂性時，一般認(rèn)可：數(shù)值方法本身計算Def4.1若一種數(shù)值方法在節(jié)點xn處的數(shù)值解yn的擾動δn≠0，而在以后的各節(jié)點值ym(m>n)上有擾動δm.當(dāng)|δm|≤|δn|，(m=n+1,n+2,……)，則稱該數(shù)值算法是穩(wěn)定的。分析某算法的數(shù)值穩(wěn)定性，通?？疾炷Ｐ头匠?/p>

y’=λy，(λ<0)Def4.1：設(shè)在節(jié)點xn處用數(shù)值法得到的理想數(shù)值解為yn，而實際計算得到的近似值為，稱為第n步數(shù)值解的擾動Def4.1若一種數(shù)值方法在節(jié)點xn處的數(shù)值解yn的擾動δnEuler法的穩(wěn)定性Euler法：yn+1=yn+hf(xn,yn)考察模型方程

y’=λy，(λ<0)

即yn+1=(1+hλ)yn假設(shè)在節(jié)點值yn上有擾動δn，在yn+1上有擾動δn+1，且δn+1僅由δn引起(計算過程不再引進(jìn)新的誤差)Euler法穩(wěn)定Euler法的穩(wěn)定的條件為0<h≤-2/λEuler法的穩(wěn)定性Euler法：yn+1=yn+hf(xn隱式Euler法穩(wěn)定性隱式Euler法：yn+1=yn+hf(xn+1,yn+1)對于模型方程y’=λy(λ<0)則yn+1=yn+hλyn+1=>

yn+1=yn/(1-hλ)yn+1的擾動∵λ<0，∴h>0

上式均成立，隱式Euler法無條件穩(wěn)定210ReImg隱式Euler法穩(wěn)定性隱式Euler法：yn+1=yn+hf梯形公式穩(wěn)定性梯形公式

yn+1=yn+h[f(xn,yn)+f(xn+1,yn+1)]/2，設(shè)模型方程為y’=λy(λ<0)則yn+1=yn+h[λyn+λyn+1]/2

λ<0上式對任意h>0時恒成立

梯形公式恒穩(wěn)定。

yn+1的擾動梯形公式穩(wěn)定性梯形公式y(tǒng)n+1=yn+h[f(xn,yn3階Runge－Kutta顯式

1~4階方法的絕對穩(wěn)定區(qū)域為k=1k=2k=3k=4-1-2-3---123ReImg3階Runge－Kutta顯式1~4階方法的絕對Taylor公式公式單步否顯式否精度單步顯式一階單步隱式一階二步顯式二階數(shù)值積分法單步隱式二階單步顯式二階Runge－Kutta方法局部截斷誤差整體截斷誤差收斂與穩(wěn)定性Taylor公式公式單步否顯式否精度單步顯式一階單步隱式一4.6

多步法某一步解和前若干步解的值有關(guān)m步線性多步法的一般形式其中a0,…,am-1,b0,…,bm為和h無關(guān)的常數(shù)初值為y0,y1,…,ym-1若bm≠0，方法是隱式的(implicit)否則是顯式的(explicit)(4.7)4.6多步法某一步解和前若干步解的值有關(guān)m步線性多步法的一Adams－Bashforth四階方法Adams－Moulton四階方法Adams－Bashforth四階方法Adams－Moult多步法的截斷誤差多步法的截斷誤差4.7

常微分方程組和高階微分方程的數(shù)值解法形如歐拉方法的計算公式隱式歐拉方法的計算公式4.7常微分方程組和高階微分方程的數(shù)值解法形如歐拉方法的計梯形公式Adams公式R-K方法梯形公式Adams公式R-K方法對于高階微分方程，總可以化成方程組的形式，例如可以化成對于高階微分方程，總可以化成方程組的形式，例如可以化成例令y’(x)=z，則有用4階RK方法，有例令y’(x)=z，則有用4階RK方法，有計算方法4_常微分方程數(shù)值解法課件4.8

常微分方程邊值問題的數(shù)值解法形如第一類邊界條件(4.8)第二類邊界條件第三類邊界條件4.8常微分方程邊值問題的數(shù)值解法形如第一類邊界條件(4.定理假定問題（4.8）式中的函數(shù)f在集合上連續(xù)，且偏導(dǎo)數(shù)fy和fy’也在D上連續(xù)。若有則邊值問題有唯一解。定理假定問題（4.8）式中的函數(shù)f在集合上連續(xù)，且偏導(dǎo)數(shù)f例邊值問題有因此邊值問題有唯一解。定理若線性邊值問題滿足則邊值問題有唯一解。例邊值問題有因此邊值問題有唯一解。定理若線性邊值問題滿足打靶法（ShootingMethod）其中，y1(x)是初值問題的解的解，y2(x)是初值問題打靶法（ShootingMethod）其中，y1(x)是初有限差分法(Finite-DifferenceMethods)以差商代替導(dǎo)數(shù)，將微分方程離散成為差分方程用Taylor展開方法有限差分法(Finite-DifferenceMethod兩式相加得到由中值定理得線性邊值問題可寫成兩式相加得到由中值定理得線性邊值問題可寫成然后可得到有限差分公式上式可寫成寫成矩陣形式然后可得到有限差分公式上式可寫成寫成矩陣形式其中其中常微分方程的數(shù)值解法NumericalSolutionstoOrdinaryDifferentialEquations

常微分方程的數(shù)值解法對象一階常微分方程初值問題:一階常微分方程組初值問題:高階常微分方程初值問題:(4.1)對象一階常微分方程初值問題:一階常微分方程組初值問題:高階常一階常微分方程初值問題:實際工程技術(shù)、生產(chǎn)、科研上會出現(xiàn)大量的微分方程問題很難得到其解析解，有的甚至無法用解析表達(dá)式來表示，因此只能依賴于數(shù)值方法去獲得微分方程的數(shù)值解。一階常微分方程初值問題:實際工程技術(shù)、生產(chǎn)、科研上會出現(xiàn)大量ab

x0x1x2...xn-1xn用數(shù)值方法，求得y(x)在每個節(jié)點xk的值y(xk)的近似值，用yk表示，即yk

≈y(xk)，這樣y0,y1,...,yn稱為微分方程的數(shù)值解求y(x)————>求y0,y1,...,yn

？微分方程的數(shù)值解法：不求y=y(x)的精確表達(dá)式,而求離散點x0,x1,…xn處的函數(shù)值設(shè)(4.1)

的解y(x)的存在區(qū)間是[a,b]，初始點x0=a，取[a,b]內(nèi)的一系列節(jié)點x0,x1,...,xn。a=x0<x1<…<xn=b，一般采用等距步長abx0x1x2...思路計算過程：方法：采用步進(jìn)式和遞推法將[a,b]n等分,a=x0<x1<…<xn=b,步長h=(b-a)/n,xk=a+kh

怎樣建立遞推公式？Taylor公式數(shù)值積分法思路計算過程：方法：采用步進(jìn)式和遞推法將[a,b]n等分,4.1Euler

公式思想:用向前差商近似代替微商.(4.2)歐拉公式(EulerScheme）4.1Euler公式思想:用向前差商近似代替微商.(幾何意義 y(x)過點P0(x0,y0)且在任意點(x,y)的切線斜率為f(x,y)y(x)在點P0(x0,y0)的切線方程為：

y=y0+f(x0,y0)(x-x0)在切線上取點P1(x1,y1)

y1=y0+f(x0,y0)(x1-x0)y1正是Euler公式所求4.類似2，過P1以f(x1,y1)為斜率作y(x)的切線，在其上取點

P2(x2,y2)，依此類推…5.折線P0P1P2…Pn…作為曲線y(x)的近似

——歐拉折線法xp0p1p2p3p4x0x1x2x3x4yy(x)幾何意義 y(x)過點P0(x0,y0)且在任意點(x,y)思想:用向后差商近似代替微商.歐拉法（續(xù)）用隱式歐拉法，每一步都需解方程（或先解出yn+1的顯式表達(dá)式），但其穩(wěn)定性好。隱式歐拉公式(4.3)思想:用向后差商近似代替微商.歐拉法（續(xù)）用隱式歐拉法，計算方法4_常微分方程數(shù)值解法課件整體誤差ek=y(xk)-yk，下面對其加以分析y1=y0+hf(x0,y0)=1+0.1×(1-0/1)=1.1y2=y1+hf(x1,y1)=1.1+0.1×(1.1-2×0.1/1.1)=1.191818y3=y2+hf(x2,y2)=1.277438…其精確解為整體誤差ek=y(xk)-yk，下面對其加以分析y1=y0+計算方法4_常微分方程數(shù)值解法課件歐拉法（續(xù)）思想:用中心差商近似代替微商.注：計算時，先用歐拉法求出y1，以后再用二步歐拉法計算。二步歐拉法(4.4)歐拉法（續(xù)）思想:用中心差商近似代替微商.注：計算時，先歐拉法（續(xù)）公式單步否顯式否單步顯式單步隱式二步顯式截斷誤差y(xn+1)-yn+1

歐拉法（續(xù)）公式單步否顯式否單步顯式單步隱式二步顯式截斷誤差截斷誤差Def4.1設(shè)y(xn)

是(4.1)式的精確解，yn是(4.2)式歐拉法得到的近似解，稱y(xn)-yn為歐拉法的整體截斷誤差.

Def4.3若某算法的局部截斷誤差為O(hp+1),稱該算法有p階精度.Def4.2假設(shè)yn=y(xn)

,即第n步計算是精確的前提下,稱Rn+1=y(xn+1)-yn+1為歐拉法的局部截斷誤差.

截斷誤差Def4.1設(shè)y(xn)是(4.1)式的精確解，分析：證明其局部截斷誤差為O(h2),可通過Taylor展開式分析。證明：Euler

公式為令yn=y(xn)，下證：

y(xn+1)-yn+1=

O(h2)由y’(x)

=f(x,y(x))

定理4.4

歐拉法的精度是一階。分析：證明其局部截斷誤差為O(h2),可通過Taylor展開二步歐拉法的局部截斷誤差Def4.5假設(shè)yn=y(xn),yn－1=y(xn－1),稱Rn+1=y(xn+1)-yn+1為二步歐拉法的局部截斷誤差.

定理4.6

隱式歐拉法的精度是一階,二步歐拉法的精度是二階。證明：對二步歐拉法進(jìn)行證明，考慮其局部截斷誤差，令yn=y(xn),yn－1=y(xn－1),將上兩式左右兩端同時相減：∴二步歐拉法的局部截斷誤差為O(h3)，其精度是二階。二步歐拉法的局部截斷誤差Def4.5假設(shè)yn=y(xn)數(shù)值積分法對右端的定積分用數(shù)值積分公式求近似值：（1）用左矩形數(shù)值積分公式：（2）用梯形公式：——梯形公式梯形公式:將顯示歐拉公式,隱式歐拉公式平均可得梯形公式是隱式、單步公式,其精度為二階數(shù)值積分法對右端的定積分用數(shù)值積分公式求近似值：（1）用左矩證：令yn=y(xn),由Talor公式有分析：梯形公式是隱式公式，證明其局部截斷誤差為O(h3)

要用到二元函數(shù)的Taylor公式。f(xn+1,yn+1)=f(xn+1,y(xn+1)+(yn+1-y(xn+1))=f(xn+1,y(xn+1))+fy(xn+1,η)(yn+1-y(xn+1))，η∈(xn，xn+1)

=y’(xn+1)+fy(xn+1,η)(yn+1-y(xn+1))=y’(xn)+hy”(xn)+O(h2)+fy(xn+1,η)(yn+1-y(xn+1))=f(xn,yn)+hy”(xn)+fy(xn+1,η)(yn+1-y(xn+1))+O(h2)又y(xn+1)=y(xn+h)=y(xn)+hy’(xn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)/2+h[f(xn,yn)+hy”(xn)]/2+O(h3)定理4.7

梯形公式的精度是二階。證：令yn=y(xn),由Talor公式有分析：梯形公式是隱從而y(xn+1)=yn+1+hfy(xn+1,η)[y(xn+1)-yn+1]/2+O(h3)∴y(xn+1)-yn+1=

hfy(xn+1,η)[y(xn+1)-yn+1]/2+O(h3)∴y(xn+1)-yn+1=O(h3)/[1-hfy(xn+1,η)/2]=O(h3)∴梯形公式的截斷誤差為O(h3)，其精度是2階。f(xn+1,yn+1)=f(xn,yn)+hy”(xn)+fy(xn+1,η)(yn+1-y(xn+1))+O(h2)y(xn+1)=yn+hf(xn,yn)/2+h[f(xn,yn)+hy”(xn)]/2+O(h3)=yn+hf(xn,yn)/2+h[f(xn+1,yn+1)-fy(xn+1,η)(yn+1-y(xn+1))+O(h2)]/2+O(h3)=yn+h[f(xn,yn)+f(xn+1,yn+1)]/2+hfy(xn+1,η)(y(xn+1)-yn+1)/2+O(h3)從而y(xn+1)=yn+1+hfy(xn+1,η)[y(解:

取h=0.01x0=0y0=y(0)=1則

y(0.01)≈y1

f(x,y)=y，由梯形公式,基于穩(wěn)定性考慮解析解解:取h=0.01x0=0y0=y(0)=歐拉公式的比較歐拉法最簡單,精度低隱式歐拉法穩(wěn)定性好二步歐拉法顯式,但需要兩步初值,且第2個初值只能由其他方法給出,可能對后面的遞推精度有影響梯形公式法精度有所提高,但為隱式,需要迭代求解,計算量大歐拉公式的比較歐拉法最簡單,精度低隱式歐拉法穩(wěn)定性好二步歐拉4.2

改進(jìn)的Euler法

Euler公式計算量小，精度低梯形公式計算量大，精度高綜合兩個公式，提出預(yù)報—校正公式：

預(yù)報

校正——改進(jìn)的Euler法寫成嵌套公式：平均化形式：4.2改進(jìn)的Euler法Euler公式計算量小[例4.4]用改進(jìn)的Euler法解初值問題在區(qū)間[0,0.4]上，步長h=0.1的解，并比較與精確解(y=1/(1-x))的差異。解：Euler法的具體形式為：yn+1=yn+hyn2，改進(jìn)的Euler法的具體形式為：∵x0=0,h=0.1,則

x1=0.1,x2=0.2,x3=0.3,x4=0.4

計算y1:

yp=y0+0.1y02=1+0.1·12=1.1yc=1+0.1X1.12=1.121y1=(1.1+1.121)/2≈1.1118同樣可求y2,、y3、,y4[例4.4]用改進(jìn)的Euler法解初值問題在區(qū)間[0,0.注：（1）令y(xn)=yn，可推導(dǎo)改進(jìn)的Euler法的局部截斷誤差

為O(h3),具有二階精度。（2）

改進(jìn)的Euler法也可寫成如下平均化形式注：（2）改進(jìn)的Euler法也可寫成如下平均化形式4.3龍格—庫塔方法如何構(gòu)造高階的方法？(4.5)為了構(gòu)造函數(shù)φ使得(4.5)式成為高階方法，Taylor對于一般的顯式單步法若討論其精度（階數(shù)），即考慮4.3龍格—庫塔方法如何構(gòu)造高階的方法？(4.5)為了構(gòu)造令則有考慮到上述方法求高階微分，實際上不實用令則有考慮到上述方法求高階微分，實際上不實用改進(jìn)的Euler公式

:Euler公式：yn+1=yn+hf(xn,yn)寫成精度：一階

精度：二階

由Lagrange中值定理，稱為y(x)在[xn,xn+1]上的平均斜率取—Euler公式

取—改進(jìn)Euler公式

Euler公式用一個點的值k1作為k*的近似值，而改進(jìn)的Euler公式用二個點的值k1和k2的平均值作為k*近似值，改進(jìn)的Euler法比Euler法精度高。改進(jìn)的Euler公式:Euler公式：yn+1=yn+hf(4.6)Runge-Kutta法的思想：在［xn,xn+1］內(nèi)多預(yù)報幾個點的ki值并用其加權(quán)平均值作為k*近似值。而構(gòu)造出具有更高精度的計算公式。多個斜率加權(quán)平均可提高精度Runge-Kutta法一般形式：(4.6)Runge-Kutta法的思想：在［xn,xn+1以k1與k2的加權(quán)平均來近似k*，設(shè)其中c1,c2,α2,b21為待定參數(shù)。適當(dāng)選取參數(shù)，使(*)式的精度為二階，即使其局部截斷誤差為O(h3)令y(xn)=yn，由泰勒公式：二階龍格—庫塔方法(*)以k1與k2的加權(quán)平均來近似k*，設(shè)其中c1,c2,α2,b由多元函數(shù)的泰勒公式和(*)式則有

比較(a)與(b)要使y(xn+1)-yn+1=O(h3)(a)(b)由多元函數(shù)的泰勒公式和(*)式則有比較(a)與(b)要使上述方程組有四個未知量，只有三個方程，有無窮多組解。取任意一組解便得一種二階龍－庫公式。

當(dāng)c1=c2=1/2,a2=b21=1時二階Runge-Kutta公式為yn+1=yn+k1/2+k2/2k1=hf(xn,yn)k2=hf(xn+h,yn+k1)此即改進(jìn)Euler法

取c2=0，c2=1，a2=1/2，b21=1/2yn+1=yn+k2k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)此為中點法或變形的

Euler公式上述方程組有四個未知量，只有三個方程，有無窮多組解。當(dāng)c1三階龍格－庫塔法是用三個值k1,k2,k3的加權(quán)平均來近似k*，即有yn+1=yn+c1k1+c2k2+c3k3k1=hf(xn,yn)k2=hf(xn+a2h,yn+b21k1)k3=hf(xn+a3h,yn+b31k1+b32k2)要使其具有三階精度，必須使局部截斷誤差為O(h4)類似二階龍格－庫塔法的推導(dǎo)，c1,c2,c3,a2,a3,b21,b31,b32應(yīng)滿足c1+c2+c3=1a2=b21a3=b31+b32c2a2+c3a3=1/2c2a22+c3a32=1/3c3a32=1/6由該方程組任意解可得三階龍格-庫塔公式例：Kutta公式kn+1=yn+(k1+4k2+k3)/6k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)k3=hf(xn+h,yn-k1+2k2)

三階龍格－庫塔法是用三個值k1,k2,k3的加權(quán)平均來近似k類似可推出四階龍格-庫塔公式，常用的有例：經(jīng)典Runge-Kutta法

yn+1=yn+(k1+2k2+2k3+k4)/6k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)k3=hf(xn+h/2,yn+k2/2)k4=hf(xn+h,yn+k3)局部截斷誤差O(h5)

還有:Gill公式及m(m>4)階龍格－庫塔法。

m>4時：計算量太大，精確度不一定提高，有時會降低。類似可推出四階龍格-庫塔公式，常用的有yn+1=yn+(k1Gill公式節(jié)省存儲單元控制舍入誤差Gill公式節(jié)省存儲單元對于經(jīng)典的四階Runge-Kutta法給出如下算法：[算法4.2]求解：

dy/dx=f(x,y)a≤x≤by(a)=y0

Step

1:

輸入a,b,y0

及NStep

2:

(b-a)/N=>h,a=>x,y0=>yStep

3:

輸出

(x,y)Step4:ForI=1T0Nhf(x,y)=>k1hf(x+h/2,y+k1/2)=>k2hf(x+h/2,y+k2/2)=>k3hf(x+h,y＋k3)=>k4y+(k1+2k2+2k3+k4)/6=>yx+h=>x輸出(x,y)Step5:停止對于經(jīng)典的四階Runge-Kutta法給出如下算法：[算法4[例4.3]用四階經(jīng)典Runge－Kutta方法解初值問題：

(1)求，

[例4.3]用四階經(jīng)典Runge－Kutta方法解初值問題：(2)求，

(2)求，自適應(yīng)龍格－庫塔法用戶提出問題I：問題：①：如何判斷|y(xn)-yn|<ε∵精度值y(xn)未知。

②：如何取h=?解①：如用p階龍格－庫塔法計算，局部截斷誤差為O(hp+1)

y’=f(x,y)y(x0)=y0

要求誤差<ε＝10-8求數(shù)值解｛｝xnh/2xn+1h如

xn-----------------xn+1

令yn=y(xn)yn+1(h)

則y(xn+1)-yn+1(h)≈chp+1步長折半xnxn+h/2xn+1分兩步計算y(xn+1)的近似值yn+1(h/2)。則y(xn+1)-yn+1(h/2)≈2c(h/2)p+1

自適應(yīng)龍格－庫塔法用戶提出問題I：問題：①：如何判斷|y(定理：對于問題I若用P階龍格－庫塔法計算y(xn+1)在步長折半前后的近似值分別為yn+1(h),yn+1(h/2)則有誤差公式

注：10

誤差的事后估計法

20

停機(jī)準(zhǔn)則：△<ε(可保證|y(xn+1)-yn+1(h/2)|<ε)

解②：⑴h取大，局部截斷誤差chp+1大，不精確⑵h取小，運算量大（步多），舍入誤差積累大解決策略：變步長龍格－庫塔法

If（△>ε）將步長折半反復(fù)計算，直至△<ε為止，取h為最后一次的步長，yn+1為最后一次計算的結(jié)果。

Elseif(△<ε)將步長加倍反復(fù)計算，直至△>ε為止，最后一次運算的前一次計算結(jié)果即為所需。定理：對于問題I若用P階龍格－庫塔法計算y(xn+1)在步長4.5收斂與穩(wěn)定性對于常微分方程初值問題(4.1)求y=y(x)--------求y(xn)------------

求yn

xn=x0+nh離散化某種公式例：Euler法，改進(jìn)的Euler法，龍格－庫塔法都是單步法。數(shù)值解法：單步法：計算yn+1時只用到前一步的結(jié)果yn。顯式單步法：

yn+1=yn+hφ(xn,yn,h)

φ(x,y,h)為增量函數(shù)，它依賴于f,僅是xn,yn,h的函數(shù)Def:若某數(shù)值方法對于任意固定的xn=x0+nh,當(dāng)

h->0(n->∞)時，

yn->y(xn)，則稱該法是收斂的。即(xn=x0+nh為固定值)4.5收斂與穩(wěn)定性對于常微分方程初值問題(4.1)求y=改進(jìn)Euler法的收斂性

判定改進(jìn)Euler法的收斂性：

yn+1=yn+h[f(xn,yn)+f(xn+h,yn+hf(xn,yn))]/2

則φ(x,y,h)=[f(x,y)+f(x+h,y+hf(x,y))]/2若：①yo=y(x0)

②f(x,y)關(guān)于y滿足李--條件:則：

限定h≤h0，則即Φ(x,y,h)

滿足李普希茲條件,由Th改進(jìn)的Euler法收斂同樣可驗證，若f(x,y)關(guān)于y滿足李普希茲條件,且y0=y(x0)時，Euler法，標(biāo)準(zhǔn)四階龍格－庫法也收斂。改進(jìn)Euler法的收斂性判定改進(jìn)Euler法的收斂性：若：4.5.2

穩(wěn)定性在討論收斂性時，一般認(rèn)可：數(shù)值方法本身計算過程是準(zhǔn)確的。實際上，并非如此：①

初始值y0有誤差δ0＝