第八章常微分方程組的數(shù)值解ppt課件

1、 在工程和科學(xué)技術(shù)的實(shí)際問題中，常需求解微分方程，但常微分方程中往往只有少數(shù)較簡單和典型的常微分方程例如線性常系數(shù)常微分方程等可求出其解析解,對(duì)于變系數(shù)常微分方程的解析求解就比較困難，而一般的非線性常微分方程的求解困難就更不用說了。大多數(shù)情況下，常微分方程只能用近似方法求解。這種近似解法可分為兩大類：一類是近似解析法，如級(jí)數(shù)解法、逐次逼近法等;另一類是數(shù)值解法，它給出方程在一些離散點(diǎn)上的近似值。其中 x 是質(zhì)量，m是離開平衡位o的距離，t為時(shí)間，c為彈簧系數(shù)。 例如：彈簧一質(zhì)量系統(tǒng)的振動(dòng)問題經(jīng)一定的簡化后可用一個(gè)二階常微分方程oxmcdtxd22來描述。 在具體求解微分方程時(shí)，需具備某種定解條

2、件，微分方程和定解條件合在一起組成定解問題。定解條件有兩種：一種是給出積分曲線在初始點(diǎn)的狀態(tài)，稱為初始條件，相應(yīng)的定解問題稱為初值問題。另一類是給出積分曲線首尾兩端的狀態(tài)，稱為邊界條件，相應(yīng)的定解問題稱為邊值問題。mxxoc :)()(000000成一個(gè)初值問題唯一確定。這就可以寫律都確定時(shí)，彈簧振動(dòng)規(guī)和初始速度時(shí)的初始位置當(dāng)彈簧在振動(dòng)開始時(shí)刻txxtxdtdxxtxtttt 我們現(xiàn)在討論常微分方程的數(shù)值解法。先從最簡單的一階常微分方程的初值問題出發(fā)開始討論。) 1 ()(),(),(0yaybaxyxfdxdy 由常微分方程理論可知：只要上式中的函數(shù)f(x,y)在區(qū)域G=axb,y內(nèi)連續(xù)，且

3、關(guān)于 y 滿足Lipschitz條件，即存在與 x, y 無關(guān)的常數(shù)L，使000022)()(0 xtxxtxxmcdtxd1, (0,1,1)iiihxx in兩相鄰點(diǎn)間的距離稱為步長。)1.(0niihxxhhii，有（常值）時(shí)稱為等步長當(dāng) 至于初值問題1的數(shù)值解法，常采用差分方法，即把一個(gè)連續(xù)的初值問題離散化為一個(gè)差分方程來求解。具體地，將1離散化后，求找其解 y=y(x)在一系列離散點(diǎn) bxxxxxani210niyyyyy210上的近似值 下面分析均假定滿足上述條件。2121),(),(yyLyxfyxf12,()Gyyyy x對(duì)內(nèi) 任 意 兩 個(gè)都 成 立 ， 則 方 程 (1)的

4、 解必 定 存 在 且 唯 一 。 對(duì)于初值問題1），先將其離散化，即把a(bǔ),b區(qū)間n等分，.), 2 , 1 , 0(nabhniihaxi其中得各離散節(jié)點(diǎn) 因?yàn)槌踔祮栴}中的初始條件 知，即可利用已知的 來求出下一節(jié)點(diǎn)處 的近似值 ；再從 來求 如此繼續(xù)，直到求出 為止。這種用按節(jié)點(diǎn)的排列順序一步一步地向前推進(jìn)的方式求解的差分算法稱為“步進(jìn)式或“遞推式算法，它是初值問題數(shù)值解法的各種差分格式的共同特點(diǎn)。因而，只要能寫出由前幾步已知信息 來計(jì)算 的遞推公式即差分格式），即可完全表達(dá)這種算法。0)(yay0y)(1xy1y1ynyiyyy10 1 iy.2y若將 和 的近似值分別記為 和 ，則得：

5、 （3）)(ixyiy) 1, 2 , 1 , 0(),(1niyxhfyyiiii)(,()()(1iiiixyxhfxyxy或) 1, 2 , 1 , 0().()()(1nixhyxyxyiii將上式改寫成：無窮小量充分小時(shí)，可忽略高階有界且當(dāng)),( 2)(2iiyhhy 展式為：處在為方程的解，則設(shè)Tayloryxxyxyyiii),()( 1)2()1,2,1 ,0()(),( 2)( )()(121nixxyhxhyxyxyiiiiiii)(1ixy1iy 這就是Euler公式格式）。 利用它可由初值 出發(fā)逐步算出 。 這類形式的方法也稱為差分方法。nyyy21iy)(iixyy

6、1)(1ixyyi) 1, 1 , 0(),(),(2)(1211nixxyhyxyRiiiiiii)4()(2hO定義：如果局部截?cái)嗾`差為定義：如果局部截?cái)嗾`差為 ，則這種數(shù)，則這種數(shù)值算法的精度為值算法的精度為p p階，故階，故EulerEuler格式的精度為一階。格式的精度為一階。 從幾何意義上來看，如圖，從幾何意義上來看，如圖， 0y 當(dāng)假定 為準(zhǔn)確值，即在 的前提下來估計(jì)誤差 ，這種截?cái)嗾`差稱為局部截?cái)嗾`差。ix)(1phO 由2）、（3知Euler公式在 處的局部截?cái)嗾`差為：由方程1知，其積分曲線 y=f(x)上任一點(diǎn)(x, y)的切線斜率 都等于函數(shù)f(x, y)的值。從初始點(diǎn)

7、（即 點(diǎn)動(dòng)身，作積分曲線 y=y(x) 在 點(diǎn)上的切線 （其斜率為 ）與直線 相交與 點(diǎn)即 點(diǎn)），得到 作為 的近似值，則有dxdy0p),(00yx10pp),(0oyxf),(11yx1y)(1xyynpOxnniixxxxxxx112100p1pip1ip1np1p2pip1ip1npnpy=f(x)0p1xx 1p2p),()(0001010yxhfydxdyxxyyoxx10)(,()()(01xxdxxyxfxyxy與相比較知，這時(shí)用切線 近似代替了曲線段 點(diǎn)近似代替了 點(diǎn)， 近似代替了 近似代替了 。 遞推繼續(xù)從 點(diǎn)出發(fā)，作一斜率為 的直線與直線 相交于 點(diǎn)即 點(diǎn)），得到 作為

8、的近似值 。如此直到 點(diǎn)。這樣得出一條折線 近似代替積分曲線 ，當(dāng)步數(shù)越多時(shí)，由于誤差的積累，折線 可能會(huì)越10pp110, ppp1p1y),(),(001yxhfxy10)(,(xxdxxyxf1p),(11yxf21pp2p),(22yx2y)(2xynpnippppp210210nipppppnippppp2101xx 曲線。來越偏離真解210nippppp1)0() 10(2yxyxyy例：求解初值問題 解：為便于進(jìn)行比較，我們后面將用多種數(shù)值方法求解上述初值問題。 這里先用Euler 公式，此處具體格式為：取步長為h=0.1，計(jì)算結(jié)果略。 由結(jié)果可見Euler方法的精度很差。1 .

9、 0, 1)0(, 0)2(001hyyxyxyhyynnnnn即為Euler格式3）。 )5(), 2 , 1 , 0()(,()(nixyxfxyiii) 1, 2 , 1 , 0()()()(1nixyhxyxyiii用向前微商) 1, 2 , 1 , 0(),(,()()(51nixyxhfxyxyiiii）可得：代入（則得分別用其近似值代替，)(),(1iixyxy) 1, 2 , 1 , 0(),(1niyxhfyyiiiiix 因?yàn)椴钌淌俏⒎值慕?，所以Euler格式也可用差商近似代替導(dǎo)數(shù)的離散方法來得到。在節(jié)點(diǎn) 處有： 顯然Euler格式具有遞推性，在計(jì)算 時(shí)只要用到前一步所得

10、的結(jié)果 一個(gè)信息就夠了，因此是一種單步格式或稱一步格式。 若用不同的數(shù)值微分計(jì)算方法也可導(dǎo)出其它形式的算法。 例如：用向后差商表示的數(shù)值微分公式1iyiy) 1, 2 , 1 , 0()(2)()(1)(11niyhxyxyhxyiiii)(21)()()()(21111iiiiiiiixxyxxxyxyxyTaylor展式由可得：代入)(,()(111iiixyxfxy12111),(2)(,()()(iiiiiiiixxfhxyxhfxyxy和局部截?cái)嗾`差：則可得近似計(jì)算公式：)6()1,2, 1 ,0(),(111niyxhfyyiiii)7()1,2, 1 ,0()()(),(2212

11、nihoxxyhRiiiii （6稱為向后Euler公式，又稱為隱式Euler公式(后退Euler格式)。后退Euler公式與Euler公式有著本質(zhì)的區(qū)別，后者是關(guān)于 的一個(gè)直接計(jì)算公式，這類公式稱作顯式的；而前者，即6中右端含有未知的 ，它實(shí)際上是關(guān)于 的一個(gè)函數(shù)方程。這類公式稱作隱式的。1ny1ny1ny 顯式與隱式兩類方法各有特點(diǎn)，使用顯式算法遠(yuǎn)比隱式算法方便，但考慮數(shù)值穩(wěn)定性等因素，人們常選用隱式算法。 隱試算法6常用迭代法來實(shí)現(xiàn)，而迭代過程實(shí)質(zhì)上是逐步顯式化。 設(shè)用顯式Euler格式算出 作迭代初值 ，以此代入6右端，使之轉(zhuǎn)化為顯示，直接計(jì)算得： ，再用 代入6右端又有：),()0(

12、1iiiiyxhfyy)0(1iy(0)(0)111(,)iiiiyyhf xy)1(1iy),()1(11)2(1iiiiyxhfyy若迭代過程收斂，則極限值 必為隱式方程6）的解，從而可獲得后退Euler方法的解。如此反復(fù)進(jìn)行可得序列 。)(1hiy, 2 , 1 , 0),()(11)1(1kyxhfyykiiiki1)(1limikikyy) 8 (),(),(2111iiiiiiyxfyxfhyy 從幾何上看，梯形公式是取 區(qū)間兩端點(diǎn)處斜率的平均斜率。,1iixx)(2iyh比較Euler公式和后退歐拉公式的截?cái)嗾`差公式4），（7易見兩者為同階，但符號(hào)相反，因此就想到將這兩種算法進(jìn)行

13、算術(shù)平均，其結(jié)果可能消除誤差的主要部分 ，而獲得更高精度的結(jié)果。這種平均方法稱為梯形方法，公式為：1ixx Euler方法是過 點(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í)改善了精度。),(iiyxf),(112iiyxQ),(11iiyxf1Q1ixx2Q1iixxy=f(x)xyoAB2Q1QP),(1iiyxQ)2 , 1 , 0()(,()(nixyxfxyiii代入)(3)(,(2)()(311iiiiiyhxyxhfxyxy得從而也可

14、得二步Euler公式及其 截?cái)嗾`差為：)10() 1, 1 , 0()()(3)9() 1, 1 , 0(),(23311nihoyhRniyxhfyyiiiiii), 2 , 1 , 0()(6)()(21)(211niyhxyxyhxyiiii也可以由導(dǎo)數(shù)的中心差分近似式得到：將區(qū)間a,bn等分，得子區(qū)間在第i+1 個(gè)子區(qū)間上，積分 1, 1 , 0,1nixxii11),()()()(1iiiixxixxiidxyxfdxxyxyxyEuler公式、后退Euler公式以及梯形公式，二步Euler公式均可利用前面的數(shù)值積分技術(shù)得到。axyxybaxyxfdxdy000)(,),(例如：對(duì)于

15、初值問題：對(duì)右端利用左矩形公式可得即 Euler格式)(,()()(1iiiixyxhfxyxy),(1iiiiyxhfyy式可得：上進(jìn)行，利用中矩形公間如果將積分改成在子區(qū),11iixx格式。二步Euleryxhfyyiiii),(211格式隱式即可得：對(duì)右端利用右矩形公式Euleryxhfyyxyxhfxyxyiiiiiiii),().(,()()(111111梯形公式式可得：對(duì)右端利用梯形積分公),(),(2111iiiiiiyxfyxfhyy各公式的截?cái)嗾`差可直接利用數(shù)值積分截?cái)嗾`差估計(jì)而得。從而可知梯形公式8的截?cái)嗾`差為：)11()() 11 , 0()(12313honixxyhR

16、iiiii梯形公式也是隱式的，可用迭代法求解，與后退Euler方法一樣，仍用Euler方法提供迭代初值，其迭代格式為：為分析迭代過程的收斂性，將12與8相減得：（12）k=0，1，2，.L為f(x,y)關(guān)于y的Lipschitz常數(shù).如果選取h充分小使得那么 . 時(shí)有 這表明迭代過程12是收斂于8的解的。 kiiiiikiiiiiyxfyxfhyyyxhfyy111101,2,12hlk 11ikiyy 011111111111122,2iikkiikiiiikiiyyhlyyhlyxfyxfhyy這一公式也可改寫為： ),(,(),(21iiiiiiiiyxhfyhxfyxfhyy（15）)

17、,(),(2111iiiiiiyxfyxfhyy(14)預(yù)測 （13） 校正),(1iiiiyxhfyy 上面已看到Euler公式計(jì)算量小但精度差，梯形方法雖然提高了計(jì)算精度，但算法復(fù)雜計(jì)算量大，在應(yīng)用12進(jìn)行迭代時(shí)，每次均要計(jì)算函數(shù)f的值，而迭代又要反復(fù)進(jìn)行多次，計(jì)算量很大，難以預(yù)測。為了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入 下一步計(jì)算，這就簡化了計(jì)算。 具體地，先用Euler公式求得一個(gè)初步的近似值 稱之為預(yù)測值。預(yù)測值 的精度可能很差。再用梯形公式8將它校正一次，即按12式迭代一次得 ，這個(gè)結(jié)果稱為校正值，這樣建立的 預(yù)測校正系統(tǒng)稱為改進(jìn)的Euler公式。1iy1iy1iy 式：上式表為下

18、列平均化形 ),(),()(2111iiippiiccpiyxhfyyyxhfyyyyy（16） 歐拉公式和改進(jìn)歐拉公式分iiipipicicpiyyyyyyyyyyyy9 .0)(1 .091.0)(1 .0905.0211及兩種格式的計(jì)算結(jié)果分別列表如下別為：000,1xy解：取步長h=0.1 iiiyyyyi9 . 0)( 1 . 01ydxdyy1)0( 1 , 0 x并比較兩法所得計(jì)算結(jié)果的精度。例：試分別用函數(shù)公式和改進(jìn)的歐拉公式求解： 由表可見，與精確解 相比，改進(jìn)的Euler公式的精度較Euler公式有明顯的提高。 下面再看兩步Euler公式9），除了給出初值 外，還需要借助其

19、它單步法如Euler公式，后退Euler公式及梯形公式等再提供一個(gè)xey oyixEuler公式改進(jìn)的Euler公式精確解iyiy iiyxy iiyxy ixy01110.10.20.310.90000000.81000000.72900000.34867843108374. 43107308. 821018182. 121092010. 10.90500000.81902500.74121760.36854100.90483740.81873080.74081820.3678794410626. 1410942. 2410994. 3410616. 6啟動(dòng)值 然后才能啟動(dòng)計(jì)算公式依次計(jì)算1

20、y32yy 用兩步Euler公式與梯形公式相匹配，又可得到下面預(yù)測-校正系統(tǒng)： 校正：),(),(2111iiiiiiyxfyxfhyyiiiiyxhfyy,211（17）（18）兩步法優(yōu)美是由于它調(diào)用了兩個(gè)節(jié)點(diǎn)上的信息，從而能以較少的計(jì)算量獲得較高的精度。 預(yù)測：與改進(jìn)的 Euler公式13）（14相比較易見17）（18的一個(gè)突出特點(diǎn)是它的預(yù)測公式與校正公式具有同階精度。據(jù)此可以比較方便的估計(jì)截?cái)嗾`差，并基于這種估計(jì)，可以提供一種提高精度的簡易方法。 )19()(3)( 311iiixyhyxy)(11iixyyiy1iy iixyy )20()(12)( 311iiixyhyxy11)(i

21、iyxy11)(iiyxy41)()(1111iiiiyxyyxy再由梯形公式截?cái)嗾`差公式11知：校正公式18的截?cái)嗾`差為： 比較19）（20可見，校正值的誤差 大約只是預(yù)測值 的誤差 的1/4(符號(hào)相反).即 ,由此導(dǎo)出誤差的事后估計(jì): 若預(yù)測公式17中的 和 都是準(zhǔn)確的。即 則由兩步Euler公式的截?cái)嗾`差公式10知：)21()(51)()(54)(11111111iiiiiiiiyyyxyyyyxy 1預(yù)測: 2改進(jìn): 3計(jì)算: *112iiihyyp)(5411iiiicppm),(11*1iiimxfm可以期望利用這樣估計(jì)出的誤差作為計(jì)算結(jié)果的一種 補(bǔ)償,有可能使精度得到改善. 以

22、和 分別表示第I步的預(yù)測值和校正值,按估計(jì)式(21), 和 分別作為 和 的改進(jìn)值,在校正值 尚未求出之前,可用上一步的偏差 替代 來改進(jìn)預(yù)測值 .這樣設(shè)計(jì)的計(jì)算方案有如下六個(gè)環(huán)節(jié):ipic)(54111iiicpp)(51111iiicpc1ip1ic11iicp1ip1iciicp 4校正: 5改進(jìn):6計(jì)算：)(2*11iiiiymhyc111151iiiicpcy),(11*1iiiyxfy有解y=yx）。Taylor展開有：00)(),(yxyyxfy(*)運(yùn)用上述方案計(jì)算 時(shí)，要用到先一步的信息 ， ， 和 更前一步的 。因此啟動(dòng)算法之步必須給出啟動(dòng)值 和 。 可用其它單步法例如改進(jìn)

23、的Euler方法來計(jì)算， 則一般取為0。計(jì)算結(jié)果表明，這種簡單的處理方法通?？梢垣@得令人滿意的效果。1iy1iy1y11cp 1y11cp *iiiiyypc) 1 ()( ! 3)( ! 2)( )()(321nnnnnxyhxyhxhyxyxy而具體則有: )1()2()2(jjjjfxffxfy)2()1()1()1()0()0()0( fyffxfyfyffxfyffy, 3 , 2j),()(yxfi式中yx的各階導(dǎo)數(shù)可由方程（*）用函數(shù)f來表達(dá)。引進(jìn)函數(shù)序列 來描述求導(dǎo)過程：)2()(2 222222yfxfyfyffyxffxfyyffxfyfy),(nnyx)()(njxy若在

24、(1)中右端取前面若干項(xiàng),并且在 處按(2)計(jì)算系數(shù) 的 近似值 的近似值 , 可導(dǎo)出如下Taylor公式: 其中p=1時(shí)一階Taylor展開式即為 而提高Taylor公式的階p即可提高計(jì)算結(jié)果精度.P階Taylor展開的局部截?cái)嗾`差為:)3(! 2)( 21pnpnnnnyphyhhyyy),(1nnnnnnyxhfyhyyyEuler公式 jny)()()!1()(1)1(111pnpnnhoyphyxy)(1nnxx因而具有P階精度: 對(duì)于一階常微分方程1的解y=y(x)，可利用中值定理得： iiiiiixxhxyxyxy11，01i即 hxhyxyxyiiii1也 即 )4(1hKxy

25、xyii式中K=hxyhxfhxyiiiiii,（5）1)0(1020yyxyxyy例:用Taylor公式求解初值問題:K可看作是 y=y(x)在區(qū)間 上的平均斜率。從而Euler公式相當(dāng)于取 點(diǎn)上的斜率 作為平均斜率K的近似值，這當(dāng)然十分粗糙，因而精度必然很低。1,iiyxiiyx ,iixxyxfdxdyi,再考慮改進(jìn)Euler公式15可改寫成 ：1, 2 , 1 , 0,21121211nihKyxfKyxfKKKhyyiiiiii 和 兩個(gè)點(diǎn)上的斜率 和 的算術(shù)平均值作為4中平均斜率K的近似值。其中 是 通過 知與4比較可見，它相當(dāng)于把iiyx ,11,iiyx1K2K2K信息 來近似

26、地預(yù)測的 。iy這個(gè)過程啟示我們，如果設(shè)法在 內(nèi)多預(yù)測幾個(gè)點(diǎn)的斜率值，然后將它們加權(quán)平均作平均斜率K的近似值，就有可能構(gòu)造出更高精度的計(jì)算公式。這就是Runge-Kutta方法的基本思想。1,iixx的如今，設(shè)想取區(qū)間 內(nèi)某一節(jié)點(diǎn) 上斜率 與 點(diǎn)上的 斜率 作線性組合即加權(quán)平均化為平均斜率K的近似，即4化為：1,iixx10pphxxipi2Kiiyx ,1K)(22111KKhyyii為了要得到 點(diǎn)上的斜率 ，需先預(yù)測Pipiyx,2K1phKyyipi根據(jù)預(yù)測值 再來算出 由此構(gòu)造出計(jì)算格式：12,phKyphxfyxfKiipipipiy1, 1 , 0,12122111niphKyph

27、xfKyxfKKKhyyiiiiIi（6）上式中含三個(gè)待定系數(shù) 和p，適當(dāng)選定它們以使算法的局部誤差為 ，從而具有二階精度。21)(3ho假定 ，分別將 和 作Taylor展開得：)(iixyy 1K2K)(),(),(),(),()( ),(21121hoyxfKyxfphyxfphKyphxfKxyyxfKiiyiixiiiiiii由2得：)()( )( 2hoxphyxyiiix將此式與 y(x) 在處的二階Taylor展開在 處的取值相比較：)()( 2)( )()(321hoxyhxhyxyxyiiii1ix代入6得：)()( )( )()(322211hoxyphxyhxyyiii

28、i成立，那么6的局部截?cái)嗾`差就等于 ，從而能具有二階精度。由系數(shù)比較知，只要：12211(7)2p)(3ho （7中三個(gè)待定 參數(shù) P ，但只有兩個(gè)方程，因此還有一個(gè)自由度。凡滿足條件7的一族格式統(tǒng)稱為二階Runge-Kutta格式。 21 當(dāng)p=1 ， 時(shí)，二階Runge-Kutta格式6即為改進(jìn)的Euler格式15）。 如取p= 1/2 ，那么 ，二階R-K格式6成為：21211021)8() 1, 1 , 0()2,2(),(12121niKhyhxfKyxfKhKyyiiiiii稱之為變形的Euler格式 。 由于8中的 是由Euler格式預(yù)測出來的區(qū)間 中的點(diǎn) 的近似解， 就近似地等

29、于此中點(diǎn)的斜率 ，因而8就相當(dāng)于用中點(diǎn) 的 斜率作為4）中平均斜率K的近似值，故格式8也稱為中點(diǎn)格式。1212Khyyii1,iixx21ix2K)(,(2121iixyxf21ix 總之，二階R-K格式用多算一次函數(shù)值f的辦法，避開了二階Taylor級(jí)數(shù)法所要求計(jì)算的 f的導(dǎo)數(shù)值，在這種意義上可以說，R-K方法實(shí)質(zhì)上是Taylor方法的變形。并用三個(gè)點(diǎn) 的斜率值 線性組合得到平均斜率K，這時(shí)計(jì)算公式為：qipiixxx321KKK3322111KKKhyyii其中 仍用公式6所取的形式。 為了預(yù)測點(diǎn) 的 斜率值 ，在區(qū)間 內(nèi)有兩個(gè)21KKpix3K,piixx1i qixx qh p q p

30、iixx為了進(jìn)一步提高精度，設(shè)除了 外，再考慮一點(diǎn) 表面上， 只含一個(gè)斜率值 ， 但 要通過 21hKyyii2K2K1K才能算出來，因此式中隱含著 ，每完成一步仍然需計(jì)算f的兩次函數(shù)值，其工作量仍與改進(jìn)的Euler格式一樣。1K斜率值 和 可以利用。我們用 和 線性組合給出區(qū)間 上的平均斜率，從而得到 的預(yù)測值1K2K1K2K,piixx)(qixy)(21KrKqhyyiqi于是，再通過計(jì)算函數(shù)值f得到：)(,(),(213KrKqhyqhxfyxfKiiqiqi這樣設(shè)計(jì)出的計(jì)算格式具有形式：)(,(),(),()(2131213322111KrKqhyqhxfKphKyphxfKyxfK

31、KKKhyyiiiiiiii（9）我們希望適當(dāng)選擇系數(shù) 和p、q、r、 使上述公式具有三階精度。為便于數(shù)學(xué)演算，引進(jìn)算子：321yfxD22222222yfyxfxD，則根據(jù)2有：DfyffDyDfyfy2 ，于是三階展開 3 2162iiiiiyhyhhyyy可表為：)10()(622321iyiiiiDfffDhDfhhfyy式中的下標(biāo) i 均表示在處取值。iiDff),(iiyxiiifyxfk,1另一方面，由于.)2(2.),(122232iyiiifDfpqfDqhqhDffKDDyfDfhpsDDr則則若取2231212()(,()2()iiiiiqhkf xqh yqh rkqk

32、fqhDfD fDrkqkxy其中的算子.2,2212iiiiiifDphphDfyxfphkyqhxfk將3,21,kkk的展開式帶入9），再令1321可得 ：iiiiDfqphfhyy3221322222331122yihp Dfq Dpqf D f此式與10式比較系數(shù)可得：2322233121316pqpqp q以及1231,1r稱滿足11）.的公式族9為三階Runge-Kutta公式（11）132121,1,1,2632pqrR ungeK uttaK utta 當(dāng)時(shí)三 階公 式 稱 為格 式 。（12）)2,()2,2(),()4(62131213211hkhkyhxfkkhyhxf

33、kyxfkkkkhyyiiiiiiii其形式為：最常用的一種經(jīng)典Rung-Kutta格式具體形式如下：公式。各種四階導(dǎo)出復(fù)雜的數(shù)學(xué)演算。可以繼續(xù)上述過程，經(jīng)過較KuttaRunge),()2,2()2,2(),()22(6342312143211hkyhxfkkhyhxfkkhyhxfkyxfkkkkkhyyiiiiiiiiii(13)四階Rung-Kutta方法每一步需要四次計(jì)算函數(shù)值f,可以證明其截?cái)嗾`差為50 h不過證明復(fù)雜。例 ： 試分別用Euler方法h=0.025改進(jìn)的Eulerh=0.05及經(jīng)典R-K方法h=0.1求解下列初值問題，比較三種方法的結(jié)果的精度。1)0(1,0,yxy

34、dxdy解：三種方法如下： Euler格式：)025. 0(,975. 0025. 01hiyiyiyiy改進(jìn)的Euler格式：icpiipiciiipyyyyyyyyyyyy95125.0)(219525.0)(05.095.0)(05.01(h=0.05)經(jīng)典的R-K格式：iiiiiiiiiiiykkkkyyykykykykyyykyk9048375.0)22(61.090475.0)1.0(9525.021.095.0)(21.043211342321(h=0.1)Euler法h=0.025改進(jìn)Eulerh=0.05R-K法h=0.1 精確解yiyiiyxy)(iyiiyxy)(iiyx

35、y)(iyx=0.1x=0.2x=0.331.1495 100.9036878900.8166518030.7379983450.9048765620.8188015930.740914370.9048375000.8187309010.7408180.904837480.818730750.7408182232 .0790 1032.820 105109144. 3510.084. 75106151. 98102 . 871048. 171002. 2x這里采用了不同的步長h值，是為了使他們所耗的計(jì)算工作量大致相同，以便于比較。由上表可見，經(jīng)典的R-K方法的精確度較改進(jìn)的Euler方法又有很

36、大的提高。這一結(jié)論也可以從理論上大致的分析出來。析出來：Euler方法的局部截?cái)嗾`差為：/22111016EECRC h計(jì)算四步后的212110410164CCRE而經(jīng)典R-K方法的局部截?cái)嗾`差則為：525210ChCRRR可見，當(dāng)2,1CC為大致相同數(shù)量級(jí)的常數(shù)時(shí)有：ERRR但要注意的是：R-K方法的導(dǎo)出利用了Taylor展開，因此要求所求的解有教好的光滑性，如果解的光滑性差，則采用經(jīng)典的R-K方法所的數(shù)值解，其精度有可能反而不及改進(jìn)的Euler方法，因此在實(shí)際計(jì)算中應(yīng)根據(jù)問題的具體情況來選擇適合的算法。 在應(yīng)用數(shù)值法求解微分方程中，選擇適當(dāng)?shù)牟介L是至關(guān)重要的。步長太大則達(dá)不到要求，步長太小

37、則步數(shù)增多，不但增加計(jì)算工作量，還可能導(dǎo)致舍入誤差的嚴(yán)重積累。尤其是當(dāng)微分方程的解y(x)變化激烈時(shí)，步長的合理取法是在變化激烈處步長取小些，在變化平緩時(shí)取大些，也就是采取自動(dòng)變步長的方法，即根據(jù)精度的要求先估計(jì)出下一步長的合理大小，然后按此計(jì)算。 從節(jié)點(diǎn) 動(dòng)身，先以h為 步長跨一步到節(jié)點(diǎn)1ix 值 如果計(jì)算公式是P階的，那么 當(dāng)h的變化不大時(shí)，式中的系數(shù)C可以近似的看作為常數(shù)。 然后將步長減半，即以 為步長，從節(jié)點(diǎn) 動(dòng)身，跨兩步到節(jié)點(diǎn) ，再求得一個(gè)近似值 。 其中每跨一步的截?cái)嗾`差為：)(1hiy21)(11pphixhochyyi1ix12phc故有：1(/ 2 )2112()2phhph

38、iiyxyco(14)（15）ixix，求一個(gè)近似 下面介紹一種Richrdson外推法：2h( /2)1hiy122111()22112111(1 5 )21 4221()221221iPhhppPiiihhppiixphhpiiipyyyxOhyyyohyyy得及若 取作為近似值，那么 hiiyy11 顯然比 的精確度都要高。 當(dāng)p=4時(shí)，可以取2()( )1111615hhiiiyyy這種修正方法與Romberg的數(shù)值積分的思路是一樣的。（15）除以14得：211()1()112hiiixhpixyyyy由此可以得出誤差事后估計(jì)式：22()()( )1111112hhhiiiipy xy

39、yy)(12hiy和 由上面的分析可見，微分方程數(shù)值解的基本思想是，通過（1） 假如 ，則反復(fù)加倍步長進(jìn)行計(jì)算，直到 時(shí)為止，并以上一次步長的計(jì)算結(jié)果作為 （2假設(shè) ， 則反復(fù)減半步長進(jìn)行計(jì)算，直到 時(shí)為止，并取其最后一次步長的計(jì)算作為這樣做時(shí) ，為了選擇步長，每一步都要反復(fù)判別， 增加了工作量，但在方程的解 y(x) 變化劇烈的情況時(shí) ，總的計(jì)算工作量可以得到減少，結(jié)果還是合算的。1iy1iy這樣就可以從步長減半前后的兩次計(jì)算結(jié)果的偏差 來判斷步長選的是否適當(dāng)，當(dāng)要求的數(shù)值精度為 時(shí)：2()()11hhiiyy110000,1,01,0.ihhiixiixxxxihyyhhhehyy e由于

40、:則而時(shí)差分方程的解時(shí)確實(shí)收斂到原微分方程的準(zhǔn)確解：下面進(jìn)一步考慮一般的單步法的收斂性。本例的Euler公式為 由此式遞推可得： iiyhy111011iihyy定義定義 ： 若一種數(shù)值方法對(duì)任意固定的若一種數(shù)值方法對(duì)任意固定的 當(dāng)當(dāng)（同時(shí)（同時(shí) ） 時(shí)時(shí) 有有 ，則稱方法是收斂的。，則稱方法是收斂的。nhxxn0)(nnxyy 0hn某種離散化手段，將微分方程轉(zhuǎn)化為差分方程來求解。這種轉(zhuǎn)化是否合理，還要看差分方程問題的解 當(dāng) 時(shí)是否會(huì)收斂到點(diǎn) 對(duì)固定的i將趨向 ，這時(shí)討論收斂是沒有意義的。ihxxi00hiy0 xoyyyy)0(考察Euler方法的收斂性 例如：對(duì)初值問題所謂的單步法，就是

41、在計(jì)算 時(shí)只用到它前一步的信息 。 Taylor級(jí)數(shù)法，Runger-Kutta方法等都是單步法的例子。顯然單步法的共同特征是，他們都是將 加上某種形式的增量得出 ，其計(jì)算公式形如： 式中的 稱作增量函數(shù)。例如： 對(duì)Euler公式，有對(duì)改進(jìn)的Euler公式有：關(guān)于單步法有下述的收斂定理： ),(1hyxhyyiiiihyx,iy1iyyxfhyx,),(,(,21,yxhfyhxfyxfhyx(1)(2)定理：設(shè)單步法定理：設(shè)單步法1具有具有p階的精度階的精度 ，且增量函數(shù)，且增量函數(shù) 關(guān)于關(guān)于y的滿足的滿足Lipschitz條件條件又設(shè)初值又設(shè)初值 是準(zhǔn)確的是準(zhǔn)確的 ，即，即y0=y(x0)

42、,則其的整體截?cái)嗾`差為則其的整體截?cái)嗾`差為:pnxhoyyn0y,pxy hxy hLyy(3)證明證明:設(shè)以設(shè)以 表示取表示取 用公式用公式1求的結(jié)果，求的結(jié)果，即即 則由于所給的方法具有則由于所給的方法具有P級(jí)精度，根據(jù)精度的定義級(jí)精度，根據(jù)精度的定義可知，有可知，有常數(shù)常數(shù)C使使 由由4與與1得：得：1iy oixyy hyxhxyyixiii),(1,11)(1pixchyyi(4)0(1)(1)1 ,(6)piiipppChehLehLL由此式遞推可得：111111111(,(),)(,)31)()()(1)()iiiiiiiipiiiiiiiippiiyyyxyhxy xhxyhh

43、Lyxyy xyyyy xyhLy xych由 假 設(shè) 條 件 （ ） 可 知 : （從 而 有 ：有如下的遞推關(guān)系：即整體截?cái)嗾`差：iiiyxye)(11111221(1)(1)(1)1(1) (1)(1)(1)1 ,(5)ppipipippppppippippChehLeChhLehLLhLhLeChChChhLehLL.故結(jié)論的證時(shí)則上式變?yōu)椋簳r(shí)，即由此可知，當(dāng)初值準(zhǔn)確, 00e1pT LppipceehohL0(1)(,1,101)ppiihLT LipnxnxxxihThLeexRxexxe 注 意 到 當(dāng) ：時(shí) 有 ：有時(shí) , 01 ,(7)pppTLTLipchee eeL從而有

44、： 此定理表明，判斷單步法1的收斂性，歸結(jié)為驗(yàn)證增量函數(shù) 能否滿足Lipschitz3） 對(duì)于Euler方法，由于其增量函數(shù) 故當(dāng)f滿足 Lipschitz 條件時(shí)它是收斂的。( , , )( , )x y hf x y(12oophh hyLipschitzLipschitzhLLL若限定為定數(shù)）則上式表明關(guān)于 滿足條件，且常數(shù)為因此改進(jìn)的Euler方法 也是收斂 的。類似的可以證明其 它的Runger-Kutta方法的收斂性。( , , )( , )(,( , )Eulerx y hf x yf xh yhf x y對(duì)于改進(jìn)的法，其增量函數(shù):則：yxhfyhxfyxhfyhxfyxfyxf

45、Ryxhyx,21),(,則由上式可得：常數(shù)為條件，且滿足關(guān)于若LLipschitzLipschitzyf( , , )( , , )12hLx y hx yhLyy上面關(guān)于收斂性的 討論有前提條件，即必須假定數(shù)值方法本身的計(jì)算是準(zhǔn)確的，實(shí)際情況并非這樣，差分方程的求解還會(huì)有計(jì)算的誤差。如由于數(shù)字的舍入的誤差而引起的小擾動(dòng)。這類小擾動(dòng)在傳播的過程中會(huì)不會(huì)惡性的增長？以至于“吞沒了差分方程的“真解呢？這就是差分方程的穩(wěn)定問題。在實(shí)際計(jì)算中，某一步產(chǎn)生的擾動(dòng)值在后面的計(jì)算中能夠被控制，甚至是逐步衰減的。 穩(wěn)定性問題比較復(fù)雜，為簡化討論，僅考察下面的模型方程 為保證微分方程本身的穩(wěn)定性，假設(shè) 先考慮

46、Euler方法的穩(wěn)定性。模型方程 的Euler公式為：設(shè)在節(jié)點(diǎn)值 上有一擾動(dòng)值 , 它的傳播使節(jié)點(diǎn)值 產(chǎn)生大小為： 的擾動(dòng)值，假設(shè)用 按Euler 公式得出：yyyyiiyhy)1(1iy1iy1iiiiyy*i0定義：若一種數(shù)值方法在節(jié)點(diǎn)值定義：若一種數(shù)值方法在節(jié)點(diǎn)值 上有擾動(dòng)上有擾動(dòng) ，而對(duì)于，而對(duì)于yi后的后的各個(gè)各個(gè)節(jié)點(diǎn)值節(jié)點(diǎn)值 上產(chǎn)生的上產(chǎn)生的 偏差均不超過偏差均不超過 ,則稱該方法是穩(wěn)則稱該方法是穩(wěn)定的。定的。ijyiiy21,2 ,(10)EulerohoR 這表明方法是條件穩(wěn)定的，其條件為記則穩(wěn)定條件可記為：*1111,18iiiiiyyh 的計(jì)算過程不再有新的誤差，則擾動(dòng)值滿

47、足,可見擾動(dòng)值滿足原來的差分方程（ ） .如果差分方程的解是不增長的，即有：1()8,11,(9)iiyyhh若則算法就是穩(wěn)定的。這一結(jié)論對(duì)后面的研究其他方法的穩(wěn)定性也同樣有效。顯然，為保證差分方程( )的解是不增長的，只要選取的 充分小 使nnnnnyhyyhyy11111xy x若自變量 表征時(shí)間，則 為具有時(shí)間量綱的量，工程中將它稱為時(shí)間常數(shù)。時(shí)間常數(shù) 可用來刻劃原方程的解（ ）的衰減速度，事實(shí)上，在時(shí)間 內(nèi)解)()()(1)1()(oyeeoyeoyyxi12ey xEulerhh衰減了倍。因此。 越小，解（ ）衰減的越快。方法的穩(wěn)定性條件表明，時(shí)間常數(shù) 越小，則穩(wěn)定性對(duì)步長的限制越苛

48、刻。）方法。再考慮后退的歐拉（Euler公式為：其后退對(duì)于模型方程Euleryy,（或稱為無條件穩(wěn)定）方法恒穩(wěn)定。，因而后退的則恒有恒有由于Euleryyhnn111100.0250.0500.075 0.100123450-1-2-3xy其方程時(shí)間常數(shù)為10.01因此有10知，要使Euler法穩(wěn)定則步長02. 02 ch如果取步長h=0.025則Euler格式為：nnyy5 . 11曲線函數(shù)。如圖：是一個(gè)衰減很快的指數(shù)其準(zhǔn)確解為：）（例：考慮初值問題：xeyyyy10010100其結(jié)果在準(zhǔn)確解)(nxy上下波動(dòng)不穩(wěn)定，再看后退Euler格式H=0.025時(shí)，形式為：nnyy5.311計(jì)算結(jié)果

49、穩(wěn)定。具體結(jié)果如下：節(jié)點(diǎn)Eule方法后退Euler方法0.025-1.50.28570.050.2.250.08160.075-3.3750.02330.1005.06250.0067 前面介紹的幾種步進(jìn)方法，在計(jì)算 時(shí)大多只用到前一個(gè)節(jié)點(diǎn)上的近似值 ，而沒有用到前幾步的計(jì)算所得出的信息，故稱為單步法 。實(shí)際上經(jīng)過多次單步法計(jì)算后，已得出一系列近似值 等。為了充分利用這些信息來計(jì)算 以減少計(jì)算工作量和獲得教高的精度，可采用如下計(jì)算公式：1iyiyiyyy,10iyyy,10和1iy11(1)(),.iiy xy是用前面若干節(jié)點(diǎn)處的函數(shù)值與導(dǎo)數(shù)的線性組合來計(jì)算的近似值所以稱為線性多步法0011,

50、0,0.kkk當(dāng)時(shí)為單步法，時(shí)為多步法步法 當(dāng)時(shí)為顯式時(shí)為隱式法）格式即為多步法（二步上面介紹的有關(guān)Euler確定。展開后的待定系數(shù)法來）中的系數(shù)可用（Taylor1),(.,111jijijijjyxff為系數(shù)其中11110kkijijjijjjyyhf （1）1(2) )(,()()(:,x ),( 111iiixxiiidxxyxfxyxyxyxfy區(qū)間上積分有在對(duì)方程由數(shù)值積分法來導(dǎo)出）中的某些特殊公式可（.,),(算方法導(dǎo)出一系列計(jì)但可利用數(shù)值積分公式無法直接算出來因此等式右邊的積分知函數(shù)由于右端積分中含有未xy(2),NewtonfAdams現(xiàn)用向后插值多項(xiàng)式來逼近中從而導(dǎo)出一種線

51、性多步法方法ikiiiikkkiiiifkktttfttftfthxNxNNewtonkkffff!) 1).(1(.! 2) 1()(:).(1.,221由第二章可知向后插值多項(xiàng)式階個(gè)數(shù)據(jù)來構(gòu)造一個(gè)共設(shè)由)(321)(,1)(1110110111iiiiiiiiiiiffhytdtfdttfhyytfftftfNk則時(shí)，當(dāng)1(1)1(1)()()( ),( ,)(1)!kkniiit ttkRxthhfx xk 其余項(xiàng)為：11100(1)2()( )()(),(3)( )1iikikiky xy xhNxth dthRxth dtfxh代入（ ）的：當(dāng)有界且時(shí)則120110()(),(4)(

52、),(5)kkiiikihRxth dto hyyhNxth dt此即為截?cái)嗾`差項(xiàng)。于是可得：012110,( ,)1()2iiiiiiiiikNf yytfyhf x yEulerRh y當(dāng)時(shí)，此即為式格式，局部截?cái)嗾`差為3111011002(2 )115()12,(6)(1),(1)(1)(1)(1)!()iikiikjijjkjjkjmmmjmtmmmkkikiRhyyyhfCtttmCdtdtmRB hy局 部 截 斷 誤 差 為 ：一 般 的 有 ：其 中局 部 截 斷 誤 差 為 ：j012345kB01211232112521223121612583324552459243724

53、972025147201901720277472026167201274720251288955144042771440792314409982144072981440287714404766048010987對(duì)應(yīng)不同k的kj和kB可計(jì)算出來分別列表如下：kkj例如 k=3 時(shí)有：)(720251937595524)5(513211fhRffffhyyiiiiiii （7） （7稱為Adams四步顯式法。它用到了四個(gè)節(jié)點(diǎn)上的f值，是一種最常用的多步法，其精度為四階。如果利用fffkiii11,共k+1個(gè)數(shù)據(jù)來構(gòu)造一個(gè)Newton內(nèi)插多次式 xNk，則與上面類似推導(dǎo)可得：fyyjikjkjiih

54、10*1局部截?cái)嗾`差為：122*1ikkkiyhBR當(dāng)k=0時(shí)，fyyiiih11即為隱式Euler格式。局部截?cái)嗾`差為：1 2121iiyhR當(dāng)k=1時(shí)，ffyyiiiih112即為梯形格式。局部截?cái)嗾`差為：1 31121iiyhR對(duì)于不同k值下的*kj和Bk*可算出，列表如下：012345Bk*012112121121212512812124132492419245241720194720251720646720261720106720191603514404751440142714407981440482144017314402760480863k*kjj11125(5)113(9195

55、(8)2419()7203iiiiiiiihkyyffffRh yk 當(dāng)時(shí),時(shí)的(8)稱為三步四階Admas隱式算法,它只用到前面三個(gè)點(diǎn)的函數(shù)f的值,但卻具有四階精度. 在同一階數(shù)下，隱式的局部截?cái)嗾`差的常數(shù)的絕對(duì)值 比顯式的|Bk|要小。*|kB 顯式的計(jì)算工作量比隱式小。 隱式的穩(wěn)定范圍比顯式大。當(dāng)k=2時(shí)，fffyyiiiiih1118512 1441241iiyhR解：取h=0.1，兩種方法具體算式如下：)9 . 07 . 39 . 55 .18(24/ 1)( 9)(37)(59)(5524/ 1 . 03213211iiiiiiiiiiyyyyyyyyyy四步顯式：1112120.

56、1/249() 19() 5() ()1/24.9(22.10.50.1)iiiiiiiiiyyyyyyyyy 三步隱式：算法啟動(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é)果部分如下表xey例：例： 試分別用Adams四步顯式和三步隱式法求解下列初值問題，并比較兩者所得結(jié)果精度。,0,1 , (0)1dyy xydx 610610610710710710710如下列定解問題：其Adams三步隱式格式為：)5199 (24/2111iiiiyyyyiie

57、eeehyy很難化為yi+1的顯式表達(dá)式，只能用迭代發(fā)求解，這就增加了計(jì)算工作量。在實(shí)際計(jì)算中，往往仿照改進(jìn)的Euler格式的構(gòu)造方程，把顯式和隱式倆種Adams方法結(jié)合起來構(gòu)成預(yù)測一校正系統(tǒng)。1)(25.0 ,0oyxedxdyyAdams隱式方法在計(jì)算過程中,對(duì)于一般的情形往往要解超越方程。以四階Adams方法為例。先由顯式方法算出近似值，作為隱式方法的預(yù)測值，然后再作校正。 11),()9375955(24/1113211iiiiiiiiiyxffffffhyy預(yù)測：1112111/24 9195(,)iiiiiiiiiyyhffffff xy校正：（）上式計(jì)算yi+1既要用到它前一步的

58、信息yi和fi,還要用到更前面三步的信息fi-1 , fi-2 , fi-3 。它是一種一步法，無法自動(dòng)啟動(dòng)，需要用其他四階單步法如四階R-K法先從y0求出y1， y2， y3作為初值。然后按(11)進(jìn)行迭代。例：用Admas預(yù)測校正方法求解下列初值問題：20,1(0)1dyxyxdxyy解：取h=0.1，按11進(jìn)行計(jì)算。結(jié)果如下表。表中 和 分別為預(yù)測值、校正值和精確解iy ,iy)(ixy( )1yxx 0110.11.09541.09540.21.18321.18320.31.26491.26490.41.34151.34161.34160.51.41411.41421.41420.61.48321.48321.48320.71.54911.54921.54920.81.61241.61241.61250.91.67331.67331.673311.7321.7321.731ixiyiy)(ixyAdams預(yù)測一校正方法與Richardson外推法相結(jié)合，可以提高計(jì)算精度而不會(huì)增加過大的計(jì)算工作量。預(yù)測公式的局部截?cái)嗾`差為：校正公式的局部截?cái)嗾`差為：5(5)61115(5)251()()720251()(12)