常微分方程組數(shù)值解

常微分方程組數(shù)值解其中x

是質量，m是離開平衡位o的距離，t為時間，c為彈簧系數(shù)。

例如：彈簧一質量系統(tǒng)的振動問題經一定的簡化后可用一個二階常微分方程來描述。

在具體求解微分方程時，需具備某種定解條件，微分方程和定解條件合在一起組成定解問題。定解條件有兩種：一種是給出積分曲線在初始點的狀態(tài)，稱為初始條件，相應的定解問題稱為初值問題。另一類是給出積分曲線首尾兩端的狀態(tài)，稱為，相應的定解問題稱為邊值問題。mxxoc

我們現(xiàn)在討論常微分方程的數(shù)值解法。先從最簡單的一階常微分方程的初值問題出發(fā)開始討論。

由常微分方程理論可知：只要上式中的函數(shù)f(x,y)在區(qū)域G={a≤x≤b,－∞＜y＜∞}內連續(xù)，且關于y滿足Lipschitz條件，即存在與x,y無關的常數(shù)L，使

至于初值問題（1）的數(shù)值解法，常采用差分方法，即把一個連續(xù)的初值問題離散化為一個差分方程來求解。具體地，將（1）離散化后，求找其解y=y(x)在一系列離散點

下面分析均假定滿足上述條件。

對于初值問題（1），先將其離散化，即把[a,b]區(qū)間n等分，得各離散節(jié)點一、Euler公式§2Euler方法

因為初值問題中的初始條件已知，即可利用已知的來求出下一節(jié)點處的近似值；再從來求如此繼續(xù)，直到求出為止。這種用按節(jié)點的排列順序一步一步地向前推進的方式求解的差分算法稱為“步進式”或“遞推式”算法，它是初值問題數(shù)值解法的各種差分格式的共同特點。因此，只要能寫出由前幾步已知信息來計算的遞推公式（即差分格式），即可完全表達這種算法。若將和的近似值分別記為和，則得：（3）

這就是Euler公式（格式）。利用它可由初值出發(fā)逐步算出。這類形式的方法也稱為差分方法。定義：如果局部截斷誤差為，則這種數(shù)值算法的精度為p階，故Euler格式的精度為一階。從幾何意義上來看，如圖，當假定為準確值，即在

的前提下來估計誤差，這種截斷誤差稱為局部截斷誤差。由（2）、（3）知Euler公式在處的局部截斷誤差為：由方程（1）知，其積分曲線y=f(x)上任一點(x,y)的切線斜率都等于函數(shù)f(x,y)的值。從初始點（即點）出發(fā)，作積分曲線y=y(x)在點上的切線（其斜率為）與直線相交與點（即點），得到作為的近似值，則有yOxy=f(x)相比較知，這時用切線近似代替了曲線段點近似代替了點，近似代替了近似代替了。遞推繼續(xù)從點出發(fā)，作一斜率為的直線與直線相交于點（即點），得到作為的近似值?！绱酥钡近c。這樣得出一條折線近似代替積分曲線，當步數(shù)越多時，由于誤差的積累，折線可能會越

解：為便于進行比較，我們后面將用多種數(shù)值方法求解上述初值問題。這里先用Euler公式，此處具體格式為：取步長為h=0.1，計算結果略。

由結果可見Euler方法的精度很差。即為Euler格式（3）。

因為差商是微分的近似，所以Euler格式也可用差商近似代替導數(shù)的離散方法來得到。在節(jié)點處有：二后退Euler格式

顯然Euler格式具有遞推性，在計算時只要用到前一步所得的結果一個信息就夠了，因此是一種單步格式或稱一步格式。若用不同的數(shù)值微分計算方法也可導出其它形式的算法。例如：用向后差商表示的數(shù)值微分公式

（6）稱為向后Euler公式，又稱為隱式Euler公式(后退Euler格式)。后退Euler公式與Euler公式有著本質的區(qū)別，后者是關于的一個直接計算公式，這類公式稱作顯式的；而前者，即（6）中右端含有未知的，它實際上是關于的一個函數(shù)方程。這類公式稱作隱式的。顯式與隱式兩類方法各有特點，使用顯式算法遠比隱式算法方便，但考慮數(shù)值穩(wěn)定性等因素，人們常選用隱式算法。

隱試算法（6）常用迭代法來實現(xiàn)，而迭代過程實質上是逐步顯式化。

設用顯式Euler格式算出作迭代初值，以此代入（6）右端，使之轉化為顯示，直接計算得：，再用代入（6）右端又有：若迭代過程收斂，則極限值必為隱式方程（6）的解，從而可獲得后退Euler方法的解。如此反復進行可得序列

。從幾何上看，梯形公式是取區(qū)間兩端點處斜率的平均斜率。比較Euler公式和后退歐拉公式的截斷誤差公式（4），（7）易見兩者為同階，但符號相反，因此就想到將這兩種算法進行算術平均，其結果可能消除誤差的主要部分，而獲得更高精度的結果。這種平均方法稱為梯形方法，公式為：三梯形公式Euler方法是過點以斜率引直線交的點A。后退Euler方法是以點處的斜率為斜率，從點引直線交于另一點B。

A、B兩點都偏離點點，梯形方法就是取A、B兩點的中點P作為Q2近似，上圖表明梯形方法確實改善了精度。y=f(x)xyoABP從而也可得二步Euler公式及其截斷誤差為：也可以由導數(shù)的中心差分近似式得到：四二步歐拉公式將區(qū)間[a,b]n等分，得子區(qū)間在第i+1個子區(qū)間上，積分

Euler公式、后退Euler公式以及梯形公式，二步Euler公式均可利用前面的數(shù)值積分技術得到。例如：對于初值問題：對右端利用左矩形公式可得即Euler格式各公式的截斷誤差可直接利用數(shù)值積分截斷誤差估計而得。從而可知梯形公式（8）的截斷誤差為：梯形公式也是隱式的，可用迭代法求解，與后退Euler方法一樣，仍用Euler方法提供迭代初值，其迭代格式為：為分析迭代過程的收斂性，將（12）與（8）相減得：（12）k=0，1，2，…..L為f(x,y)關于y的Lipschitz常數(shù).如果選取h充分小使得則.時有這表明迭代過程（12）是收斂于（8）的解的。

五改進的Euler公式這一公式也可改寫為：

（15）(14)預測（13）

校正

上面已看到Euler公式計算量小但精度差，梯形方法雖然提高了計算精度，但算法復雜計算量大，在應用（12）進行迭代時，每次均要計算函數(shù)f的值，而迭代又要反復進行多次，計算量很大，難以預測。為了控制計算量，通常只迭代一兩次就轉入下一步計算，這就簡化了計算。

具體地，先用Euler公式求得一個初步的近似值稱之為預測值。預測值的精度可能很差。再用梯形公式（8）將它校正一次，即按（12）式迭代一次得，這個結果稱為校正值，這樣建立的預測校正系統(tǒng)稱為改進的Euler公式。

式：上式表為下列平均化形（16）

歐拉公式和改進歐拉公式分及兩種格式的計算結果分別列表如下別為：解：取步長h=0.1并比較兩法所得計算結果的精度。例：試分別用函數(shù)公式和改進的歐拉公式求解：

由表可見，與精確解相比，改進的Euler公式的精度較Euler公式有明顯的提高。下面再看兩步Euler公式（9），除了給出初值外，還需要借助其它單步法（如Euler公式，后退Euler公式及梯形公式等）再提供一個Euler公式改進的Euler公式精確解01110.10.20.310.90000000.81000000.72900000.34867840.90500000.81902500.74121760.36854100.90483740.81873080.74081820.3678794啟動值然后才能啟動計算公式依次計算

用兩步Euler公式與梯形公式相匹配，又可得到下面預測-校正系統(tǒng)：校正：（17）（18）兩步法優(yōu)美是由于它調用了兩個節(jié)點上的信息，從而能以較少的計算量獲得較高的精度。

預測：與改進的Euler公式（13）（14）相比較易見（17）（18）的一個突出特點是它的預測公式與校正公式具有同階精度。據此可以比較方便的估計截斷誤差，并基于這種估計，可以提供一種提高精度的簡易方法。再由梯形公式截斷誤差公式（11）知：校正公式（18）的截斷誤差為：

比較（19）（20）可見，校正值的誤差大約只是預測值的誤差的1/4(符號相反).即

,由此導出誤差的事后估計:

若預測公式（17）中的和都是準確的。即則由兩步Euler公式的截斷誤差公式（10）知：1預測:

2改進:

3計算:

可以期望利用這樣估計出的誤差作為計算結果的一種補償,有可能使精度得到改善.

以和分別表示第I步的預測值和校正值,按估計式(21),

和分別作為和的改進值,在校正值尚未求出之前,可用上一步的偏差替代來改進預測值.這樣設計的計算方案有如下六個環(huán)節(jié):4校正:

5改進:6計算：§3Runge-Kutta方法有解y=y（x）。Taylor展開有：(*)一Taylor級數(shù)法運用上述方案計算時，要用到先一步的信息，，和更前一步的。因此啟動算法之步必須給出啟動值和?？捎闷渌鼏尾椒ǎɡ绺倪M的Euler方法）來計算，則一般取為0。計算結果表明，這種簡單的處理方法通?？梢垣@得令人滿意的效果。而具體則有:

式中y（x）的各階導數(shù)可由方程（*）用函數(shù)f來表達。引進函數(shù)序列來描述求導過程：若在(1)中右端取前面若干項,并且在處按(2)計算系數(shù)的近似值的近似值,可導出如下Taylor公式:

其中p=1時一階Taylor展開式即為

而提高Taylor公式的階p即可提高計算結果精度.P階Taylor展開的局部截斷誤差為:Euler公式因而具有P階精度:

對于一階常微分方程（1）的解y=y(x)，可利用中值定理得：，即也即式中K=（5）二Runge-Kutta公式的導出例:用Taylor公式求解初值問題:K可看作是y=y(x)在區(qū)間上的平均斜率。從而Euler公式相當于取點上的斜率作為平均斜率K的近似值，這當然十分粗糙，因而精度必然很低。再考慮改進Euler公式（15）可改寫成：

和兩個點上的斜率和的算術平均值作為（4）中平均斜率K的近似值。其中是通過已知與（4）比較可見，它相當于把信息來近似地預測的。這個過程啟示我們，如果設法在內多預測幾個點的斜率值，然后將它們加權平均作平均斜率K的近似值，就有可能構造出更高精度的計算公式。這就是Runge-Kutta方法的基本思想。的現(xiàn)在，設想取區(qū)間內某一節(jié)點上斜率與點上的斜率作線性組合（即加權平均）化為平均斜率K的近似，即（4）化為：為了要得到點上的斜率，需先預測根據預測值再來算出由此構造出計算格式：（6）上式中含三個待定系數(shù)和p，適當選定它們以使算法的局部誤差為，從而具有二階精度。假定，分別將和作Taylor展開得：由（2）得：將此式與y(x)在處的二階Taylor展開在處的取值相比較：代入（6）得：成立，則（6）的局部截斷誤差就等于，從而能具有二階精度。由系數(shù)比較知，只要：

（7）中三個待定參數(shù)P，但只有兩個方程，因此還有一個自由度。凡滿足條件（7）的一族格式統(tǒng)稱為二階Runge-Kutta格式。

當p=1，時，二階Runge-Kutta格式（6）即為改進的Euler格式（15）。如取p=1/2，則，二階R-K格式（6）成為：稱之為變形的Euler格式。

由于（8）中的是由Euler格式預測出來的區(qū)間中的點的近似解，就近似地等于此中點的斜率，因此（8）就相當于用中點的斜率作為（4）中平均斜率K的近似值，故格式（8）也稱為中點格式。

總之，二階R-K格式用多算一次函數(shù)值f的辦法，避開了二階Taylor級數(shù)法所要求計算的f的導數(shù)值，在這種意義上可以說，R-K方法實質上是Taylor方法的變形。并用三個點的斜率值線性組合得到平均斜率K，這時計算公式為：其中仍用公式（6）所取的形式。

為了預測點的斜率值，在區(qū)間內有兩個為了進一步提高精度，設除了外，再考慮一點

三三階Runge-Kutta方法

表面上，只含一個斜率值，但要通過才能算出來，因此式中隱含著，每完成一步仍然需計算f的兩次函數(shù)值，其工作量仍與改進的Euler格式一樣。斜率值和可以利用。我們用和線性組合給出區(qū)間上的平均斜率，從而得到的預測值于是，再通過計算函數(shù)值f得到：這樣設計出的計算格式具有形式：（9）我們希望適當選擇系數(shù)和p、q、r、使上述公式具有三階精度。為便于數(shù)學演算，引進算子：，則根據（2）有：，，于是三階Ｔａｙｌｏｒ展開可表為：式中的下標i均表示在處取值。將的展開式帶入（9），再令可得：此式與（10）式比較系數(shù)可得：以及稱滿足（11）.的公式族（9）為三階Runge-Kutta公式（11）（12）四四階Rung-Kutta方法最常用的一種經典Rung-Kutta格式具體形式如下：(13)四階Rung-Kutta方法每一步需要四次計算函數(shù)值f,可以證明其截斷誤差為不過證明復雜。例：試分別用Euler方法（h=0.025）改進的Euler（h=0.05）及經典R-K方法（h=0.1）求解下列初值問題，比較三種方法的結果的精度。解：三種方法如下：Euler格式：改進的Euler格式：(h=0.05)經典的R-K格式：(h=0.1)Euler法h=0.025改進Eulerh=0.05R-K法h=0.1精確解yx=0.1x=0.2x=0.30.9036878900.8166518030.7379983450.9048765620.8188015930.740914370.9048375000.8187309010.7408180.904837480.818730750.74081822x這里采用了不同的步長h值，是為了使他們所耗的計算工作量大致相同，以便于比較。由上表可見，經典的R-K方法的精確度較改進的Euler方法又有很大的提高。這一結論也可以從理論上大致的分析出來。析出來：Euler方法的局部截斷誤差為：計算四步后的而經典R-K方法的局部截斷誤差則為：可見，當為大致相同數(shù)量級的常數(shù)時有：但要注意的是：R-K方法的導出利用了Taylor展開，因此要求所求的解有教好的光滑性，如果解的光滑性差，則采用經典的R-K方法所的數(shù)值解，其精度有可能反而不及改進的Euler方法，因此在實際計算中應根據問題的具體情況來選擇適合的算法。

五步長的自動選擇---變步長的Runge-Kutta方法

在應用數(shù)值法求解微分方程中，選擇適當?shù)牟介L是至關重要的。步長太大則達不到要求，步長太小則步數(shù)增多，不但增加計算工作量，還可能導致舍入誤差的嚴重積累。尤其是當微分方程的解y(x)變化激烈時，步長的合理取法是在變化激烈處步長取小些，在變化平緩時取大些，也就是采取自動變步長的方法，即根據精度的要求先估計出下一步長的合理大小，然后按此計算。

從節(jié)點出發(fā)，先以h為步長跨一步到節(jié)點

值如果計算公式是P階的，則當h的變化不大時，式中的系數(shù)C可以近似的看作為常數(shù)。然后將步長減半，即以為步長，從節(jié)點出發(fā)，跨兩步到節(jié)點，再求得一個近似值。其中每跨一步的截斷誤差為：故有：(14)（15），求一個近似

下面介紹一種Richrdson外推法：作為近似值，則

的精確度都要高。當p=4時，可以取這種修正方法與Romberg的數(shù)值積分的思路是一樣的。（15）除以（14）得：由此可以得出誤差事后估計式：

由上面的分析可見，微分方程數(shù)值解的基本思想是，通過（1）如果，則反復加倍步長進行計算，直到時為止，并以上一次步長的計算結果作為（2）若，則反復減半步長進行計算，直到時為止，并取其最后一次步長的計算作為這樣做時，為了選擇步長，每一步都要反復判別，增加了工作量，但在方程的解y(x)變化劇烈的情況時，總的計算工作量可以得到減少，結果還是合算的。這樣就可以從步長減半前后的兩次計算結果的偏差來判斷步長選的是否適當，當要求的數(shù)值精度為時：§4單步法的收斂性和穩(wěn)定性一單步法的收斂性本例的Euler公式為由此式遞推可得：定義：若一種數(shù)值方法對任意固定的當（同時）時有，則稱方法是收斂的。某種離散化手段，將微分方程轉化為差分方程來求解。這種轉化是否合理，還要看差分方程問題的解當時是否會收斂到點對固定的i將趨向，這時討論收斂是沒有意義的。考察Euler方法的收斂性

例如：對初值問題所謂的單步法，就是在計算時只用到它前一步的信息。Taylor級數(shù)法，Runger-Kutta方法等都是單步法的例子。顯然單步法的共同特征是，他們都是將加上某種形式的增量得出，其計算公式形如：式中的稱作增量函數(shù)。例如：對Euler公式，有對改進的Euler公式有：關于單步法有下述的收斂定理：(1)(2)定理：設單步法（1）具有p階的精度，且增量函數(shù)關于y的滿足Lipschitz條件又設初值是準確的，即y0=y(x0),則其的整體截斷誤差為:(3)證明:設以表示取用公式（1）求的結果，即

則由于所給的方法具有P級精度，根據精度的定義可知，有常數(shù)C使由（4）與（1）得：(4)

此定理表明，判斷單步法（1）的收斂性，歸結為驗證增量函數(shù)能否滿足Lipschitz（3）對于Euler方法，由于其增量函數(shù)故當f滿足Lipschitz條件時它是收斂的。因此改進的Euler方法也是收斂的。類似的可以證明其它的Runger-Kutta方法的收斂性。上面關于收斂性的討論有前提條件，即必須假定數(shù)值方法本身的計算是準確的，實際情況并非這樣，差分方程的求解還會有計算的誤差。如由于數(shù)字的舍入的誤差而引起的小擾動。這類小擾動在傳播的過程中會不會惡性的增長？以至于“淹沒”了差分方程的“真解”呢？這就是差分方程的穩(wěn)定問題。在實際計算中，某一步產生的擾動值在后面的計算中能夠被控制，甚至是逐步衰減的。

穩(wěn)定性問題比較復雜，為簡化討論，僅考察下面的模型方程為保證微分方程本身的穩(wěn)定性，假設先考慮Euler方法的穩(wěn)定性。模型方程的Euler公式為：設在節(jié)點值上有一擾動值,它的傳播使節(jié)點值產生大小為：的擾動值，假設用按Euler公式得出：定義：若一種數(shù)值方法在節(jié)點值上有擾動，而對于yi后的各個節(jié)點值上產生的偏差均不超過,則稱該方法是穩(wěn)定的。二單步法的穩(wěn)定性0.0250.0500.0750.100123450-1-2-3xy其方程時間常數(shù)為因此有（10）知，要使Euler法穩(wěn)定則步長如果取步長h=0.025則Euler格式為：其結果在準確解上下波動不穩(wěn)定，再看后退Euler格式H=0.025時，形式為：計算結果穩(wěn)定。具體結果如下：節(jié)點Eule方法后退Euler方法0.025-1.50.28570.050.2.250.08160.075-3.3750.02330.1005.06250.0067

前面介紹的幾種步進方法，在計算時大多只用到前一個節(jié)點上的近似值，而沒有用到前幾步的計算所得出的信息，故稱為單步法。實際上經過多次單步法計算后，已得出一系列近似值等。為了充分利用這些信息來計算以減少計算工作量和獲得教高的精度，可采用如下計算公式：§5線性多步法一顯式Adams法（1）j對應不同k的和可計算出來分別列表如下：k例如k=3時有：