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

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

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

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

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

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

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

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

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

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

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

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

——?dú)W拉折線法xp0p1p2p3p4x0x1x2x3x4yy(x)幾何意義 y(x)過點(diǎn)P0(x0,y0)且在任意點(diǎn)(x,y)思想:用向后差商近似代替微商.歐拉法（續(xù)）用隱式歐拉法，每一步都需解方程（或先解出yn+1的顯式表達(dá)式），但其穩(wěn)定性好。隱式歐拉公式(4.3)思想:用向后差商近似代替微商.歐拉法（續(xù)）用隱式歐拉法，計(jì)算方法4_常微分方程數(shù)值解法課件整體誤差ek=y(xk)-yk，下面對(duì)其加以分析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，下面對(duì)其加以分析y1=y0+計(jì)算方法4_常微分方程數(shù)值解法課件歐拉法（續(xù)）思想:用中心差商近似代替微商.注：計(jì)算時(shí)，先用歐拉法求出y1，以后再用二步歐拉法計(jì)算。二步歐拉法(4.4)歐拉法（續(xù)）思想:用中心差商近似代替微商.注：計(jì)算時(shí)，先歐拉法（續(xù)）公式單步否顯式否單步顯式單步隱式二步顯式截?cái)嗾`差y(xn+1)-yn+1

歐拉法（續(xù)）公式單步否顯式否單步顯式單步隱式二步顯式截?cái)嗾`差截?cái)嗾`差Def4.1設(shè)y(xn)

是(4.1)式的精確解，yn是(4.2)式歐拉法得到的近似解，稱y(xn)-yn為歐拉法的整體截?cái)嗾`差.

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

,即第n步計(jì)算是精確的前提下,稱Rn+1=y(xn+1)-yn+1為歐拉法的局部截?cái)嗾`差.

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

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

y(xn+1)-yn+1=

O(h2)由y’(x)

=f(x,y(x))

定理4.4

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

定理4.6

隱式歐拉法的精度是一階,二步歐拉法的精度是二階。證明：對(duì)二步歐拉法進(jìn)行證明，考慮其局部截?cái)嗾`差，令yn=y(xn),yn－1=y(xn－1),將上兩式左右兩端同時(shí)相減：∴二步歐拉法的局部截?cái)嗾`差為O(h3)，其精度是二階。二步歐拉法的局部截?cái)嗾`差Def4.5假設(shè)yn=y(xn)數(shù)值積分法對(duì)右端的定積分用數(shù)值積分公式求近似值：（1）用左矩形數(shù)值積分公式：（2）用梯形公式：——梯形公式梯形公式:將顯示歐拉公式,隱式歐拉公式平均可得梯形公式是隱式、單步公式,其精度為二階數(shù)值積分法對(duì)右端的定積分用數(shù)值積分公式求近似值：（1）用左矩證：令yn=y(xn),由Talor公式有分析：梯形公式是隱式公式，證明其局部截?cái)嗾`差為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)∴梯形公式的截?cái)嗾`差為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個(gè)初值只能由其他方法給出,可能對(duì)后面的遞推精度有影響梯形公式法精度有所提高,但為隱式,需要迭代求解,計(jì)算量大歐拉公式的比較歐拉法最簡單,精度低隱式歐拉法穩(wěn)定性好二步歐拉4.2

改進(jìn)的Euler法

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

預(yù)報(bào)

校正——改進(jìn)的Euler法寫成嵌套公式：平均化形式：4.2改進(jìn)的Euler法Euler公式計(jì)算量小[例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

計(jì)算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法的局部截?cái)嗾`差

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

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

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

精度：二階

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

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

Euler公式用一個(gè)點(diǎn)的值k1作為k*的近似值，而改進(jìn)的Euler公式用二個(gè)點(diǎn)的值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ù)報(bào)幾個(gè)點(diǎn)的ki值并用其加權(quán)平均值作為k*近似值。而構(gòu)造出具有更高精度的計(jì)算公式。多個(gè)斜率加權(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ù)，使(*)式的精度為二階，即使其局部截?cái)嗾`差為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)要使上述方程組有四個(gè)未知量，只有三個(gè)方程，有無窮多組解。取任意一組解便得一種二階龍－庫公式。

當(dāng)c1=c2=1/2,a2=b21=1時(shí)二階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)此為中點(diǎn)法或變形的

Euler公式上述方程組有四個(gè)未知量，只有三個(gè)方程，有無窮多組解。當(dāng)c1三階龍格－庫塔法是用三個(gè)值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)要使其具有三階精度，必須使局部截?cái)嗾`差為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)

三階龍格－庫塔法是用三個(gè)值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)局部截?cái)嗾`差O(h5)

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

m>4時(shí)：計(jì)算量太大，精確度不一定提高，有時(shí)會(huì)降低。類似可推出四階龍格-庫塔公式，常用的有yn+1=yn+(k1Gill公式節(jié)省存儲(chǔ)單元控制舍入誤差Gill公式節(jié)省存儲(chǔ)單元對(duì)于經(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:停止對(duì)于經(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階龍格－庫塔法計(jì)算，局部截?cái)嗾`差為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分兩步計(jì)算y(xn+1)的近似值yn+1(h/2)。則y(xn+1)-yn+1(h/2)≈2c(h/2)p+1

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

注：10

誤差的事后估計(jì)法

20

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

解②：⑴h取大，局部截?cái)嗾`差chp+1大，不精確⑵h取小，運(yùn)算量大（步多），舍入誤差積累大解決策略：變步長龍格－庫塔法

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

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

求yn

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

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

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

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

yn->y(xn)，則稱該法是收斂的。即(xn=x0+nh為固定值)4.5收斂與穩(wěn)定性對(duì)于常微分方程初值問題(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法收斂同樣可驗(yàn)證，若f(x,y)關(guān)于y滿足李普希茲條件,且y0=y(x0)時(shí)，Euler法，標(biāo)準(zhǔn)四階龍格－庫法也收斂。改進(jìn)Euler法的收斂性判定改進(jìn)Euler法的收斂性：若：4.5.2

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

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

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

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

y’=λy，(λ<0)

即yn+1=(1+hλ)yn假設(shè)在節(jié)點(diǎn)值yn上有擾動(dòng)δn，在yn+1上有擾動(dòng)δn+1，且δn+1僅由δn引起(計(jì)算過程不再引進(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)對(duì)于模型方程y’=λy(λ<0)則yn+1=yn+hλyn+1=>

yn+1=yn/(1-hλ)yn+1的擾動(dòng)∵λ<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上式對(duì)任意h>0時(shí)恒成立

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

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

1~4階方法的絕對(duì)穩(wěn)定區(qū)域?yàn)閗=1k=2k=3k=4-1-2-3---123ReImg3階Runge－Kutta顯式1~4階方法的絕對(duì)Taylor公式公式單步否顯式否精度單步顯式一階單步隱式一階二步顯式二階數(shù)值積分法單步隱式二階單步顯式二階Runge－Kutta方法局部截?cái)嗾`差整體截?cái)嗾`差收斂與穩(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多步法的截?cái)嗾`差多步法的截?cái)嗾`差4.7

常微分方程組和高階微分方程的數(shù)值解法形如歐拉方法的計(jì)算公式隱式歐拉方法的計(jì)算公式4.7常微分方程組和高階微分方程的數(shù)值解法形如歐拉方法的計(jì)梯形公式Adams公式R-K方法梯形公式Adams公式R-K方法對(duì)于高階微分方程，總可以化成方程組的形式，例如可以化成對(duì)于高階微分方程，總可以化成方程組的形式，例如可以化成例令y’(x)=z，則有用4階RK方法，有例令y’(x)=z，則有用4階RK方法，有計(jì)算方法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ù)值解法對(duì)象一階常微分方程初值問題:一階常微分方程組初值問題:高階常微分方程初值問題:(4.1)對(duì)象一階常微分方程初值問題:一階常微分方程組初值問題:高階常一階常微分方程初值問題:實(shí)際工程技術(shù)、生產(chǎn)、科研上會(huì)出現(xiàn)大量的微分方程問題很難得到其解析解，有的甚至無法用解析表達(dá)式來表示，因此只能依賴于數(shù)值方法去獲得微分方程的數(shù)值解。一階常微分方程初值問題:實(shí)際工程技術(shù)、生產(chǎn)、科研上會(huì)出現(xiàn)大量ab

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

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

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

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

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

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

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

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

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

——?dú)W拉折線法xp0p1p2p3p4x0x1x2x3x4yy(x)幾何意義 y(x)過點(diǎn)P0(x0,y0)且在任意點(diǎn)(x,y)思想:用向后差商近似代替微商.歐拉法（續(xù)）用隱式歐拉法，每一步都需解方程（或先解出yn+1的顯式表達(dá)式），但其穩(wěn)定性好。隱式歐拉公式(4.3)思想:用向后差商近似代替微商.歐拉法（續(xù)）用隱式歐拉法，計(jì)算方法4_常微分方程數(shù)值解法課件整體誤差ek=y(xk)-yk，下面對(duì)其加以分析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，下面對(duì)其加以分析y1=y0+計(jì)算方法4_常微分方程數(shù)值解法課件歐拉法（續(xù)）思想:用中心差商近似代替微商.注：計(jì)算時(shí)，先用歐拉法求出y1，以后再用二步歐拉法計(jì)算。二步歐拉法(4.4)歐拉法（續(xù)）思想:用中心差商近似代替微商.注：計(jì)算時(shí)，先歐拉法（續(xù)）公式單步否顯式否單步顯式單步隱式二步顯式截?cái)嗾`差y(xn+1)-yn+1

歐拉法（續(xù)）公式單步否顯式否單步顯式單步隱式二步顯式截?cái)嗾`差截?cái)嗾`差Def4.1設(shè)y(xn)

是(4.1)式的精確解，yn是(4.2)式歐拉法得到的近似解，稱y(xn)-yn為歐拉法的整體截?cái)嗾`差.

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

,即第n步計(jì)算是精確的前提下,稱Rn+1=y(xn+1)-yn+1為歐拉法的局部截?cái)嗾`差.

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

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

y(xn+1)-yn+1=

O(h2)由y’(x)

=f(x,y(x))

定理4.4

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

定理4.6

隱式歐拉法的精度是一階,二步歐拉法的精度是二階。證明：對(duì)二步歐拉法進(jìn)行證明，考慮其局部截?cái)嗾`差，令yn=y(xn),yn－1=y(xn－1),將上兩式左右兩端同時(shí)相減：∴二步歐拉法的局部截?cái)嗾`差為O(h3)，其精度是二階。二步歐拉法的局部截?cái)嗾`差Def4.5假設(shè)yn=y(xn)數(shù)值積分法對(duì)右端的定積分用數(shù)值積分公式求近似值：（1）用左矩形數(shù)值積分公式：（2）用梯形公式：——梯形公式梯形公式:將顯示歐拉公式,隱式歐拉公式平均可得梯形公式是隱式、單步公式,其精度為二階數(shù)值積分法對(duì)右端的定積分用數(shù)值積分公式求近似值：（1）用左矩證：令yn=y(xn),由Talor公式有分析：梯形公式是隱式公式，證明其局部截?cái)嗾`差為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)∴梯形公式的截?cái)嗾`差為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個(gè)初值只能由其他方法給出,可能對(duì)后面的遞推精度有影響梯形公式法精度有所提高,但為隱式,需要迭代求解,計(jì)算量大歐拉公式的比較歐拉法最簡單,精度低隱式歐拉法穩(wěn)定性好二步歐拉4.2

改進(jìn)的Euler法

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

預(yù)報(bào)

校正——改進(jìn)的Euler法寫成嵌套公式：平均化形式：4.2改進(jìn)的Euler法Euler公式計(jì)算量小[例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

計(jì)算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法的局部截?cái)嗾`差

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

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

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

精度：二階

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

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

Euler公式用一個(gè)點(diǎn)的值k1作為k*的近似值，而改進(jìn)的Euler公式用二個(gè)點(diǎn)的值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ù)報(bào)幾個(gè)點(diǎn)的ki值并用其加權(quán)平均值作為k*近似值。而構(gòu)造出具有更高精度的計(jì)算公式。多個(gè)斜率加權(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ù)，使(*)式的精度為二階，即使其局部截?cái)嗾`差為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)要使上述方程組有四個(gè)未知量，只有三個(gè)方程，有無窮多組解。取任意一組解便得一種二階龍－庫公式。

當(dāng)c1=c2=1/2,a2=b21=1時(shí)二階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)此為中點(diǎn)法或變形的

Euler公式上述方程組有四個(gè)未知量，只有三個(gè)方程，有無窮多組解。當(dāng)c1三階龍格－庫塔法是用三個(gè)值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)要使其具有三階精度，必須使局部截?cái)嗾`差為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)

三階龍格－庫塔法是用三個(gè)值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)局部截?cái)嗾`差O(h5)

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

m>4時(shí)：計(jì)算量太大，精確度不一定提高，有時(shí)會(huì)降低。類似可推出四階龍格-庫塔公式，常用的有yn+1=yn+(k1Gill公式節(jié)省存儲(chǔ)單元控制舍入誤差Gill公式節(jié)省存儲(chǔ)單元對(duì)于經(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:停止對(duì)于經(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階龍格－庫塔法計(jì)算，局部截?cái)嗾`差為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分兩步計(jì)算y(xn+1)的近似值yn+1(h/2)。則y(xn+1)-yn+1(h/2)≈2c(h/2)p+1

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

注：10

誤差的事后估計(jì)法

20

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

解②：⑴h取大，局部截?cái)嗾`差chp+1大，不精確⑵h取小，運(yùn)算量大（步多），舍入誤差積累大解決策略：變步長龍格－庫塔法

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

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

求yn

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

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

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

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

yn->y(xn)，則稱該法是收斂的。即(xn=x0+nh為固定值)4.5收斂與穩(wěn)定性對(duì)于常微分方程初值問題(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法收斂同樣可驗(yàn)證，若f(x,y)關(guān)于y滿足李普希茲條件,且y0=y(x0)時(shí)，Euler法，標(biāo)準(zhǔn)四階龍格－庫法也收斂。改進(jìn)Euler法的收斂性判定改進(jìn)Euler法的收斂性：若：4.5.2

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

初始值y0有誤差δ0＝