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

在工程和科學技術的實際問題中，常需求解微分方程，但常微分方程中往往只有少數(shù)較簡單和典型的常微分方程（例如線性常系數(shù)常微分方程等）可求出其解析解,對于變系數(shù)常微分方程的解析求解就比較困難，而一般的非線性常微分方程的求解困難就更不用說了。大多數(shù)情況下，常微分方程只能用近似方法求解。這種近似解法可分為兩大類：一類是近似解析法，如級數(shù)解法、逐次逼近法等;另一類是數(shù)值解法，它給出方程在一些離散點上的近,第八章 常微分方程數(shù)值解法,1、 引 言,似值。,其中 x 是質量，m是離開平衡位o的距離，t為時間，c為彈簧系數(shù)。,在具體求解微分方程時，需具備某種定解條件，微分方程和定解條件合在一起組成定解問題。定解條件有兩種：一種是給出積分曲線在初始點的狀態(tài)，稱為初始條件，相應的定解問題稱為初值問題。另一類是給出積分曲線首尾兩端的狀態(tài)，稱為邊界條件，相應的定解問題稱為邊值問題。,我們現(xiàn)在討論常微分方程的數(shù)值解法。先從最簡單的一階常微分方程的初值問題出發(fā)開始討論。,由常微分方程理論可知：只要上式中的函數(shù)f(x,y)在區(qū)域 G=axb,y內連續(xù)，且關于 y 滿足Lipschitz條件，即存在與 x, y 無關的常數(shù)L，使,至于初值問題（1）的數(shù)值解法，常采用差分方法，即把一個連續(xù)的初值問題離散化為一個差分方程來求解。具體地，將（1）離散化后，求找其解 y=y(x)在一系列離散點,下面分析均假定滿足上述條件。,一、Euler公式,2 Euler方法,因為初值問題中的初始條件 已知，即可利用已知的 來求出下一節(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) 在 點上的切線 （其斜率為 ）與直線 相交與 點（即 點），得到 作為 的近似值，則有,相比較知，這時用切線 近似代替了曲線段 點近似 代替了 點， 近似代替了 近似代替了 。 遞推繼續(xù)從 點出發(fā)，作一斜率為 的直線 與直線 相交于 點（即 點），得到 作為 的近似值 。如此直到 點。這樣得出一條折線 近似代替積分曲線 ，當 步數(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）右端又有：,從幾何上看，梯形公式是取 區(qū)間兩端點處斜率的 平均斜率。,三 梯形公式,Euler方法是過 點以斜率 引直線交 的點A。 后退Euler方法是以點 處的斜率 為斜率，從 點引直線交 于另一點B。 A、B兩點都偏離點 點，梯形方法就是取A、B兩點的中點P作為Q2近似，上圖表明梯形方法確實改善了精度。,從而也可得二步Euler公式及其 截斷誤差為：,也可以由導數(shù)的中心差分近似式得到：,四 二步歐拉公式,將區(qū)間a,bn等分，得子區(qū)間 在第i+1 個子區(qū)間上，積分,例如：對于初值問題：,對右端利用左矩形公式可得 即 Euler格式,各公式的截斷誤差可直接利用數(shù)值積分截斷誤差估計而得。從而可知梯形公式（8）的截斷誤差為：,梯形公式也是隱式的，可用迭代法求解，與后退Euler方法一樣，仍用Euler方法提供迭代初值，其迭代格式為：,為分析迭代過程的收斂性，將（12）與（8）相減得：,（12）,k=0，1，2，,L為f(x,y)關于y的Lipschitz常數(shù).如果選取h充分小使得 則 . 時有 這表明迭代過程（12）是收斂于（8）的解的。,五 改進的Euler公式,上面已看到Euler公式計算量小但精度差，梯形方法雖然提高了 計算精度，但算法復雜計算量大，在應用（12）進行迭代時，每次 均要計算函數(shù)f的值，而迭代又要反復進行多次，計算量很大，難以 預測。為了控制計算量，通常只迭代一兩次就轉入 下一步計算， 這就簡化了計算。,具體地，先用Euler公式求得一個初步的近似值 稱之為預測 值。預測值 的精度可能很差。再用梯形公式（8）將它校正 一次，即按（12）式迭代一次得 ，這個結果稱為校正值，這 樣建立的 預測校正系統(tǒng)稱為改進的Euler公式。,由表可見，與精確解 相比，改進的Euler公式的精度較Euler 公式有明顯的提高。 下面再看兩步Euler公式（9），除了給出初值 外，還需要借助 其它單步法（如Euler公式，后退Euler公式及梯形公式等）再提供一個,Euler公式,改進的Euler公式,精確解,0,1,1,1,0.1,0.2,0.3,1,0.9000000,0.8100000,0.7290000,0.3486784,0.9050000,0.8190250,0.7412176,0.3685410,0.9048374,0.8187308,0.7408182,0.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é):,3 Runge-Kutta方法,運用上述方案計算 時，要用到先一步的信息 ， ， 和 更前一步的 。因此啟動算法之步必須給出啟動值 和 。 可用其它單步法（例如改進的Euler方法）來計算， 則一般取為0。計算結果表明，這種簡單的處理方法通常 可以獲得令人滿意的效果。,而具體則有:,式中y（x）的各階導數(shù)可由方程（*）用函數(shù)f來表達。引進函數(shù) 序列 來描述求導過程：,二 Runge-Kutta公式的導出,例:用Taylor公式求解初值問題:,K可看作是 y=y(x)在區(qū)間 上的平均斜率。從而Euler公式 相當于取 點上的斜率 作為平均斜率K的 近似值，這當然十分粗糙，因而精度必然很低。,這個過程啟示我們，如果設法在 內多預測幾個點的斜率 值，然后將它們加權平均作平均斜率K的近似值，就有可能構造 出更高精度的計算公式。這就是Runge-Kutta方法的基本思想。,根據預測值 再來算出 由此構造出計算格式：,假定 ，分別將 和 作Taylor展開得：,由（2）得：,（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方法的變形。,三 三階Runge-Kutta方法,斜率值 和 可以利用。我們用 和 線性組合給出區(qū)間 上的平均斜率，從而得到 的預測值,于是，再通過計算函數(shù)值f得到：,這樣設計出的計算格式具有形式：,（9）,我們希望適當選擇系數(shù) 和p、q、r、 使上述公式具有三階精 度。為便于數(shù)學演算，引進算子：,，,則根據（2）有：,于是三階展開,可表為：,式中 的下標 i 均表示在 處取值。,（12）,四 四階Rung-Kutta方法,最常用的一種經典Rung-Kutta格式具體形式如下：,(13),解：三種方法如下：,Euler格式：,改進的Euler格式：,(h=0.05),經典的R-K格式：,(h=0.1),Euler法h=0.025,改進Eulerh=0.05,R-K法h=0.1,精確解y,x=0.1,x=0.2,x=0.3,0.903687890,0.816651803,0.737998345,0.904876562,0.818801593,0.74091437,0.904837500,0.818730901,0.740818,0.90483748,0.81873075,0.74081822,x,這里采用了不同的步長h值，是為了使他們所耗的計算工作量 大致相同，以便于比較。由上表可見，經典的R-K方法的精確度 較改進的Euler方法又有很大的提高。這一結論也可以從理論上 大致的分析出來。,析出來：Euler方法的局部截斷誤差為：,計算四步后的,而經典R-K方法的局部截斷誤差則為：,可見，當,為大致相同數(shù)量級的常數(shù)時有：,但要注意的是：R-K方法的導出利用了Taylor展開，因此要求 所求的解有教好的光滑性，如果解的光滑性差，則采用經典的 R-K方法所的數(shù)值解，其精度有可能反而不及改進的Euler方法， 因此在實際計算中應根據問題的具體情況來選擇適合的算法。,五 步長的自動選擇-變步長的Runge-Kutta方法,在應用數(shù)值法求解微分方程中，選擇適當?shù)牟介L是至關重要的。步 長太大則達不到要求，步長太小則步數(shù)增多，不但增加計算工作量， 還可能導致舍入誤差的嚴重積累。尤其是當微分方程的解y(x)變化 激烈時，步長的合理取法是在變化激烈處步長取小些，在變化平緩 時取大些，也就是采取自動變步長的方法，即根據精度的要求先估 計出下一步長的合理大小，然后按此計算。,作為近似值，則,的精確度都要高。,當p=4時，可以取,這種修正方法與Romberg的數(shù)值積分的思路是一樣的。（15） 除以（14）得：,由此可以得出誤差,事后估計式：,由上面的分析可見，微分方程數(shù)值解的基本思想是，通過,（1） 如果 ，則反復加倍步長進行計算，直到 時為止，并以上一次步長的計算結果作為 （2）若 ， 則反復減半步長進行計算，直到 時為止，并取其最后一次步長的計算作為 這樣做時 ，為了選擇步長，每一步都要反復判別， 增加 了工作量，但在方程的解 y(x) 變化劇烈的情況時 ，總的計算 工作量可以得到減少，結果還是合算的。,這樣就可以從步長減半前后的兩次計算結果的偏差 來判斷步長選的是否適當，當要求的 數(shù)值精度為 時：,4 單步法的收斂性和穩(wěn)定性,一 單步法的收斂性,本例的Euler公式為 由此式遞推可得：,定義 ： 若一種數(shù)值方法對任意固定的 當 （同時 ） 時 有 ，則稱方法是收斂的。,某種離散化手段，將微分方程轉化為差分方程來求解。這種轉化是 否合理，還要看差分方程問題的解 當 時是否會收斂到點 對固定的i將趨向 ，這時討論收斂是沒有意義的。,此定理表明，判斷單步法（1）的收斂性，歸結為驗證增量函數(shù) 能否滿足Lipschitz（3） 對于Euler方法，由于其增量函數(shù) 故當f滿足 Lipschitz 條件時它是收斂的。,穩(wěn)定性問題比較復雜，為簡化討論，僅考察下面的模型方程 為保證微分方程本身的穩(wěn)定性，假設 先考慮Euler方法的穩(wěn)定性。模型方程 的Euler公式為： 設在節(jié)點值 上有一擾動值 , 它的傳播使節(jié)點值 產生大小為： 的擾動值，假設用 按Euler 公式得出：,定義：若一種數(shù)值方法在節(jié)點值 上有擾動 ，而對于yi后的各個 節(jié)點值 上產生的 偏差均不超過 ,則稱該方法是穩(wěn)定的。,二 單步法的穩(wěn)定性,0.025,0.050,0.075,0.100,1,2,3,4,5,0,-1,-2,-3,x,y,其方程時間常數(shù)為,因此有（10）知，要使Euler法穩(wěn)定 則步長,如果取步長h=0.025則Euler格式為：,其結果在準確解,上下波動不穩(wěn)定，再看后退Euler格式,H=0.025時，形式為：,計算結果穩(wěn)定。具體結果如下：,節(jié)點,Eule方法,后退Euler方法,0.025,-1.5,0.2857,0.050.,2.25,0.0816,0.075,-3.375,0.0233,0.100,5.0625,0.0067,前面介紹的幾種步進方法，在計算 時大多只用到前一個節(jié)點上 的近似值 ，而沒有用到前幾步的計算所得出的信息，故稱為單步 法 。實際上經過多次單步法計算后，已得出一系列近似值 等。為了充分利用這些信息來計算 以減少計算 工作量和獲得教高的精度，可采用如下計算公式：,5 線性多步法,一 顯式Adams法,（1）,例如 k=3 時有：,（7）,（7）稱為Adams四步顯式法。它用到了四個節(jié)點上的f值，是,一種最常用的多步法，其精度為四階。,如果利用,共k+1個數(shù)據來構造一個Newton內插多次式,，則與上面類似推導可得：,局部截斷誤差為：,二 Adams隱式公式, 在同一階數(shù)下，隱式的局部截斷誤差的常數(shù)的絕對值 比顯式的|Bk|要小。,三 Adzms 顯式與隱式的比較,當k=2時，,解：取h=0.1，兩種方法具體算式如下：,算法啟動值（對四步顯式當x1=0.1, x2 =0.2, x3 =0.3時，三步隱式，當x1 =0.1, x2 =0.2時）可用同階的RungeKutta公式計算。而本例則是由精確解 算出的。計算結果部分如下表,如下列定解問題： 其Ada