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

在工程和科學(xué)技術(shù)的實(shí)際問(wèn)題中，常需求解微分方程，但常微分方程中往往只有少數(shù)較簡(jiǎn)單和典型的常微分方程（例如線性常系數(shù)常微分方程等）可求出其解析解,對(duì)于變系數(shù)常微分方程的解析求解就比較困難，而一般的非線性常微分方程的求解困難就更不用說(shuō)了。大多數(shù)情況下，常微分方程只能用近似方法求解。這種近似解法可分為兩大類：一類是近似解析法，如級(jí)數(shù)解法、逐次逼近法等;另一類是數(shù)值解法，它給出方程在一些離散點(diǎn)上的近,第八章 常微分方程數(shù)值解法,1、 引 言,似值。,其中 x 是質(zhì)量，m是離開(kāi)平衡位o的距離，t為時(shí)間，c為彈簧系數(shù)。,在具體求解微分方程時(shí)，需具備某種定解條件，微分方程和定解條件合在一起組成定解問(wèn)題。定解條件有兩種：一種是給出積分曲線在初始點(diǎn)的狀態(tài)，稱為初始條件，相應(yīng)的定解問(wèn)題稱為初值問(wèn)題。另一類是給出積分曲線首尾兩端的狀態(tài)，稱為邊界條件，相應(yīng)的定解問(wèn)題稱為邊值問(wèn)題。,我們現(xiàn)在討論常微分方程的數(shù)值解法。先從最簡(jiǎn)單的一階常微分方程的初值問(wèn)題出發(fā)開(kāi)始討論。,由常微分方程理論可知：只要上式中的函數(shù)f(x,y)在區(qū)域 G=axb,y內(nèi)連續(xù)，且關(guān)于 y 滿足Lipschitz條件，即存在與 x, y 無(wú)關(guān)的常數(shù)L，使,至于初值問(wèn)題（1）的數(shù)值解法，常采用差分方法，即把一個(gè)連續(xù)的初值問(wèn)題離散化為一個(gè)差分方程來(lái)求解。具體地，將（1）離散化后，求找其解 y=y(x)在一系列離散點(diǎn),下面分析均假定滿足上述條件。,一、Euler公式,2 Euler方法,因?yàn)槌踔祮?wèn)題中的初始條件 已知，即可利用已知的 來(lái)求出下一節(jié)點(diǎn)處 的近似值 ；再?gòu)?來(lái)求 如此繼續(xù)，直到求出 為止。這種用按節(jié)點(diǎn)的排列順序一步一步地向前推進(jìn)的方式求解的差分算法稱為“步進(jìn)式”或“遞推式”算法，它是初值問(wèn)題數(shù)值解法的各種差分格式的共同特點(diǎn)。因此，只要能寫(xiě)出由前幾步已知信息 來(lái)計(jì)算 的遞推公式（即差分格式），即可完全表達(dá)這種算法。,若將 和 的近似值分別記為 和 ， 則得： （3）,這就是Euler公式（格式）。 利用它可由初值 出發(fā)逐步算出 。 這類形式的方法也稱為差分方法。,定義：如果局部截?cái)嗾`差為 ，則這種數(shù)值算法的精度為p階，故Euler格式的精度為一階。 從幾何意義上來(lái)看，如圖，,當(dāng)假定 為準(zhǔn)確值，即在 的前提下來(lái)估計(jì)誤差 ，這種截?cái)嗾`差稱為局部截?cái)嗾`差。,由（2）、（3）知Euler公式在 處的局部截?cái)嗾`差為：,由方程（1）知，其積分曲線 y=f(x)上任一點(diǎn)(x, y)的切線斜率 都等于函數(shù)f(x, y)的值。從初始點(diǎn) （即 點(diǎn)）出發(fā)，作積分 曲線 y=y(x) 在 點(diǎn)上的切線 （其斜率為 ）與直線 相交與 點(diǎn)（即 點(diǎn)），得到 作為 的近似值，則有,相比較知，這時(shí)用切線 近似代替了曲線段 點(diǎn)近似 代替了 點(diǎn)， 近似代替了 近似代替了 。 遞推繼續(xù)從 點(diǎn)出發(fā)，作一斜率為 的直線 與直線 相交于 點(diǎn)（即 點(diǎn)），得到 作為 的近似值 。如此直到 點(diǎn)。這樣得出一條折線 近似代替積分曲線 ，當(dāng) 步數(shù)越多時(shí)，由于誤差的積累，折線 可能會(huì)越,解：為便于進(jìn)行比較，我們后面將用多種數(shù)值方法求解上述初值問(wèn)題。 這里先用Euler 公式，此處具體格式為： 取步長(zhǎng)為h=0.1，計(jì)算結(jié)果略。 由結(jié)果可見(jiàn)Euler方法的精度很差。,即為Euler格式（3）。,因?yàn)椴钌淌俏⒎值慕?，所以Euler格式也可用差商近似代替導(dǎo)數(shù)的離散方法來(lái)得到。在節(jié)點(diǎn) 處有：,二 后退Euler格式,顯然Euler格式具有遞推性，在計(jì)算 時(shí)只要用到前一步所得的結(jié)果 一個(gè)信息就夠了，因此是一種單步格式或稱一步格式。 若用不同的數(shù)值微分計(jì)算方法也可導(dǎo)出其它形式的算法。 例如：用向后差商表示的數(shù)值微分公式,（6）稱為向后Euler公式，又稱為隱式Euler公式(后退Euler格式)。后退Euler公式與Euler公式有著本質(zhì)的區(qū)別，后者是關(guān)于 的一個(gè)直接計(jì)算公式，這類公式稱作顯式的；而前者，即（6）中右端含有未知的 ，它實(shí)際上是關(guān)于 的一個(gè)函數(shù)方程。這類公式稱作隱式的。,顯式與隱式兩類方法各有特點(diǎn)，使用顯式算法遠(yuǎn)比隱式算法方便，但考慮數(shù)值穩(wěn)定性等因素，人們常選用隱式算法。,隱試算法（6）常用迭代法來(lái)實(shí)現(xiàn)，而迭代過(guò)程實(shí)質(zhì)上是逐步 顯式化。,設(shè)用顯式Euler格式算出 作迭代初值 ，以此代入（6）右端，使之轉(zhuǎn)化為顯示，直接計(jì)算得： ，再用 代入（6）右端又有：,從幾何上看，梯形公式是取 區(qū)間兩端點(diǎn)處斜率的 平均斜率。,三 梯形公式,Euler方法是過(guò) 點(diǎn)以斜率 引直線交 的點(diǎn)A。 后退Euler方法是以點(diǎn) 處的斜率 為斜率，從 點(diǎn)引直線交 于另一點(diǎn)B。 A、B兩點(diǎn)都偏離點(diǎn) 點(diǎn)，梯形方法就是取A、B兩點(diǎn)的中點(diǎn)P作為Q2近似，上圖表明梯形方法確實(shí)改善了精度。,從而也可得二步Euler公式及其 截?cái)嗾`差為：,也可以由導(dǎo)數(shù)的中心差分近似式得到：,四 二步歐拉公式,將區(qū)間a,bn等分，得子區(qū)間 在第i+1 個(gè)子區(qū)間上，積分,例如：對(duì)于初值問(wèn)題：,對(duì)右端利用左矩形公式可得 即 Euler格式,各公式的截?cái)嗾`差可直接利用數(shù)值積分截?cái)嗾`差估計(jì)而得。從而可知梯形公式（8）的截?cái)嗾`差為：,梯形公式也是隱式的，可用迭代法求解，與后退Euler方法一樣，仍用Euler方法提供迭代初值，其迭代格式為：,為分析迭代過(guò)程的收斂性，將（12）與（8）相減得：,（12）,k=0，1，2，,L為f(x,y)關(guān)于y的Lipschitz常數(shù).如果選取h充分小使得 則 . 時(shí)有 這表明迭代過(guò)程（12）是收斂于（8）的解的。,五 改進(jìn)的Euler公式,上面已看到Euler公式計(jì)算量小但精度差，梯形方法雖然提高了 計(jì)算精度，但算法復(fù)雜計(jì)算量大，在應(yīng)用（12）進(jìn)行迭代時(shí)，每次 均要計(jì)算函數(shù)f的值，而迭代又要反復(fù)進(jìn)行多次，計(jì)算量很大，難以 預(yù)測(cè)。為了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入 下一步計(jì)算， 這就簡(jiǎn)化了計(jì)算。,具體地，先用Euler公式求得一個(gè)初步的近似值 稱之為預(yù)測(cè) 值。預(yù)測(cè)值 的精度可能很差。再用梯形公式（8）將它校正 一次，即按（12）式迭代一次得 ，這個(gè)結(jié)果稱為校正值，這 樣建立的 預(yù)測(cè)校正系統(tǒng)稱為改進(jìn)的Euler公式。,由表可見(jiàn)，與精確解 相比，改進(jìn)的Euler公式的精度較Euler 公式有明顯的提高。 下面再看兩步Euler公式（9），除了給出初值 外，還需要借助 其它單步法（如Euler公式，后退Euler公式及梯形公式等）再提供一個(gè),Euler公式,改進(jìn)的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,啟動(dòng)值 然后才能啟動(dòng)計(jì)算公式依次計(jì)算,用兩步Euler公式與梯形公式相匹配，又可得到下面預(yù)測(cè)-校正 系統(tǒng)：,校正：,（17）,（18）,兩步法優(yōu)美是由于它調(diào)用了兩個(gè)節(jié)點(diǎn)上的信息，從而能以較少的 計(jì)算量獲得較高的精度。,預(yù)測(cè)：,與改進(jìn)的,Euler公式（13）（14）相比較易見(jiàn)（17）（18）的一個(gè),突出特點(diǎn)是它的預(yù)測(cè)公式與校正公式具有同階精度。據(jù)此可以比較 方便的估計(jì)截?cái)嗾`差，并基于這種估計(jì)，可以提供一種提高精度的 簡(jiǎn)易方法。,再由梯形公式截?cái)嗾`差公式（11）知：校正公式（18）的截?cái)嗾`差 為： 比較（19）（20）可見(jiàn)，校正值的誤差 大約只是預(yù)測(cè)值 的誤差 的1/4(符號(hào)相反).即 ,由此導(dǎo)出誤差的事后估計(jì):,若預(yù)測(cè)公式（17）中的 和 都是準(zhǔn)確的。即 則由兩步Euler公式的截?cái)嗾`差公式（10）知：,1預(yù)測(cè): 2改進(jìn): 3計(jì)算:,可以期望利用這樣估計(jì)出的誤差作為計(jì)算結(jié)果的一種 補(bǔ)償,有 可能使精度得到改善. 以 和 分別表示第I步的預(yù)測(cè)值和校正值,按估計(jì)式(21), 和 分別作為 和 的改進(jìn) 值,在校正值 尚未求出之前,可用上一步的偏差 替代 來(lái)改進(jìn)預(yù)測(cè)值 .這樣設(shè)計(jì)的計(jì)算方案有如下六個(gè) 環(huán)節(jié):,3 Runge-Kutta方法,運(yùn)用上述方案計(jì)算 時(shí)，要用到先一步的信息 ， ， 和 更前一步的 。因此啟動(dòng)算法之步必須給出啟動(dòng)值 和 。 可用其它單步法（例如改進(jìn)的Euler方法）來(lái)計(jì)算， 則一般取為0。計(jì)算結(jié)果表明，這種簡(jiǎn)單的處理方法通常 可以獲得令人滿意的效果。,而具體則有:,式中y（x）的各階導(dǎo)數(shù)可由方程（*）用函數(shù)f來(lái)表達(dá)。引進(jìn)函數(shù) 序列 來(lái)描述求導(dǎo)過(guò)程：,二 Runge-Kutta公式的導(dǎo)出,例:用Taylor公式求解初值問(wèn)題:,K可看作是 y=y(x)在區(qū)間 上的平均斜率。從而Euler公式 相當(dāng)于取 點(diǎn)上的斜率 作為平均斜率K的 近似值，這當(dāng)然十分粗糙，因而精度必然很低。,這個(gè)過(guò)程啟示我們，如果設(shè)法在 內(nèi)多預(yù)測(cè)幾個(gè)點(diǎn)的斜率 值，然后將它們加權(quán)平均作平均斜率K的近似值，就有可能構(gòu)造 出更高精度的計(jì)算公式。這就是Runge-Kutta方法的基本思想。,根據(jù)預(yù)測(cè)值 再來(lái)算出 由此構(gòu)造出計(jì)算格式：,假定 ，分別將 和 作Taylor展開(kāi)得：,由（2）得：,（7）中三個(gè)待定 參數(shù) P ，但只有兩個(gè)方程，因此還有一個(gè)自由度。凡滿足條件（7）的一族格式統(tǒng)稱為二階Runge-Kutta格式。,當(dāng)p=1 ， 時(shí)，二階Runge-Kutta格式（6）即為改進(jìn)的 Euler格式（15）。 如取p= 1/2 ，則 ，二階R-K格式（6）成為：,稱之為變形的Euler格式 。,由于（8）中的 是由Euler格式預(yù)測(cè)出來(lái)的區(qū) 間 中的點(diǎn) 的近似解， 就近似地等于此中點(diǎn)的斜率 ，因此（8）就相當(dāng)于用中點(diǎn) 的 斜率作為（4） 中平均斜率K的近似值，故格式（8）也稱為中點(diǎn)格式。,總之，二階R-K格式用多算一次函數(shù)值f的辦法，避開(kāi)了二 階Taylor級(jí)數(shù)法所要求計(jì)算的 f的導(dǎo)數(shù)值，在這種意義上可以說(shuō)， R-K方法實(shí)質(zhì)上是Taylor方法的變形。,三 三階Runge-Kutta方法,斜率值 和 可以利用。我們用 和 線性組合給出區(qū)間 上的平均斜率，從而得到 的預(yù)測(cè)值,于是，再通過(guò)計(jì)算函數(shù)值f得到：,這樣設(shè)計(jì)出的計(jì)算格式具有形式：,（9）,我們希望適當(dāng)選擇系數(shù) 和p、q、r、 使上述公式具有三階精 度。為便于數(shù)學(xué)演算，引進(jìn)算子：,，,則根據(jù)（2）有：,于是三階展開(kāi),可表為：,式中 的下標(biāo) i 均表示在 處取值。,（12）,四 四階Rung-Kutta方法,最常用的一種經(jīng)典Rung-Kutta格式具體形式如下：,(13),解：三種方法如下：,Euler格式：,改進(jìn)的Euler格式：,(h=0.05),經(jīng)典的R-K格式：,(h=0.1),Euler法h=0.025,改進(jìn)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,這里采用了不同的步長(zhǎng)h值，是為了使他們所耗的計(jì)算工作量 大致相同，以便于比較。由上表可見(jiàn)，經(jīng)典的R-K方法的精確度 較改進(jìn)的Euler方法又有很大的提高。這一結(jié)論也可以從理論上 大致的分析出來(lái)。,析出來(lái)：Euler方法的局部截?cái)嗾`差為：,計(jì)算四步后的,而經(jīng)典R-K方法的局部截?cái)嗾`差則為：,可見(jiàn)，當(dāng),為大致相同數(shù)量級(jí)的常數(shù)時(shí)有：,但要注意的是：R-K方法的導(dǎo)出利用了Taylor展開(kāi)，因此要求 所求的解有教好的光滑性，如果解的光滑性差，則采用經(jīng)典的 R-K方法所的數(shù)值解，其精度有可能反而不及改進(jìn)的Euler方法， 因此在實(shí)際計(jì)算中應(yīng)根據(jù)問(wèn)題的具體情況來(lái)選擇適合的算法。,五 步長(zhǎng)的自動(dòng)選擇-變步長(zhǎng)的Runge-Kutta方法,在應(yīng)用數(shù)值法求解微分方程中，選擇適當(dāng)?shù)牟介L(zhǎng)是至關(guān)重要的。步 長(zhǎng)太大則達(dá)不到要求，步長(zhǎng)太小則步數(shù)增多，不但增加計(jì)算工作量， 還可能導(dǎo)致舍入誤差的嚴(yán)重積累。尤其是當(dāng)微分方程的解y(x)變化 激烈時(shí)，步長(zhǎng)的合理取法是在變化激烈處步長(zhǎng)取小些，在變化平緩 時(shí)取大些，也就是采取自動(dòng)變步長(zhǎng)的方法，即根據(jù)精度的要求先估 計(jì)出下一步長(zhǎng)的合理大小，然后按此計(jì)算。,作為近似值，則,的精確度都要高。,當(dāng)p=4時(shí)，可以取,這種修正方法與Romberg的數(shù)值積分的思路是一樣的。（15） 除以（14）得：,由此可以得出誤差,事后估計(jì)式：,由上面的分析可見(jiàn)，微分方程數(shù)值解的基本思想是，通過(guò),（1） 如果 ，則反復(fù)加倍步長(zhǎng)進(jìn)行計(jì)算，直到 時(shí)為止，并以上一次步長(zhǎng)的計(jì)算結(jié)果作為 （2）若 ， 則反復(fù)減半步長(zhǎng)進(jìn)行計(jì)算，直到 時(shí)為止，并取其最后一次步長(zhǎng)的計(jì)算作為 這樣做時(shí) ，為了選擇步長(zhǎng)，每一步都要反復(fù)判別， 增加 了工作量，但在方程的解 y(x) 變化劇烈的情況時(shí) ，總的計(jì)算 工作量可以得到減少，結(jié)果還是合算的。,這樣就可以從步長(zhǎng)減半前后的兩次計(jì)算結(jié)果的偏差 來(lái)判斷步長(zhǎng)選的是否適當(dāng)，當(dāng)要求的 數(shù)值精度為 時(shí)：,4 單步法的收斂性和穩(wěn)定性,一 單步法的收斂性,本例的Euler公式為 由此式遞推可得：,定義 ： 若一種數(shù)值方法對(duì)任意固定的 當(dāng) （同時(shí) ） 時(shí) 有 ，則稱方法是收斂的。,某種離散化手段，將微分方程轉(zhuǎn)化為差分方程來(lái)求解。這種轉(zhuǎn)化是 否合理，還要看差分方程問(wèn)題的解 當(dāng) 時(shí)是否會(huì)收斂到點(diǎn) 對(duì)固定的i將趨向 ，這時(shí)討論收斂是沒(méi)有意義的。,此定理表明，判斷單步法（1）的收斂性，歸結(jié)為驗(yàn)證增量函數(shù) 能否滿足Lipschitz（3） 對(duì)于Euler方法，由于其增量函數(shù) 故當(dāng)f滿足 Lipschitz 條件時(shí)它是收斂的。,穩(wěn)定性問(wèn)題比較復(fù)雜，為簡(jiǎn)化討論，僅考察下面的模型方程 為保證微分方程本身的穩(wěn)定性，假設(shè) 先考慮Euler方法的穩(wěn)定性。模型方程 的Euler公式為： 設(shè)在節(jié)點(diǎn)值 上有一擾動(dòng)值 , 它的傳播使節(jié)點(diǎn)值 產(chǎn)生大小為： 的擾動(dòng)值，假設(shè)用 按Euler 公式得出：,定義：若一種數(shù)值方法在節(jié)點(diǎn)值 上有擾動(dòng) ，而對(duì)于yi后的各個(gè) 節(jié)點(diǎn)值 上產(chǎn)生的 偏差均不超過(guò) ,則稱該方法是穩(wěn)定的。,二 單步法的穩(wěn)定性,0.025,0.050,0.075,0.100,1,2,3,4,5,0,-1,-2,-3,x,y,其方程時(shí)間常數(shù)為,因此有（10）知，要使Euler法穩(wěn)定 則步長(zhǎng),如果取步長(zhǎng)h=0.025則Euler格式為：,其結(jié)果在準(zhǔn)確解,上下波動(dòng)不穩(wěn)定，再看后退Euler格式,H=0.025時(shí)，形式為：,計(jì)算結(jié)果穩(wěn)定。具體結(jié)果如下：,節(jié)點(diǎn),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,前面介紹的幾種步進(jìn)方法，在計(jì)算 時(shí)大多只用到前一個(gè)節(jié)點(diǎn)上 的近似值 ，而沒(méi)有用到前幾步的計(jì)算所得出的信息，故稱為單步 法 。實(shí)際上經(jīng)過(guò)多次單步法計(jì)算后，已得出一系列近似值 等。為了充分利用這些信息來(lái)計(jì)算 以減少計(jì)算 工作量和獲得教高的精度，可采用如下計(jì)算公式：,5 線性多步法,一 顯式Adams法,（1）,例如 k=3 時(shí)有：,（7）,（7）稱為Adams四步顯式法。它用到了四個(gè)節(jié)點(diǎn)上的f值，是,一種最常用的多步法，其精度為四階。,如果利用,共k+1個(gè)數(shù)據(jù)來(lái)構(gòu)造一個(gè)Newton內(nèi)插多次式,，則與上面類似推導(dǎo)可得：,局部截?cái)嗾`差為：,二 Adams隱式公式, 在同一階數(shù)下，隱式的局部截?cái)嗾`差的常數(shù)的絕對(duì)值 比顯式的|Bk|要小。,三 Adzms 顯式與隱式的比較,當(dāng)k=2時(shí)，,解：取h=0.1，兩種方法具體算式如下：,算法啟動(dòng)值（對(duì)四步顯式當(dāng)x1=0.1, x2 =0.2, x3 =0.3時(shí)，三步隱式，當(dāng)x1 =0.1, x2 =0.2時(shí)）可用同階的RungeKutta公式計(jì)算。而本例則是由精確解 算出的。計(jì)算結(jié)果部分如下表,如下列定解問(wèn)題： 其Ada