數(shù)值分析：第九章 常微分方程數(shù)值解

1第九章

常微分方程數(shù)值解2常微分方程初始問題

許多科學(xué)技術(shù)問題，例如天文學(xué)中的星體運動，空間技術(shù)中的物體飛行，自動控制中的系統(tǒng)分析，力學(xué)中的振動，工程問題中的電路分析等，都可歸結(jié)為常微分方程的初值問題。

常微分方程中往往只有少數(shù)較簡單和典型的常微分方程（線性常系數(shù)常微分方程等）可求出其解析解。本章考慮的是一階微分方程（組）的數(shù)值解。高階的常微分方程可以轉(zhuǎn)化成一階微分方程組。3舉例例：單擺運動問題OAm可以轉(zhuǎn)化成一階微分方程組4常微分方程解的存在唯一性

定理

考慮一階常微分方程的初值問題對任意定義在上的都成立，則上述方程存在唯一解。只要在上連續(xù),且關(guān)于y滿足Lipschitz條件，即存在與無關(guān)的常數(shù)L使5常微分方程的數(shù)值解

數(shù)值解的基本思想滿足微分方程

常微分方程等價于下列積分方程

常微分方程初值問題的數(shù)值解法，是要尋求解函數(shù)y(x)在一系列節(jié)點

處的近似值

關(guān)鍵在于導(dǎo)數(shù)如何近似。

關(guān)鍵在于積分如何近似。6初始問題的歐拉方法

初值問題的數(shù)值解法一般按節(jié)點從左至右的順序依次求出的近似，所以稱為步進法

單步法

從初值開始,依次求出,后一步的值只依靠前一步的來計算典型的單步法是Euler(歐拉)方法,其計算格式是:

多步法在計算時需要前面若干步值來計算。7歐拉方法的推導(dǎo)

以Euler法為例說明初值問題數(shù)值方法的三種基本途徑

Taylor展開法

忽略高階項,取近似值可得到Euler公式

數(shù)值微分法,用差商代替導(dǎo)數(shù)

8歐拉方法的推導(dǎo)則

數(shù)值積分法

將區(qū)間積分右端采用左矩形公式數(shù)值積分，得9幾個簡單的數(shù)值方法

顯式歐拉法

隱式歐拉法方法上是向后差商近似導(dǎo)數(shù)，或右矩形數(shù)值積分得到

由于未知數(shù)yn+1同時出現(xiàn)在等式的兩邊，不能直接得到，故稱為隱式/*implicit*/歐拉公式.

一般先用顯式計算一個初值，再迭代求解。10幾個簡單的數(shù)值方法

梯形公式

顯、隱式兩種算法的平均或梯形公式數(shù)值積分得到

改進的歐拉公式

先用顯式歐拉公式作預(yù)測，算出),(nnnyxfhy+=1+ny

再將代入隱式梯形公式的右邊作校正，得到1+ny)],(),([211++++=nnnnnxfyxfhyy1+ny11單步法的誤差分析

整體截斷誤差

數(shù)值解和精確解之差

整體截斷誤差除與步計算有關(guān)外,還與的計算有關(guān)，所以一般很難分析清楚。

局部截斷誤差一般單步法有增量形式

在假設(shè)yn=y(xn)，即第n步計算是精確的前提下，考慮截斷誤差Tn+1=y(xn+1)

yn+1稱為局部截斷誤差

12算法的精度若某算法的局部截斷誤差為O(hp+1),則稱該算法有p階精度。顯式歐拉，隱式歐拉為1階精度

梯形公式和改進歐拉公式為2階精度

顯式歐拉局部截斷誤差梯形公式局部截斷誤差13例例：考察初值問題在區(qū)間[0,0.5]上的解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解。0.00.10.20.30.40.5精確解改進歐拉法

歐拉隱式歐拉顯式

節(jié)點xi1.0000

2.00004.0000

8.00001.6000101

3.2000101

1.00002.5000101

6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107從表中可以看出，這幾種方法計算誤差比較大，說明這些方法精度不夠，需尋找更高精度的方法14高精度方法

基于Taylor展開式可以建立任意精度的單步法

而所以，可以構(gòu)造格式

這種格式使用到了各階偏導(dǎo)數(shù)，使用不便。15Runge-Kutta方法

Runge-Kutta方法的基本思想

只要對平均斜率提供一種算法，便相應(yīng)地得到一種微分方程的數(shù)值計算公式。16Runge-Kutta方法

考察改進的歐拉法可以將其改寫為：

改進歐拉公式卻是利用了xn與xn-1兩個點的斜率值k1=f(xn,yn)與k2=f(xn+1,yn+hk1)取算術(shù)平均作為平均斜率的近似值。17Runge-Kutta方法啟示設(shè)法在[xk,xk+1]上多預(yù)報幾個點的斜率值，然后將它們加權(quán)平均作為平均斜率的近似值，則有可能構(gòu)造出更高精度的計算公式推廣其中

由f(x,y)在一些點的值加權(quán)平均所得這種方法稱為Runge-Kutta方法,簡記為R-K公式.若

由r個f(x,y)的值加權(quán)平均所得的公式稱為r級R-K公式

18二級Runge-Kutta方法

二級Runge-Kutta方法),(),(][12122111phKyphxfKyxfKKKhyyiiiiii++==++=+ll

希望能確定系數(shù)

1、

2、p，使得到的算法格式有2階精度，即在的前提假設(shè)下，使得

將K2在(xi,yi)點作Taylor展開19二級Runge-Kutta方法

將K2代入第1式，得到

將yi+1與y(xi+1)在xi點的泰勒展開作比較20二級Runge-Kutta方法要求，則必須有：

這里有3個未知數(shù)，2個方程,解不唯一

就是改進的歐拉法。

就是中點公式。

21三階與四階顯式R-K方法

三級Runge-Kutta方法其中及均為待定參數(shù).

局部截斷誤差為

只要將按二元函數(shù)泰勒展開，使，22三階與四階顯式R-K方法可得待定參數(shù)滿足方程

這是8個未知數(shù)6個方程的方程組，解也不是惟一的.23三階與四階顯式R-K方法一個特殊的方法

同樣的方法，可以得到經(jīng)典4階R-K方法24注更一般R-K方法

龍格-庫塔法的主要運算在于計算的值，即計算f的值。Butcher于1965年給出了計算量與可達到的最高精度階數(shù)的關(guān)系：25注753可達到的最高精度642每步須算Ki的個數(shù)

由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用低階算法而將步長h

取小。26

例

解設(shè)取步長，從到用四階龍格-庫塔方法求解初值問題這里，經(jīng)典的四階龍格-庫塔公式為2728F='y-2*x/y';a=0;b=1;h=0.1;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=1;fori=1:nx=X(i);y=Y(i);

K1=h*eval(F);

x=x+h/2;y=y+K1/2;K2=h*eval(F);x=x;y=Y(i)+K2/2;K3=h*eval(F);x=X(i)+h;y=Y(i)+K3;K4=h*eval(F);Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;end2930

求微分方程（組）的解析解命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始條件’,‘自變量’)

求微分方程（組）的數(shù)值解命令

[t，x]=solver（’f’,ts,x0,options）ode45ode23ode113ode15sode23s由待解方程寫成的m-文件名ts=[t0，tf]，t0、tf為自變量的初值和終值函數(shù)的初值函數(shù)值用于設(shè)定誤差限(缺省時設(shè)定相對誤差10-3,絕對誤差10-6),命令為：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分別為設(shè)定的相對誤差和絕對誤差.31單步法的收斂性定義若某算法對于任意固定的x=x0+nh，當(dāng)h0

(即n

)時有yn

y(xn)，則稱該算法是收斂的。

例

就初值問題考察歐拉顯式格式的收斂性。解：該問題的精確解為歐拉公式為對任意固定的x=xn=nh

，有

32定理:對初值問題的單步法,若局部截斷誤差為,且函數(shù)

對y滿足Lipschitz條件,即存在L>0,使得對一切成立,則該方法收斂,且有單步法的收斂性收斂性的條件歐拉，改進歐拉以及Runge-Kutta方法都是收斂的

整體截斷誤差比局部截斷誤差低一階33

對于歐拉方法，由于其增量函數(shù)就是，故當(dāng)關(guān)于滿足利普希茨條件時它是收斂的.

改進的歐拉方法，其增量函數(shù)由歐拉方法的收斂性

改進歐拉方法的收斂性

假設(shè)關(guān)于滿足利普希茨條件，記利普希茨常數(shù)為,34單步法的穩(wěn)定性

考慮收斂性時，當(dāng)h0時有yn

y(xn)，此時yn為理論數(shù)值解，也就沒有考慮到計算過程中的舍入誤差

當(dāng)考慮計算過程中舍入誤差對最后結(jié)果是否有影響，所以引入數(shù)值算法絕對穩(wěn)定的概念定義

若一種數(shù)值方法在節(jié)點值上大小為擾動，于以后各節(jié)點值上產(chǎn)生的偏差均不超過，則稱該方法是絕對穩(wěn)定的.由于f(x,y)比較復(fù)雜，給穩(wěn)定性研究帶來困難，所以考慮試驗方程35單步法的穩(wěn)定性

試驗方程

該試驗方程為線性方程，其精確解是穩(wěn)定的。若對該方程數(shù)值算法都穩(wěn)定，則對其它方程也靠不住。對一般非線性方程，局部可以化成線性的來近似

單步法的穩(wěn)定性

單步法用于試驗方程，若得到的解為滿足

則是絕對穩(wěn)定的。36單步法的穩(wěn)定性記，，的范圍叫絕對穩(wěn)定區(qū)域

hlh=h絕對穩(wěn)定區(qū)域與實軸之交叫絕對穩(wěn)定區(qū)間

方法A比方法B的絕對穩(wěn)定區(qū)域大稱A比B穩(wěn)定顯式歐拉法的穩(wěn)定性0-1-2ReImgh在復(fù)平面上是以(-1,0)為圓心的單位圓絕對穩(wěn)定區(qū)間為【-2,0】37隱式歐拉法的穩(wěn)定性210ReImg所以絕對穩(wěn)定區(qū)域為左半平面絕對穩(wěn)定區(qū)間隱式歐拉比顯式歐拉穩(wěn)定性好38梯形公式的穩(wěn)定性所以絕對穩(wěn)定區(qū)域為左半平面絕對穩(wěn)定區(qū)間39經(jīng)典4階龍格-庫塔方法的穩(wěn)定性40

研究一般方程的穩(wěn)定性時，相當(dāng)于模型方程中的λ

41線性多步法

一個特例對于初值問題若用中心差商離散導(dǎo)數(shù)稱為二步法42線性多步法

推廣當(dāng)

10時，為隱式公式;

1=0則為顯式公式。

43線性多步方法的構(gòu)造

基于Taylor展開的構(gòu)造方法: