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

1、編輯ppt1 第九章第九章 常微分方程的數(shù)值解法常微分方程的數(shù)值解法 1、引言引言 2、初值問題的數(shù)值解法、初值問題的數(shù)值解法-單步法單步法 3、龍格、龍格-庫塔方法庫塔方法 4、收斂性與穩(wěn)定性收斂性與穩(wěn)定性 5、初值問題的數(shù)值解法、初值問題的數(shù)值解法多步法多步法 6、方程組和剛性方程、方程組和剛性方程 7、習(xí)題和、習(xí)題和總結(jié)總結(jié)主主 要要 內(nèi)內(nèi) 容容編輯ppt2主主 要要 內(nèi)內(nèi) 容容研究的問題研究的問題數(shù)值解法的意義數(shù)值解法的意義1、 引引 言言編輯ppt3現(xiàn)實世界中大多數(shù)事物內(nèi)部聯(lián)系非常復(fù)雜找出其狀態(tài)和狀態(tài)變化規(guī)律之間的相互聯(lián)系，也即一個或一些函數(shù)與他們的導(dǎo)數(shù)之間的關(guān)系此種關(guān)系的數(shù)學(xué)表達(dá)就

2、為1.1.什么是微分方程什么是微分方程 ? ?編輯ppt4如何建立數(shù)學(xué)模型已在建模課程中得到討論，各類微分方程本身和他們的解所具有的特性已在常微分方程及數(shù)學(xué)物理方程中得以解釋，編輯ppt5尋找解析解的過程稱為求解微分方程組。尋找解析解的過程稱為求解微分方程組。 一個或一組具有所要求階連續(xù)導(dǎo)數(shù)的解析函數(shù)，將它代入微分方程（組），恰使其所有條件都得到滿足的解稱為解析解解析解( (或古典解或古典解),),稱為真解或解。稱為真解或解。編輯ppt64.4.什么是微分方程的數(shù)值解什么是微分方程的數(shù)值解? ?雖然求解微分方程有許多解析方法,但解析方法只能夠求解一些特殊類型的方程,從實際意義上來講我們更關(guān)心的

3、是某些 特定的自變量在某一個定義范圍內(nèi)的一系列離散點一系列離散點上的近似值.尋找數(shù)值解的過程稱為數(shù)值求解微分方程。尋找數(shù)值解的過程稱為數(shù)值求解微分方程。數(shù)值解數(shù)值解編輯ppt7在大量的實際方程中出現(xiàn)的函數(shù)起碼的連續(xù)性都無法保證，更何況要求階的導(dǎo)數(shù)求解數(shù)值解很多微分方程很多微分方程根本求不到根本求不到問題的解析解！問題的解析解！重要手段。編輯ppt8常微分方程的數(shù)值解法常用來求近似解根據(jù)提供的算法通過計算機通過計算機便捷地實現(xiàn)便捷地實現(xiàn)5.常微分方程數(shù)值解法的特點常微分方程數(shù)值解法的特點編輯ppt900( , )(1.1)( ) (1.2)yf x yaxby xy 本章主要討論一階常微方程的初

4、值問題6.基本知識基本知識其中f (x,y)是已知函數(shù)，(1.2)是定解條件也稱為初值條件。初值條件。編輯ppt10則稱 f (x,y) 對y 滿足李普希茲條件，L 稱為Lipschitz常數(shù).常微分方程的理論指出：常微分方程的理論指出：當(dāng) f (x,y) 定義在區(qū)域 G=(axb,y)若存在正的常數(shù) L 使：1212|( ,)( ,)| (1.3)f x yf x yL yy12 , ,xa byy使得對任意的及都成立就可保證方程解的存在唯一性就可保證方程解的存在唯一性(Lipschitz)條件條件編輯ppt11我們以下的討論，都在滿足上述條件下進行.若 f (x,y) 在區(qū)域 G連續(xù)，關(guān)于

5、y一階常微分方程的初值問題的解存在且唯一問題的解存在且唯一.一階常微分方程組常表述為：1111( ,) ( ) ( ,)mmmmyf x yyaxbyfx yy 1010() ()mmyxyx編輯ppt12寫成向量形式為(1)()00000( , )(),(,)mTyf x yaxby xyxxx高階常微分方程定解問題如二階定解問題：(,) ()()yfxyyaxbyay a編輯ppt13這些解法都可以寫成向量形式編輯ppt14簡單的數(shù)值方法與基本概念簡單的數(shù)值方法與基本概念 1. 簡單簡單歐拉法（歐拉法（Euler） 2后退的歐拉法后退的歐拉法 3梯形法梯形法 4改進改進EulerEuler

6、法法 2、初值問題的數(shù)值解法、初值問題的數(shù)值解法單步法單步法編輯ppt151. 簡單的歐拉簡單的歐拉(Euler)方法方法考慮模型：00( , ) (1.1)() (1.2)yf x yaxby xy 在精度要求不高時通過歐拉方法的討論歐拉方法歐拉方法編輯ppt162. 歐拉方法的導(dǎo)出歐拉方法的導(dǎo)出把區(qū)間a,b分為n個小區(qū)間步長為要計算出解函數(shù) y(x) 在一系列節(jié)點()iiyy x)/()( iiixaihhhban,一般取即等距節(jié)點處的近似值01 naxxxb1(-)iiihxxN N等分等分編輯ppt1700( ,)(1.1)() (1.2)yfx yaxby xy 對微分方程(1.1)

7、兩端從1nnxx到進行積分11( ,( )nnnnxxxxy dxfx y xdx11()()( , ( )nnxnnxy xy xf x y x dx編輯ppt18右端積分用右端積分用左矩形數(shù)值左矩形數(shù)值求積公式求積公式2()( )()()()2bagg x dxba g aba( )( , ( )g xf x y x令1)1(, ()(, () nnnnxxnnf xy xnnyyf xy xh得編輯ppt19x0 x11()()()(,)nnnnnny xy xhy xyh f xy)1,., 0(),(1 niyxfhyyiiii亦稱為亦稱為歐拉折線法歐拉折線法 /* Eulers p

8、olygonal arc method*/ 每步計算只用到1ny或用向前差商近似導(dǎo)數(shù)1()()()nnny xy xyxh100021111(,)(,) (,)nnnnyyhf xyyyhf xyyyhf xy依上述公式逐次計算可得：ny例題編輯ppt203.歐拉公式有明顯的幾何意義依此類推得到一折線故也稱Euler為單步法。ny公式右端只含有已知項所以又稱為顯格式的單步法。000011122(,) ( )(,)( ) (,)(,) xyy xxyy xxxxyxy 過 點 的 曲 線 是 解 在 作 的 切 線 與 直 線 交 于 再 作 切 線 交 于 編輯ppt21也稱歐拉折線法.就是用

9、這條折線近似地代替曲線( )y x( )yy xxynx1nxnp1np1np x編輯ppt22從上述幾何意義上得知，由Euler法所得的折線明顯偏離了積分曲線，可見此方法非常粗糙。非常粗糙。4.歐拉法的局部截斷誤差：在假設(shè)第 i 步計算是精確的前提下，考慮截斷誤差稱為局部截斷誤差/* local truncation error */。1()iiiRyxy編輯ppt23定義定義若某算法的局部截斷誤差為O(hp+1)，則稱該算法有p 階精度。Ri 的的主項主項/* leading term */ 歐拉法的局部截斷誤差：232()()hiyxO h歐拉法具有 1 階精度。),()()()()()

10、(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR 定義定義在假設(shè) yi = y(xi)，即第 i 步計算是精確的前提下，考慮的截斷誤差 Ri = y(xi+1) yi+1 稱為局部截斷誤差 /* local truncation error */。編輯ppt24如果單步差分公式的局部截斷誤差為O(hp+1),則稱該公式為p p階方法階方法.這里p為非負(fù)整數(shù).顯然,階數(shù)越高,方法的精度越高.Taylor展開式,一元函數(shù)的Taylor展開式為: 321! 3)(! 2)()()()()(hxyhxyhxyxyhxyxynnnnnn23() ( )( )( )() ( ,)112

11、hRy xyy xhy xy xO hyhf x yiiiiiiiii232()()hiyxO h若某算法的局部截斷誤差為若某算法的局部截斷誤差為O(hp+1)，則稱該算法有，則稱該算法有p 階精度。階精度。Ri 的的主項主項/* leading term */編輯ppt255. 歐拉公式的改進歐拉公式的改進： 隱式歐拉法隱式歐拉法 /* implicit Euler method */向后差商近似導(dǎo)數(shù)向后差商近似導(dǎo)數(shù)hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,., 0(),(111 niyxfhyyiiii由于未知數(shù)由于未知數(shù) yi+1 同時出現(xiàn)在

12、等式的兩邊，不能直接得到，故同時出現(xiàn)在等式的兩邊，不能直接得到，故稱為稱為隱式隱式 /* implicit */ 歐拉公式，而前者稱為歐拉公式，而前者稱為顯式顯式 /* explicit */ 歐拉公式。歐拉公式。一般先用顯式計算一個初值，再一般先用顯式計算一個初值，再迭代迭代求解。求解。 隱式隱式歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：11)(iiiyxyR)()(322hOxyih 即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。編輯ppt266.梯形公式梯形公式 /* trapezoid formula */ 顯、隱式兩種算法的顯、隱式兩種算法的平均平均)1,., 0(),

13、(),(2111 niyxfyxfhyyiiiiii注：注：的確有局部截斷誤差的確有局部截斷誤差 ， 即梯形公式具有即梯形公式具有2 階精度，比歐拉方法有了進步。階精度，比歐拉方法有了進步。但注意到該公式是但注意到該公式是隱式隱式公式，計算時不得不用到公式，計算時不得不用到迭代法，其迭代收斂性與歐拉公式相似。迭代法，其迭代收斂性與歐拉公式相似。)()(311hOyxyRiii 中點歐拉公式中點歐拉公式 /* midpoint formula */中心差商近似導(dǎo)數(shù)中心差商近似導(dǎo)數(shù)hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(21

14、1 niyxfhyyiiii假設(shè)假設(shè) ，則可以導(dǎo)出，則可以導(dǎo)出即中點公式具有即中點公式具有 2 階精度。階精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 需要需要2個初值個初值 y0和和 y1來啟動遞推來啟動遞推過程，這樣的算法稱為過程，這樣的算法稱為雙步法雙步法 /* double-step method */，而前面的三種算法都是，而前面的三種算法都是單步法單步法 /* single-step method */。編輯ppt27方方 法法 顯式歐拉顯式歐拉隱式歐拉隱式歐拉梯形公式梯形公式中點公式中點公式簡單簡單精度低精度低穩(wěn)定性最好穩(wěn)定性最好精度低精度低, 計算

15、量大計算量大精度提高精度提高計算量大計算量大精度提高精度提高, 顯式顯式多一個初值多一個初值, 可能影響精度可能影響精度編輯ppt28 改進歐拉法改進歐拉法 /* modified Eulers method */Step 1: 先用先用顯式顯式歐拉公式作歐拉公式作預(yù)測預(yù)測，算出，算出),(1iiiiyxfhyy Step 2: 再將再將 代入代入隱式隱式梯形公式的右邊作梯形公式的右邊作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy注：注：此法亦稱為此法亦稱為預(yù)測預(yù)測-校正法校正法 /* predictor-corrector method */。可以證明該算

16、法具有可以證明該算法具有 2 階精度，同時可以看到它是個階精度，同時可以看到它是個單單步步遞推格式，比隱式公式的迭代求解過程遞推格式，比隱式公式的迭代求解過程簡單簡單。后面將。后面將看到，它的看到，它的穩(wěn)定性高穩(wěn)定性高于顯式歐拉法。于顯式歐拉法。 )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii編輯ppt293 龍格龍格 - 庫塔法庫塔法 /* Runge-Kutta Method */ 考察改進的歐拉法，可以將其改寫為：考察改進的歐拉法，可以將其改寫為：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1

17、K2 的的平均值平均值嗎？嗎？步長一定是一個步長一定是一個h 嗎？嗎？單步遞推法的單步遞推法的基本思想基本思想是從是從 ( xi , yi ) 點出發(fā)，以點出發(fā)，以某一斜某一斜率率沿直線達(dá)到沿直線達(dá)到 ( xi+1 , yi+1 ) 點。歐拉法及其各種變形所點。歐拉法及其各種變形所能達(dá)到的最高精度為能達(dá)到的最高精度為2階階。建立高精度的單步遞推格式。建立高精度的單步遞推格式。編輯ppt30首先希望能確定系數(shù)首先希望能確定系數(shù) 1、 2、p，使得到的算法格式有，使得到的算法格式有2階階精度，即在精度，即在 的前提假設(shè)下，使得的前提假設(shè)下，使得 )(iixyy )()(311hOyxyRiii S

18、tep 1: 將將 K2 在在 ( xi , yi ) 點作點作 Taylor 展開展開)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 將改進歐拉法推廣為：將改進歐拉法推廣為：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 將將 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxy

19、xyhyyiiiiiiii 編輯ppt31Step 3: 將將 yi+1 與與 y( xi+1 ) 在在 xi 點的點的泰勒泰勒展開作比較展開作比較)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，則必須有：，則必須有：)()(311hOyxyRiii21,1221 p 這里有這里有 個未知個未知數(shù)，數(shù)， 個方程。個方程。32存在存在無窮多個解無窮多個解。所有滿足上式的格式統(tǒng)稱為。所有滿足上式的格式統(tǒng)稱為2階龍格階龍格 - 庫庫塔格式塔格式。21, 121 p注意到，注意到， 就是改進的歐拉法。就是改進的歐拉法

20、。 Q: 為獲得更高的精度，應(yīng)該如何進一步推廣？為獲得更高的精度，應(yīng)該如何進一步推廣？編輯ppt32其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均為待定均為待定系數(shù)，確定這些系數(shù)的系數(shù)，確定這些系數(shù)的步驟與前面相似。步驟與前面相似。 ).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用為四級最常用為四級4階階經(jīng)典龍格經(jīng)典龍格-庫塔法庫塔法 /

21、* Classical Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 編輯ppt33注：注： 龍格龍格-庫塔法庫塔法的主要運算在于計算的主要運算在于計算 Ki 的值，即計算的值，即計算 f 的的值。值。Butcher 于于1965年給出了計算量與可達(dá)到的最高精年給出了計算量與可達(dá)到的最高精度階數(shù)的關(guān)系：度階數(shù)的關(guān)系：753可達(dá)到的最高精度可達(dá)到的最高精度642每步須算每步須算Ki 的個數(shù)的個數(shù))(2hO)(3hO)(4hO)(5hO)(6hO)

22、(4hO)(2nhO8n 由于龍格由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好解函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用采用低階算法低階算法而將步長而將步長h 取小取小。深入研究龍格深入研究龍格-庫塔法請看此處！庫塔法請看此處！編輯ppt344 收斂性與穩(wěn)定性收斂性與穩(wěn)定性 /* Convergency and Stability */ 1.收斂性收斂性 /* Convergency */定義定義 若某算法對于任意固定的若某算法對于任意固定的 x = xi = x0 + i h，當(dāng)，當(dāng) h0 ( 同

23、時同時 i ) 時有時有 yi y( xi )，則稱該算法是，則稱該算法是收斂收斂的。的。 例：例：就初值問題就初值問題 考察歐拉顯式格式的收斂性?？疾鞖W拉顯式格式的收斂性。 0)0(yyyy 解：解：該問題的精確解為該問題的精確解為 xeyxy 0)( 歐拉公式為歐拉公式為iiiiyhyhyy)1 (1 0)1 (yhyii 對任意固定的對任意固定的 x = xi = i h ，有，有iixhhxihyhyy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 編輯ppt35 2.穩(wěn)定性穩(wěn)定性 /* Stability */例：例：考察初值問題考察初值問題 在區(qū)間在區(qū)間

24、0, 0.5上的解。上的解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精確解精確解改進歐拉法改進歐拉法 歐拉隱式歐拉隱式歐拉顯式歐拉顯式 節(jié)點節(jié)點 xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.765

25、6 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong ?!編輯ppt36定義定義若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計算中都算中都逐步衰減逐步衰減，則稱該算法是，則稱該算法是絕對穩(wěn)定的絕對穩(wěn)定的 /*absolutely stable */。一般分析時為簡單起見，只考慮一般分析時為簡單起見，只考慮試驗方程試驗方程 /* test equation */yy 常數(shù)，可以常數(shù)，可以是復(fù)數(shù)是復(fù)數(shù)當(dāng)步長取為當(dāng)步長取為 h 時，將某算法應(yīng)用于上式

26、，并假設(shè)只在初值時，將某算法應(yīng)用于上式，并假設(shè)只在初值產(chǎn)生誤差產(chǎn)生誤差 ，則若此誤差以后逐步衰減，就稱該，則若此誤差以后逐步衰減，就稱該算法相對于算法相對于 絕對穩(wěn)定絕對穩(wěn)定， 的全體構(gòu)成的全體構(gòu)成絕對穩(wěn)定區(qū)域絕對穩(wěn)定區(qū)域。我們稱我們稱算法算法A 比算法比算法B 穩(wěn)定穩(wěn)定，就是指，就是指 A 的絕對穩(wěn)定區(qū)域比的絕對穩(wěn)定區(qū)域比 B 的的大大。000yy h h h編輯ppt37例：例：考察顯式歐拉法考察顯式歐拉法011)1(yhyhyyiiii 000yy 011)1(yhyii 01111)1( iiiihyy由此可見，要保證初始誤差由此可見，要保證初始誤差 0 以后逐步衰減，以后逐步衰減，必

27、須滿足：必須滿足：hh 1|1| h0-1-2ReImg例：例：考察隱式歐拉法考察隱式歐拉法11 iiiyhyy iiyhy 11101111 iih可見絕對穩(wěn)定區(qū)域為：可見絕對穩(wěn)定區(qū)域為：1|1| h210ReImg注：注：一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式法的好。法的好。編輯ppt38例：例：隱式龍格隱式龍格-庫塔法庫塔法 ),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii 而而顯式顯式 1 4 階方法的絕對穩(wěn)定階方法的絕對穩(wěn)定區(qū)域為區(qū)域為 )2,2(1111KhyhxfKhKyyiiii其中其中

28、2階方法階方法 的絕對穩(wěn)定區(qū)域為的絕對穩(wěn)定區(qū)域為0ReImgk=1k=2k=3k=4-1-2-3-123ReImg無條件穩(wěn)定無條件穩(wěn)定編輯ppt39例例1 用歐拉方法用歐拉方法，隱式歐拉方法隱式歐拉方法和和歐拉中點公式歐拉中點公式計算計算01(0 )1yyxy 的近似解，取步長的近似解，取步長h=0.1，并與精確值比較，并與精確值比較解解:歐拉方法的算式為：歐拉方法的算式為：10(1)0,1, 9 .1iiyyhiy歐拉隱式方法在本題可解出方程，不必迭代，公式為：歐拉隱式方法在本題可解出方程，不必迭代，公式為：100,1, 9 .11iiyyihy歐拉中點公式是兩步法，第一步歐拉中點公式是兩步

29、法，第一步y(tǒng)1用歐拉公式，以后用公式用歐拉公式，以后用公式11021, 2 , 9 .1iiiyyh yiy編輯ppt40本題精確解為本題精確解為y=e-x,計算結(jié)果見表計算結(jié)果見表9-1例例2 用用歐拉公式歐拉公式和和梯形公式梯形公式的預(yù)估校正法計算：的預(yù)估校正法計算：的數(shù)值解，取的數(shù)值解，取h=0.1，梯形公式只迭代一次，并與精確值比較，梯形公式只迭代一次，并與精確值比較.方程的解析解為方程的解析解為:201(0 )1xyyxyy xy21 解解: 本例中歐拉公式為：本例中歐拉公式為：1020.20.11.10,1,9.1iiiiiiiixxyyyyiyyy編輯ppt41梯形公式只校正一次

30、的格式為梯形公式只校正一次的格式為(0)1(0)111(0)100.21.1220.051iiiiiiiiiiiixyyyxxyyyyyyy結(jié)果列入表結(jié)果列入表9.2編輯ppt421. Runge-Kutta 1. Runge-Kutta 法的一般形式法的一般形式 2. 22. 2階階Runge-Kutta Runge-Kutta 方法方法 3. 3. 經(jīng)典經(jīng)典Runge-Kutta Runge-Kutta 方法方法 4 4R-K-Fehhlberg R-K-Fehhlberg 方法方法 5. 5. 隱式隱式R-KR-K方法方法 6. 6. 變步長方法變步長方法 編輯ppt431.Runge-

31、1.Runge-Kutta Kutta 法的一般形式法的一般形式Runge-Kutta 法是用區(qū)間上若干個點上的導(dǎo)數(shù)的線性組合得到平均斜率，以提高方法的階。其一般形式為 ：11Lnniiiyyhk1222133331132211,1,2,(,)(,)(,()nnnnnniininii jjjkfxc h yc ha kiLkf x ykf xc h yc hkkf xc h yc h a ka k編輯ppt44式(9.11) 稱L級級p階階Runge-Kutta方法方法(簡稱R-K法法)。當(dāng)L1就是歐拉法，當(dāng)L2時為改進的歐拉法。 1111 , 1 , 1 Liiii jijca。111()(

32、)LnnniiiTy xy xhk其中它的局部截斷誤差是(9.11)2. 22. 2級級2 2階階Runge-KuttaRunge-Kutta方法方法令令 L=2L=2，則，則 11 122()nnyyhkk編輯ppt4512221(,)(,)nnnnkfxykfxc h yc hk111 122()()()nnnTy xy xhkk其局部截斷誤差是其局部截斷誤差是： 這是3個未知量的兩個方程，其中有一個自由參數(shù)，方程組有無窮多解，從而得到一族2級2階R-KR-K方法 。1222112c編輯ppt463. 3. 經(jīng)典經(jīng)典Runge-KuttaRunge-Kutta方法方法我們可以構(gòu)造出一族3級

33、3階，一族4級4階和一族5級4階等R-K方法。最常用的4級4階是如下的經(jīng)典經(jīng)典R-KR-K方法方法：112341213243(22)6(,)(,)22(,)22(,)nnnnnnnnnnhyykkkkkfxyhhkfxykhhkfxykkfxh yhk編輯ppt474 4R-K-Fehhlberg R-K-Fehhlberg 方法方法R-K-FehhlbergR-K-Fehhlberg方法方法是在R-K方法的基礎(chǔ)上引進誤差和步長控制的辦法。即利用5階R-K法 112456166656285619213512825564305055iiyyKKKKK估計4階R-K的局部誤差，其中 5431151

34、410421972565104821625KKKKyyii)329323,83()41,4(),(213121KKyhxhfKKyhxhfKyxhfKiiiiii編輯ppt48)401141041859256553442278,2()410485451336808216439,()219772962197720021971932,1312(543216432153214KKKKKyhxhfKKKKKyhxhfKKKKyhxhfKiiiiii注：注：R-K-Fehhlberg比4階R-K方法具有更大的優(yōu)越性，他是計算穩(wěn)定，高精度的方法，他的顯著優(yōu)點是，每一步僅需計算f的6個值；若用4階R-K-L

35、與5階R-K-L在一起使用，則每步需要計算f的10個值。推薦使用！ 編輯ppt495. 5. 隱式隱式R-KR-K方法方法類似于顯式R-K公式，稍加改變，就得到隱式隱式R-K方法方法。11Lnniiiyyhk1,Lininii jjjkfxc h yc ha k它與顯式R-K公式的區(qū)別在于：顯式公式中對系數(shù)求和的上限是i-1，從而構(gòu)成的矩陣是一個嚴(yán)格下三角陣。而在隱式公式中對系數(shù)求和的上限是L，從而構(gòu)成的矩陣是方陣，需要用迭代法求出Ki，。推導(dǎo)隱式公式的思路和方法與顯式R-K法類似。通常，同級的隱式公式獲得比顯式公式更高的階。編輯ppt50通常，同級的隱式公式獲得比顯式公式更高的階。常用的隱式

36、通常，同級的隱式公式獲得比顯式公式更高的階。常用的隱式R-KR-K法有法有:2,2111nnnnnyyhxhfyy),(),(2111nnnnnnyxfyxfhyy112112212()2331323,6412333231,6124nnnnnnhyykkkfxh yk hk hkfxh yk hk h1級級2階中點公式階中點公式 :2級級2階梯形公式：階梯形公式： 2級級4階階R-K公式：公式： 編輯ppt516.6.變步長方法變步長方法在單步法中每一積分步步長實際上是相互獨立的，步長的選擇具有了靈活性。根據(jù)合理地選擇每一積分步的步長，既保證精度的要求，又可以減少計算量，從而減小舍入誤差。其方

37、便的控制手段是基于誤差的事后估計式。對于給定的精度 ,如果 ,反復(fù)將步長折半進行計算，直至 為止,這時取最終得到的作為結(jié)果；如果 為止，這時再將步長折半一次，就得到所要的結(jié)果。這種通過加倍或折半處理步長的計算方法稱為變步長方法。 注注：推薦使用精度好計算量低的變步長方法。用四階顯式R-K方法做變步長方法是實踐中較好的方法！ 編輯ppt52 例 分別用改進的歐拉格式和四階龍格庫塔格式解初值問題（取步長h=0.2）：2(0)1xyyyy (01)x編輯ppt53節(jié)點 改進歐拉法 四階龍格庫塔法 準(zhǔn)確解 0 1 1 1 0.2 1.186667 1.183229 1.1832160.4 1.3483

38、12 1.341667 1.3416410.6 1.493704 1.483281 1.4832400.8 1.627861 1.612514 1.612452 1 1.754205 1.732142 1.732051表（注：已指出過準(zhǔn)確解12yx）編輯ppt54單步法的相容性、收斂性和穩(wěn)定性單步法的相容性、收斂性和穩(wěn)定性 1.1.單步法的相容性單步法的相容性 2.2.單步法的收斂性單步法的收斂性 3.3.單步法的穩(wěn)定性單步法的穩(wěn)定性 4.4.相容性和收斂性的關(guān)系相容性和收斂性的關(guān)系 5.5.相容性和方法階的關(guān)系相容性和方法階的關(guān)系 6.6.穩(wěn)定性和收斂性穩(wěn)定性和收斂性 7.7.絕對穩(wěn)定性和絕

39、對穩(wěn)定域絕對穩(wěn)定性和絕對穩(wěn)定域 編輯ppt55定義一定義一：對于（：對于（9.1.1）常微分方程初值問題）常微分方程初值問題單步法的形式可以變表示為單步法的形式可以變表示為 （9.2.19）其中其中 h為步長為步長若對求解區(qū)間中任一固定的若對求解區(qū)間中任一固定的x，當(dāng)，當(dāng) 時皆成立時皆成立 100( ,; )nnnnyyh x y hyy0h00()()limlim( ,;)hhy xhy xx y hh則稱由（9.2.19）確定的單步法與微分方程初值問題是相容相容的注注意到上式左邊極限為由（由（9.1.19.1.1）知它應(yīng)等于從而由相容性定義得稱相容性條件相容性條件。 ( , ;0)( ,

40、)x yf x y編輯ppt56單步法的收斂性單步法的收斂性定義二：定義二： 設(shè) y(x) 是（9.1.1）的解， 是單步法（9.2.19）產(chǎn)生的數(shù)值解，對于每一個固定的 , ， 當(dāng) 即 。若成立 ， 則稱該方法是收斂的。nynx0h1nnxx1()nnyy x1nnxxh編輯ppt57定義三定義三： 若一個數(shù)值方法在基點 處的值有 的擾動，在 此后各基點 （mn）處的值產(chǎn)生的偏差均不超 過 ，則稱該方法是的。 單步法的穩(wěn)定性有以下定理nymy定理二定理二： 若單步法的增量函數(shù)對變量y滿足 Lipschtiz 條件， 則單步方法是的。編輯ppt58相容性和收斂性的關(guān)系相容性和收斂性的關(guān)系 若單

41、步法的增量函數(shù)對變量y滿足Lipschtiz 條件，即存在與 h , x 無關(guān)的常數(shù) L，對區(qū)域D= 任意兩點（x，y1）,(x,y2)成立，則單步法收斂的充分必要條件是相容性條件成立。（讀者自證）,axby 編輯ppt59相容性和方法階的關(guān)系相容性和方法階的關(guān)系 若單步法是若單步法是p p階方法則成立階方法則成立 若單步法滿足相容性條件，得若單步法滿足相容性條件，得 所以所以 =0=0也就是說單步法的階數(shù)一定要是正數(shù)。由于我們考慮也就是說單步法的階數(shù)一定要是正數(shù)。由于我們考慮的單步法皆為正整數(shù)，的單步法皆為正整數(shù)，p p至少為一。因此我們考慮的至少為一。因此我們考慮的單步法都滿足相容性條件。

42、單步法都滿足相容性條件。 1()( )( , ; )()py x hy xhx y hO h101( )( , )()limphy xf x yO hh0()limphO h編輯ppt60 穩(wěn)定性和收斂性的關(guān)系穩(wěn)定性和收斂性的關(guān)系若單步法的增量函數(shù)滿足定理二的條件即單若單步法的增量函數(shù)滿足定理二的條件即單步法是穩(wěn)定的則單步法收斂的充分必要條件步法是穩(wěn)定的則單步法收斂的充分必要條件是是相容性條件成立。相容性條件成立。 ( , ;0)( , )x yf x y編輯ppt61絕對穩(wěn)定性和絕對穩(wěn)定域絕對穩(wěn)定性和絕對穩(wěn)定域 穩(wěn)定性問題是一個比較復(fù)雜的問題。為了簡化討論一穩(wěn)定性問題是一個比較復(fù)雜的問題。為

43、了簡化討論一般僅對試驗方程般僅對試驗方程 進行考察。這里假定進行考察。這里假定Re0,Re0h0和的允許范圍，稱為該方法的和的允許范圍，稱為該方法的 絕對穩(wěn)定域絕對穩(wěn)定域。 yynn ky編輯ppt62絕對穩(wěn)定性和絕對穩(wěn)定域絕對穩(wěn)定性和絕對穩(wěn)定域2 2將Euler方法應(yīng)用到試驗方程得誤差方程是要求誤差不增長則nnnnyhhyyy)1 (1nnnh11|1|1hnn編輯ppt63則Euler 方法的絕對穩(wěn)定域是以 為半徑的圓的內(nèi)部。為了保證穩(wěn)定性步長有所控制。假如當(dāng) 時h應(yīng)滿足 ，當(dāng) 時， h 應(yīng)滿足 等等。同樣我們可以將試驗方程用到其它各種單步法當(dāng)中，求出其絕對穩(wěn)定域。在實際應(yīng)用中必須在這個范

44、圍內(nèi)，否則誤差傳播相當(dāng)大，即不穩(wěn)定。 1h5520 h10050110020 h編輯ppt64 1.Adams1.Adams方法方法 2.2.米爾尼方法、漢明方法及辛普森方法米爾尼方法、漢明方法及辛普森方法 3.3.預(yù)測校正方法預(yù)測校正方法 4.4.多步法的相容性、穩(wěn)定性和收斂性多步法的相容性、穩(wěn)定性和收斂性 5 初值問題的數(shù)值解法初值問題的數(shù)值解法多步法多步法編輯ppt65 考慮型如考慮型如 的的k步法，步法， 稱為稱為阿當(dāng)姆斯阿當(dāng)姆斯（Adams）方法）方法10knknkiniiyyhf112330,0,00101211211211211224133121443211.Adams1.Ada

45、ms方法方法拿其中一種為例，推導(dǎo)其公式。對上面所說的法一般形拿其中一種為例，推導(dǎo)其公式。對上面所說的法一般形式式若取若取 ，有，有方程組方程組同樣解之，得到同樣解之，得到3步步4階隱式階隱式Adams公式和它的余項。公式和它的余項。編輯ppt661112(9195)24nnnnnnhyyffff5(5)12119(),(,)720nnnnnTh yxx 其它請讀者自證。我們僅將結(jié)論列于下表。其它請讀者自證。我們僅將結(jié)論列于下表。 Adams顯式公式顯式公式1pc1nnnyyhf12211( 3)2nnnnhyyff5123221(23165)12nnnnnhyyfff3843321(55593

46、79)24nnnnnnhyyffff251720kp公式11223344編輯ppt67Adams隱式公式隱式公式1pc11()2nnnnhyyff1122121(58)12nnnnnhyyfff12432321(9195)24nnnnnnhyyffff19720434321(25164626410619 )720nnnnnnnhyyfffff3160kp公式12233445注：注：在這些方法中，在這些方法中，4階的階的Adams預(yù)測校正方法具有方法預(yù)測校正方法具有方法簡潔、使用方便、易排序、高精度等優(yōu)點。尤其當(dāng)函數(shù)簡潔、使用方便、易排序、高精度等優(yōu)點。尤其當(dāng)函數(shù)f比較復(fù)雜時更能顯出它的優(yōu)越性。

47、比較復(fù)雜時更能顯出它的優(yōu)越性。編輯ppt68 同 理 得 到同 理 得 到 5 個 待 定 參 數(shù) 方 程 組 。 解 之個 待 定 參 數(shù) 方 程 組 。 解 之得得 ， ， ， ， 。 構(gòu)成著名的構(gòu)成著名的Miline 4Miline 4步步4 4階顯式公式和它的余階顯式公式和它的余項。項。 13380341， 3820313nnyy124(22)3nnnhfff5(5)13114( ),(,)45nnnnnTh yxx2.2.米爾尼方法、漢明方法及米爾尼方法、漢明方法及辛普森方法辛普森方法編輯ppt69 同理得到同理得到5 5個參數(shù)方程組。求解后就構(gòu)成著名的個參數(shù)方程組。求解后就構(gòu)成著名

48、的3 3步步4 4階隱式階隱式HammingHamming公式和它的余項。公式和它的余項。 若取，也得到若取，也得到5個參數(shù)方程組。求解后就構(gòu)成個參數(shù)方程組。求解后就構(gòu)成Simpson隱式公式和其余項。隱式公式和其余項。 121(9)8nnnyyy113(2)8nnnhfff5(5)1211(),(,)40nnnnnTh yxx米爾尼方法、漢明方法及米爾尼方法、漢明方法及辛普森方法辛普森方法1111(4)3nnnnnhyyfff5(5)1111(),(,)90nnnnnTh yxx 編輯ppt70 不論單步法或多步法，隱式公式比顯式公式穩(wěn)定性好，但在實際使用隱式公式時，都會遇到兩個問題：一個是

49、隱式公式如何能方便地進行計算；另一個是實際計算步長取多大。如隱式梯形公式，每往前推進一步，不必進行多次迭代，而是采用一階顯式Euler公式預(yù)測，二階隱式梯形公式校正一次，構(gòu)成顯式改進Euler公式，能達(dá)到與梯形公式同階的精度，即二階精度。 3. 預(yù)測校正方法預(yù)測校正方法編輯ppt71 對于線性多步公式，一般地，不難驗證，如果 預(yù)測公式是階或階精度，校正公式為階精度， 則用預(yù)測公式提供初值，校正公式迭代一次的 效果也能達(dá)到階精度，再迭代下去，效果就不 明顯了。預(yù)測-校正技術(shù)即保證了計算精度，又 使隱式計算顯式化，克服了隱式公式要反復(fù)迭 代的困難。至于實際計算步長的選取，我們對 預(yù)測-校正公式使用

50、外推原理，得到誤差估計式 ，用來調(diào)整計算步長，使達(dá)到控制誤差要求。對于線性多步方法常用的預(yù)測-校正公式有Miline-Hamming方法和4階階Adams顯隱式顯隱式預(yù)測-校正公式編輯ppt72 將Miline Miline 公式公式和HammingHamming公式公式結(jié)合，構(gòu)成預(yù)測-校正公式 Miline公式公式和Hamming公式公式的局部截斷誤差分別為13124(22)3pnnnnnhyyfff1211113(9)( (,) 2)88pnnnnnnnhyyyf xyff5(5)611114()()()45MMnnnnTy xyh yxO h5(5)61111()( )()40HHnnn

51、nTy xyh yxO h修正修正Miline-Hamming公式公式編輯ppt73利用外推原理，即加上后消去局部截斷誤差的主項。使 說明經(jīng)過外推后的算法其精度提高一階。忽略誤差項，上式可改寫為 11611144045()()1144045MHnnnyyy xO h111114114() ()040454045MHnnny xyy11191129112() ()0360360360360MHnnny xyy編輯ppt74 由于 和是 在計算過程中獲得的數(shù)據(jù)，稱為Miline Miline 公式公式和HammingHamming公式公式的事后誤差估計式。我們可以用它們來調(diào)節(jié)計算步長的大小，即選擇

52、一個合適的的步長，使 ，其中的是要求達(dá)到的計算精度。 111121 ()91120MHnnny xyy1111112()()121MHMnnnny xyyy11119()()121HHMnnnny xyyy1Mny1Hny119()121HMnnyy編輯ppt75 又可得到 MilineMiline公式公式和HammingHamming公式公式的修正公式，它們分別是 從而構(gòu)成如下的修正從而構(gòu)成如下的修正HammingHamming預(yù)測預(yù)測- -校正公式，校正公式，簡稱修正簡稱修正HammingHamming公式：公式： 1111112()121MmMHMnnnnyyyy11119()121Hm

53、HHMnnnnyyyy編輯ppt76預(yù)測：修正：校正： 修正： 在應(yīng)用這套公式時，先由同階單步法提供初值 ， ， ， 。計算 時可取。31npnyy)22(3421nnnfffh)(12111211pncnpnpmnyyyy)2),(83)9 (8111121nnpmnnnncnffyxfhyyy)(12191111pncncnnyyyy0y1y2y3y4ypcyy33編輯ppt77 將將4步步4階顯式階顯式Adams公式作為預(yù)測公式和公式作為預(yù)測公式和3步步4階階式式Adams公式作為校正公式，構(gòu)成公式作為校正公式，構(gòu)成Adams預(yù)測預(yù)測-校正公式。校正公式。 它們的局部截斷誤差分別是它們的

54、局部截斷誤差分別是 1123(5559379)24pnnnnnnhyyffff1112(9195)24cnnnnnnhyyffff5(5)6111251()()()720ppnnnnTy xyh yxO h5(5)611119()( )()720ccnnnnTy xyh yxO hAdams預(yù)測預(yù)測-校正公式校正公式編輯ppt78利用外推原理，將上兩式作線性組合消去局部截斷誤差的主項。使計算精度至少提高一階，同時得到兩個修正公式，將它們和上兩式構(gòu)成了如下修正修正Adams公式公式：預(yù)測： 修正： 校正: 1123(5559379)24pnnnnnnhyyffff11251()270pmpcpn

55、nnnyyyy11112(9 (,) 195)24cpmnnnnnnnhyyf xyfff編輯ppt79修正： 同理，在計算時，調(diào)節(jié)計算步長 使 , 由同階單步法提供初值 , ， ， 。計算 時，可取 。 理論分析得出它們的絕對穩(wěn)定區(qū)域，修正Hamming法： ; 4階Adams預(yù)測預(yù)測-校正法校正法： 其中 )(270191111pncncnnyyyyh)(2701911pncnyy0y1y2y3y4ypcyy33069. 0h075.0h)(,(xyxyf我們也可以在教學(xué)演示上看出修正的預(yù)測校正格式的誤差非常的小。編輯ppt80 多步法的相容性條件多步法的相容性條件 k步法的一般形式為 其

56、中 由對單步法的討論可知,當(dāng) 時,方法階數(shù)至少為1。即對 y1,y=x 應(yīng)精確成立。令 y 1，所以令y=x得 101111011 ()nnnrn rnnnrn ryyyyhffff0| , 000k0 ,y100 ,kk1y4.4.多步法的相容性、穩(wěn)定性和收斂性多步法的相容性、穩(wěn)定性和收斂性編輯ppt81所以所以我們稱為我們稱為線性多步法的相容性條件。線性多步法的相容性條件。 1010101110() (1) ()(1)nkkkkkkkkkkxh kkhkk 即編輯ppt82 k步法需要k個出發(fā)值，而初值問題只提供了一個初值，在使用k步法時尚需要其它方法作補充k-1個出發(fā)值，今假定它們 是

57、，當(dāng) 這k-1個值都 應(yīng)收斂于共同的極限y（a）即 在討論多步法收斂性時我們總假定（9.3.12）成立 定義四定義四： 1,.2, 1),(klhyll0h 1,.2 , 1,)(limlim00klayhylhlhnknkknknknkknkfffhyayaya011011.編輯ppt83 的解 收斂于初值問題的解 y(x)。這里 xa+nh. 定義五定義五：稱多項式 為k步法的特征多項式。如果特征多項式的零點皆不大于1，且等 于1的零點是單重的，稱根條件成立。稱多項式 為第二特征多項式。 顯然根條件可以表示 定理三定理三：k步法收斂的充要條件為： （1）相容性條件成立。 （2）特征多項式的

58、零點皆不大于1，且等于1的零點是單重的。 （稱2為）特征根條件。 ny,baxnx011.)(kkkk011.)(kkkk) 1 () 1 (, 0) 1 (編輯ppt84多步法的穩(wěn)定性 關(guān)于線性多步法成立以下定理：關(guān)于線性多步法成立以下定理： 定理四：定理四：若函數(shù)若函數(shù)f(x,y)f(x,y)對變量對變量y滿足滿足Lipschtiz 條件在條件在與與h,x無關(guān)的常數(shù)無關(guān)的常數(shù)L，對區(qū)域?qū)^(qū)域D= D= 任意兩點（任意兩點（x x，y1y1）,(x,y2),(x,y2)成立成立 k步法的相容性、收斂性、穩(wěn)定性有以下關(guān)系步法的相容性、收斂性、穩(wěn)定性有以下關(guān)系 對于常微分方程右端函數(shù)對于常微分方

59、程右端函數(shù)f(x,y)，在相容性條件成立情，在相容性條件成立情況下，況下，k步法的收斂性和穩(wěn)定性是等價的。步法的收斂性和穩(wěn)定性是等價的。 事實上在相容性條件成立時，收斂性和穩(wěn)定性的充要條事實上在相容性條件成立時，收斂性和穩(wěn)定性的充要條件都是特征根條件成立。件都是特征根條件成立。ybxaD,| 12| ); 2,(); 1,(|yyLhyxhyx編輯ppt85多步法的絕對穩(wěn)定性和絕對穩(wěn)定域多步法的絕對穩(wěn)定性和絕對穩(wěn)定域 將線性多步法的公式應(yīng)用到試驗方程 進行考察。這里假定Re 0,即試驗方程本身是穩(wěn)定的。 記 得 是常系數(shù)差分方程，其特征方程為 記它的 k個特征根為 并設(shè)它們是互異的。顯然根與

60、有關(guān)，不妨記為 注意到當(dāng) 時 這時由特征方程得 由線性多步法的相容性條件得 是一個根。不妨設(shè), 差分方程的解為 其中系數(shù)由線性多步法的出發(fā)值確定。yyhu inkiikiinkyy000)()(kii,.2, 1, hu kii,.2 , 1),(0h00)(10, 0)(1nkknnndddy.2211編輯ppt86另一方面，y(0)=1的試驗方程的精確解為 ， 設(shè)多步法截斷誤差為 ，由此可得 ，我們稱 為主根，其它根都為增根。 定義五定義五：線性多步法的絕對穩(wěn)定區(qū)域?qū)o定的 ，如果特征方程的特征根 皆按模小于1，則線性多步法關(guān)于u是絕對穩(wěn)定的。使得 成立的 構(gòu)成絕對穩(wěn)定區(qū)域。注：從誤差角度