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

1、第第7 7章章 常微分方程數(shù)值解法常微分方程數(shù)值解法 7.1 7.1 引言引言 一、有關(guān)常微分方程一、有關(guān)常微分方程 二、數(shù)值解法二、數(shù)值解法三、數(shù)值解法的三種類型三、數(shù)值解法的三種類型 /* Numerical Methods for Ordinary Differential Equations */一、有關(guān)常微分方程一、有關(guān)常微分方程1. 1. 什么是常微分方程的初值問題？什么是常微分方程的初值問題？| ),(),(|yyLyxfyxf 00)(),(yxyyxfy理論上可以證明：只要函數(shù)理論上可以證明：只要函數(shù)f(xf(x，y)y)適當(dāng)光適當(dāng)光滑滑關(guān)于關(guān)于y y滿足李普希茲滿足李普希茲

2、(Lipschitz)(Lipschitz)條件條件常微分方程常微分方程初值問題初值問題則初值問題的解存在唯一。則初值問題的解存在唯一。例例0)0(2yxy0)0(21yxyy思考：常微分方程中的未知數(shù)是什么？思考：常微分方程中的未知數(shù)是什么？2.2.常微分方程的一般解常微分方程的一般解( (解析解解析解) ) 對(duì)一些典型的微分方程對(duì)一些典型的微分方程( (可分離變量方可分離變量方程，一階線性方程等等程，一階線性方程等等) )，有可能找出它們，有可能找出它們的一般解表達(dá)式，然后用初始條件確定表的一般解表達(dá)式，然后用初始條件確定表達(dá)式中的任意常數(shù)，這樣解即能確定。達(dá)式中的任意常數(shù)，這樣解即能確定

3、。 例如例如 求解求解0)0(2yxy解：分離變量得解：分離變量得 dy=2xdx 積分得積分得y=x2+c 由初值得由初值得c=0 故解為故解為y=x2 注：生產(chǎn)和科研中所處理的微分方注：生產(chǎn)和科研中所處理的微分方程往往很復(fù)雜且大多得不出一般解程往往很復(fù)雜且大多得不出一般解 二、常微分方程的數(shù)值解二、常微分方程的數(shù)值解 由于在實(shí)用上對(duì)初值問題，一由于在實(shí)用上對(duì)初值問題，一般是要求得到解在若干點(diǎn)上滿足規(guī)般是要求得到解在若干點(diǎn)上滿足規(guī)定精確度的近似值定精確度的近似值yi，或者是得到，或者是得到一個(gè)滿足精確度要求的便于計(jì)算的一個(gè)滿足精確度要求的便于計(jì)算的近似表達(dá)式。近似表達(dá)式。 故常微分方程的數(shù)值

4、解就是求故常微分方程的數(shù)值解就是求出在若干點(diǎn)上解的近似值。出在若干點(diǎn)上解的近似值。 定義：常微分方程初值問題的數(shù)值解定義：常微分方程初值問題的數(shù)值解 一般是指在由初始點(diǎn)一般是指在由初始點(diǎn)x0開始的若干開始的若干離散的離散的x值處，即對(duì)值處，即對(duì)x0 x1x2xn，求出準(zhǔn)確值求出準(zhǔn)確值y(x1),y(x2),，y(xn)的近似的近似值值y1，y2，yn 本章只討論本章只討論x0, x1,xn等距的情等距的情況，設(shè)況，設(shè) xi+1- -xi = h, i = 0, 1, , n-1上式中的上式中的h值稱為步長(zhǎng)值稱為步長(zhǎng) .對(duì)于常微分初值問題對(duì)于常微分初值問題一般解一般解 y = y(x)數(shù)值解數(shù)值

5、解x x0 x1 x2 xn y y0 y1 y2 yn 注：本章使用的符號(hào)注：本章使用的符號(hào)y(xi)：一般解：一般解y=y(x)在在x=xi 處的精確值處的精確值.yi： 一 般 解： 一 般 解 y = y ( x ) 在在 x = xi處處 的近似值的近似值.三、數(shù)值解法的三種類型三、數(shù)值解法的三種類型1.用差商代替導(dǎo)數(shù)用差商代替導(dǎo)數(shù) 由導(dǎo)數(shù)定義由導(dǎo)數(shù)定義 hxyhxyxyn)()(lim)( 0故若故若h的值較小，則有的值較小，則有 hxyhxyxy)()()( 代入代入y= f (x,y)可得可得 ),()()(yxfhxyhxy即即 y(x+h)y(x)+hf(x,y) 故原初值

6、問題可離散化為故原初值問題可離散化為 2 , 1 , 0)(),(001ixyyyxhfyyiiii 于是由初始值于是由初始值y(x0)=y0出發(fā)，可出發(fā)，可依次地計(jì)算出依次地計(jì)算出 y1=y0+hf(x0,y0) y2=y1+hf(x0+h,y1) yn=yn-1+hf(x0+(n-1)h,yn-1) 2.使用泰勒公式使用泰勒公式 在微分方程在微分方程y=f(x,y)中，中，y是是x及及y(x)的函數(shù)的函數(shù).由于精確值由于精確值y(x+h)在在h=0處的泰勒展式為處的泰勒展式為 )(24)(6)(2)( )()()4(432xyhxyhxyhxhyxyhxy根據(jù)根據(jù)y y=f(x,y)=f(

7、x,y)，得公式得公式若取泰勒展式的前兩項(xiàng)，則有若取泰勒展式的前兩項(xiàng)，則有y(x+h) y(x)+hy(x)2 , 1 , 0)(),(001ixyyyxhfyyiiiiyffxfdxdyyfxfy又因?yàn)橛忠驗(yàn)閥=f(x,y)，故，故 故若取泰勒展式的前三項(xiàng)，則可故若取泰勒展式的前三項(xiàng)，則可得公式得公式 nnnnyhhyyy221 ),(),(),(2),(2nnynnnnxnnnyxfyxfyxfhyxhfy 與上類似，一般可取公式為如下形與上類似，一般可取公式為如下形式式 )(21! 2pnpnnnnyphyhhyyy 注：應(yīng)用泰勒公式求數(shù)值解，從形式上注：應(yīng)用泰勒公式求數(shù)值解，從形式上看

8、簡(jiǎn)單，其實(shí)具體構(gòu)造這種公式往往是相當(dāng)看簡(jiǎn)單，其實(shí)具體構(gòu)造這種公式往往是相當(dāng)困難的，因?yàn)樗枰峁?dǎo)數(shù)值，困難的，因?yàn)樗枰峁?dǎo)數(shù)值，y y( (j j) )n n當(dāng)階當(dāng)階數(shù)提高時(shí)，求導(dǎo)過程可能很復(fù)雜，因此泰勒數(shù)提高時(shí)，求導(dǎo)過程可能很復(fù)雜，因此泰勒公式通常不直接使用，但可以用它來啟發(fā)思公式通常不直接使用，但可以用它來啟發(fā)思路。路。 3.使用數(shù)值積分的方法使用數(shù)值積分的方法對(duì)于微分方程對(duì)于微分方程 y=f(x,y)兩邊求兩邊求x0到到x的定積分的定積分 dxxyxfxxdxdxdyxx)(,(00 即即 dttytfxxxyxy)(,()()(00 或?qū)憺榛驅(qū)憺?dttytfxxyxy)(,()

9、(00 這就是與初值問題這就是與初值問題 00)(),(yxyyxfy等價(jià)的積分方程。只要用某種數(shù)值等價(jià)的積分方程。只要用某種數(shù)值積分方法便可建立起近似公式。積分方法便可建立起近似公式。 例：對(duì)積分部分應(yīng)用左矩形公式，則例：對(duì)積分部分應(yīng)用左矩形公式，則有有 dttytfxxxyxyiiii)(,()()(11)(,(21iiiixyxfxx 2 , 1 , 0)(),(001ixyyyxhfyyiiii例：對(duì)積分部分應(yīng)用梯形公式，則有例：對(duì)積分部分應(yīng)用梯形公式，則有 dttytfxxxyxyiiii)(,()()(11 )(,()(,(2111 iiiiiixyxfxyxfxx故其數(shù)值公式為故

10、其數(shù)值公式為 ),(),(2111iiiiiiyxfyxfhyy其特點(diǎn)是其特點(diǎn)是y yi i+1+1表示成表示成x xi i，x xi i+1+1及及y yi i的的隱函數(shù)隱函數(shù). . 7.2 簡(jiǎn)單數(shù)值方法與基本概念簡(jiǎn)單數(shù)值方法與基本概念一、歐拉公式一、歐拉公式二、隱式的歐拉公式二、隱式的歐拉公式 三、梯形公式三、梯形公式四、中點(diǎn)歐拉公式四、中點(diǎn)歐拉公式 五、改進(jìn)的歐拉公式五、改進(jìn)的歐拉公式 考慮一階常微分方程的初值問題考慮一階常微分方程的初值問題 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy要計(jì)算出解函數(shù)要計(jì)算出解函數(shù) y(x) 在一系列

11、節(jié)點(diǎn)在一系列節(jié)點(diǎn) a = x0 x1 xn= b 處的近似值處的近似值),., 1()(nixyyii 節(jié)點(diǎn)間距節(jié)點(diǎn)間距 為步長(zhǎng)，通常采用等距節(jié)點(diǎn)，為步長(zhǎng)，通常采用等距節(jié)點(diǎn)，即取即取 hi = h (常數(shù)常數(shù))。) 1,., 0(1 nixxhiii7.2.1 7.2.1 歐拉方法歐拉方法 歐拉公式：歐拉公式：x0 x1向前差商近似導(dǎo)數(shù)向前差商近似導(dǎo)數(shù)hxyxyxy)()()(010 ),()()()(000001yxfhyxyhxyxy 1y 記為記為) 1,., 0(),(1 niyxfhyyiiii亦稱為歐拉折線法亦稱為歐拉折線法 定義定義在假設(shè)在假設(shè) yi = y(xi)，即第，即第

12、i 步計(jì)算是精確的步計(jì)算是精確的前提下，考慮的截?cái)嗾`差前提下，考慮的截?cái)嗾`差 Ri = y(xi+1) yi+1 稱為稱為局部截?cái)嗾`差局部截?cái)嗾`差 /* local truncation error */。定義定義若某算法的局部截?cái)嗾`差為若某算法的局部截?cái)嗾`差為O(hp+1)，則稱該算法有則稱該算法有p 階精度。階精度。 歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：),()()()()()(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR )()(322hOxyih 歐拉法具有歐拉法具有 1 階精度。階精度。7.2.2 7.2.2 隱式歐拉法隱式歐拉法 /* implicit

13、 Euler method */向后差商近似導(dǎo)數(shù)向后差商近似導(dǎo)數(shù)hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,.,0(),(111 niyxfhyyiiii由于未知數(shù)由于未知數(shù) yi+1 同時(shí)出現(xiàn)在等式的兩邊，不能直同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱為隱式接得到，故稱為隱式 /* implicit */ 歐拉公式，而歐拉公式，而前者稱為顯式前者稱為顯式 /* explicit */ 歐拉公式。歐拉公式。)1,.,0(),(111 niyxfhyyiiii 隱式隱式歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：11)(iiiyxyR)()(322hO

14、xyih 即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。7.2.3 7.2.3 梯形公式梯形公式 /* trapezoid formula */ 顯、隱式兩種算法的平均顯、隱式兩種算法的平均)1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：的確有局部截?cái)嗾`差注：的確有局部截?cái)嗾`差 ， 即梯形公式具有即梯形公式具有2 階精度，比歐拉方法有了進(jìn)步。但階精度，比歐拉方法有了進(jìn)步。但注意到該公式是隱式公式，計(jì)算時(shí)不得不用到迭代法，注意到該公式是隱式公式，計(jì)算時(shí)不得不用到迭代法，其迭代收斂性與歐拉公式相似。其迭代收斂性與歐拉公式相似。)()(311hOyxyRiii

15、7.2.47.2.4中點(diǎn)歐拉公式中點(diǎn)歐拉公式 /* midpoint formula */中心差商近似導(dǎo)數(shù)中心差商近似導(dǎo)數(shù)hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假設(shè)假設(shè) ，則可以導(dǎo)出，則可以導(dǎo)出即中點(diǎn)公式具有即中點(diǎn)公式具有 2 階精度。階精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 方方 法法 各方法的比較：各方法的比較：7.2.57.2.5 改進(jìn)歐拉法改進(jìn)歐拉法 /* modified Eulers method */Step 1: 先用顯式歐拉公式作預(yù)測(cè)

16、，算出先用顯式歐拉公式作預(yù)測(cè)，算出),(1iiiiyxfhyy Step 2: 再將再將 代入隱式梯形公式的右邊作校代入隱式梯形公式的右邊作校正，得到正，得到1 iy),(),(2111 iiiiiiyxfyxfhyy )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii7.3 龍格龍格 - 庫塔法庫塔法 /* Runge-Kutta Method */一、再看一、再看Taylor方法方法二、二、Runge-Kutta 方法方法 0)(,),(yaybaxyxfdxdy復(fù)習(xí)復(fù)習(xí) 一元函數(shù)的一元函數(shù)的Taylor展開展開 )(24)(6)(2)( )()()4(4

17、32xyhxyhxyhxhyxyhxy 20000000200000000000)(,()()(,(2)(,(21)(,()(,(),(),(yyyxfyyxxyxfxxyxfyyyxfxxyxfyxfyxfyxyxyx 二元函數(shù)的二元函數(shù)的Taylor展開展開一、再看一、再看Taylor方法方法 )(24)(6)(2)( )()()4(432xyhxyhxyhxhyxyhxy一般可取公式為如下形式 )(21! 2pnpnnnnyphyhhyyy 注：應(yīng)用泰勒展式構(gòu)造公式時(shí)，希望注：應(yīng)用泰勒展式構(gòu)造公式時(shí)，希望近近似值的泰勒展式與精確值的泰勒展式有盡可似值的泰勒展式與精確值的泰勒展式有盡可能多

18、的項(xiàng)相同能多的項(xiàng)相同，即：兩者關(guān)于，即：兩者關(guān)于h h的展式中的的展式中的系數(shù)系數(shù)盡可能多的項(xiàng)相同盡可能多的項(xiàng)相同，這樣局部截?cái)嗾`差，這樣局部截?cái)嗾`差與與h h的更高次同階，從而獲得高精度的近似的更高次同階，從而獲得高精度的近似公式。公式。 建立高精度的單步遞推格式。建立高精度的單步遞推格式。單步遞推法的基本思想是從單步遞推法的基本思想是從 ( xi , yi ) 點(diǎn)出點(diǎn)出發(fā)，以某一斜率沿直線達(dá)到發(fā)，以某一斜率沿直線達(dá)到 ( xi+1 , yi+1 ) 點(diǎn)。歐拉法及其各種變形所能達(dá)到的最點(diǎn)。歐拉法及其各種變形所能達(dá)到的最高精度為高精度為2階階。二、二、Runge-Kutta 方法方法 考察改進(jìn)

19、的歐拉法，可以將其改寫為：考察改進(jìn)的歐拉法，可以將其改寫為：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 首先希望能確定系數(shù)首先希望能確定系數(shù) 1、 2、p，使得到的算法格式有，使得到的算法格式有2階階精度，即在精度，即在 的前提假設(shè)下，使得的前提假設(shè)下，使得 )(iixyy )()(311hOyxyRiii Step 1: 將將 K2 在在 ( xi , yi ) 點(diǎn)作點(diǎn)作 Taylor 展開展開(用二元函數(shù)泰勒展開用二元函數(shù)泰勒展開)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOx

20、yphxyii 將改進(jìn)歐拉法推廣為：將改進(jìn)歐拉法推廣為：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 將將 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii Step 3: 將將 yi+1 與與 y( xi+1 ) 在在 xi 點(diǎn)的點(diǎn)的泰勒泰勒展開作比較展開作比較)()()()(322211hOxyphxyhyyiiii )()(2)()()(

21、321hOxyhxyhxyxyiiii 要求要求 ，則必須有：，則必須有：)()(311hOyxyRiii21,1221 p 這里有這里有 個(gè)未知個(gè)未知數(shù)，數(shù)， 個(gè)方程。個(gè)方程。32存在無窮多個(gè)解。所有滿足上式的格式統(tǒng)稱為存在無窮多個(gè)解。所有滿足上式的格式統(tǒng)稱為2階龍格階龍格 - 庫庫塔格式。塔格式。21, 121 p注意到，注意到， 就是改進(jìn)的歐拉法。就是改進(jìn)的歐拉法。 問題問題: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, ,

22、 i 1 ) 均為待定系均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似。數(shù)，確定這些系數(shù)的步驟與前面相似。 ).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用為四級(jí)最常用為四級(jí)4階經(jīng)典龍格階經(jīng)典龍格-庫塔法庫塔法 /* Classical Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 注：注： 龍格

23、龍格-庫塔法庫塔法的主要運(yùn)算在于計(jì)算的主要運(yùn)算在于計(jì)算 Ki 的值，即的值，即計(jì)算計(jì)算 f 的值。的值。Butcher 于于1965年給出了計(jì)算量年給出了計(jì)算量與可達(dá)到的最高精度階數(shù)的關(guān)系：與可達(dá)到的最高精度階數(shù)的關(guān)系：753可達(dá)到的最高可達(dá)到的最高精度精度642每步須算每步須算Ki 的的個(gè)數(shù)個(gè)數(shù))(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n 由于龍格由于龍格-庫塔法的導(dǎo)出基于泰勒展開，庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對(duì)故精度主要受解函數(shù)的光滑性影響。對(duì)于光滑性不太好的解，最好采用低階算于光滑性不太好的解，最好采用低階算法而將步長(zhǎng)法而將步

24、長(zhǎng)h 取小。取小。7.4 收斂性與穩(wěn)定性 /*Convergency and Stability */一、收斂性一、收斂性二、穩(wěn)定性二、穩(wěn)定性定義定義 若某算法對(duì)于任意固定的若某算法對(duì)于任意固定的 x = xi = x0 + i h，當(dāng)，當(dāng) h0 時(shí)有時(shí)有 yi y( xi )，則稱該算法是收斂的。，則稱該算法是收斂的。 例：就初值問題例：就初值問題 考察歐拉顯式格式的收斂性?？疾鞖W拉顯式格式的收斂性。 0)0(yyyy 解：該問題的精確解為解：該問題的精確解為 xeyxy 0)( 歐拉公式為歐拉公式為iiiiyhyhyy)1 (1 0)1 (yhyii 對(duì)任意固定的對(duì)任意固定的 x = xi

25、 = i h ，有，有iixhhxihyhyy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 一、收斂性一、收斂性 /* Convergency */ 二、二、穩(wěn)定性穩(wěn)定性 /* Stability */例：考察初值問題例：考察初值問題 在區(qū)間在區(qū)間0, 0.5上的解。上的解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精確解精確解改進(jìn)歐拉法改進(jìn)歐拉法 歐拉隱式歐拉隱式歐拉顯式歐拉顯式 節(jié)點(diǎn)節(jié)點(diǎn) xixey30 1.0000 2.0000 4.

26、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.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong ?!定義定義若某算法在計(jì)算過程中任一步產(chǎn)生的誤差在以后的計(jì)若某算法在計(jì)算過程中任一步產(chǎn)生的誤差在以后的計(jì)算中都逐步衰減，則稱該算法是絕對(duì)穩(wěn)定的算中都逐步衰減，則稱

27、該算法是絕對(duì)穩(wěn)定的 /*absolutely stable */。一般分析時(shí)為簡(jiǎn)單起見，只考慮試驗(yàn)方程一般分析時(shí)為簡(jiǎn)單起見，只考慮試驗(yàn)方程 /* test equation */yy 當(dāng)步長(zhǎng)取為當(dāng)步長(zhǎng)取為 h 時(shí)，將某算法應(yīng)用于上式，并假設(shè)只在初值時(shí)，將某算法應(yīng)用于上式，并假設(shè)只在初值產(chǎn)生擾動(dòng)產(chǎn)生擾動(dòng) ，則若此誤差以后逐步衰減，就稱該，則若此誤差以后逐步衰減，就稱該算法相對(duì)于算法相對(duì)于 絕對(duì)穩(wěn)定，絕對(duì)穩(wěn)定， 的全體構(gòu)成絕對(duì)穩(wěn)定區(qū)域。的全體構(gòu)成絕對(duì)穩(wěn)定區(qū)域。我們稱算法我們稱算法A 比算法比算法B 穩(wěn)定，就是指穩(wěn)定，就是指 A 的絕對(duì)穩(wěn)定區(qū)域比的絕對(duì)穩(wěn)定區(qū)域比 B 的大。的大。000yy h h

28、h例：考察顯式歐拉法例：考察顯式歐拉法011)1(yhyhyyiiii 000yy 011)1(yhyii 01111)1( iiiihyy由此可見，要保證初始誤差由此可見，要保證初始誤差 0 以后逐步衰減，以后逐步衰減，必須滿足：必須滿足：hh 1|1| h0-1-2ReImg例：考察隱式歐拉法例：考察隱式歐拉法11 iiiyhyy iiyhy 11101111 iih可見絕對(duì)穩(wěn)定區(qū)域?yàn)椋嚎梢娊^對(duì)穩(wěn)定區(qū)域?yàn)椋?|1| h210ReImg注：一般來說，隱式歐拉法的絕對(duì)穩(wěn)定性比同階的顯式注：一般來說，隱式歐拉法的絕對(duì)穩(wěn)定性比同階的顯式法的好。法的好。例：隱式龍格例：隱式龍格-庫塔法庫塔法 ),.

29、, 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii 而顯式而顯式 1 4 階方法的絕對(duì)穩(wěn)定階方法的絕對(duì)穩(wěn)定區(qū)域?yàn)閰^(qū)域?yàn)?)2,2(1111KhyhxfKhKyyiiii其中其中2階方法階方法 的絕對(duì)穩(wěn)定區(qū)域?yàn)榈慕^對(duì)穩(wěn)定區(qū)域?yàn)?ReImgk=1k=2k=3k=4-1-2-3-123ReImg無條件穩(wěn)定無條件穩(wěn)定舉例：對(duì)穩(wěn)定性的理解舉例：對(duì)穩(wěn)定性的理解考察初值問題考察初值問題 1)0(100yyy其準(zhǔn)確解是其準(zhǔn)確解是xexy100)( 這是一個(gè)按指數(shù)曲線衰減得很快的函數(shù)。這是一個(gè)按指數(shù)曲線衰減得很快的函數(shù)。用用EulerEuler方法解方程，得方法解方程，得nny

30、hy)1001(1 (1)(1)若取若取h=0.025h=0.025，則，則EulerEuler公式的具體形式為公式的具體形式為nnnyyhy5 . 1)1001(1 (2)(2)若取若取h=0.005h=0.005，則，則EulerEuler公式的具體形式為公式的具體形式為nnnyyhy5 . 0)1001(1 (3)(3)若取若取h=0.025h=0.025，則后退的，則后退的EulerEuler法公式為法公式為nnyy5 . 311 將三種情況的計(jì)算結(jié)果列表如下將三種情況的計(jì)算結(jié)果列表如下節(jié)點(diǎn)節(jié)點(diǎn)EulerEuler方法方法(h=0.025)(h=0.025)EulerEuler方法方法

31、(h=0.005)(h=0.005)后退后退EulerEuler法法0.0250.025-1.5-1.50.50.50.28570.28570.0500.0502.252.250.250.250.08160.08160.0750.075-3.375-3.3750.1250.1250.02330.02330.1000.1005.06255.06250.06250.06250.00670.0067從表中可以得出結(jié)論：從表中可以得出結(jié)論：不穩(wěn)定不穩(wěn)定穩(wěn)定穩(wěn)定穩(wěn)定穩(wěn)定考察初值問題考察初值問題 1)0(100yyy舉例舉例分別用分別用EulerEuler法和后退法和后退EulerEuler法求解，從法求

32、解，從穩(wěn)定性角度出發(fā)，穩(wěn)定性角度出發(fā)，h h應(yīng)如何取值？應(yīng)如何取值？解解EulerEuler方法的穩(wěn)定域?yàn)榉椒ǖ姆€(wěn)定域?yàn)?|1| h即即1|1001| h解得解得02. 00 h則若取則若取h=0.025h=0.025，計(jì)算過程不穩(wěn)定。，計(jì)算過程不穩(wěn)定。 若取若取h=0.005h=0.005，計(jì)算過程穩(wěn)定。，計(jì)算過程穩(wěn)定。后退后退EulerEuler方法的穩(wěn)定域?yàn)榉椒ǖ姆€(wěn)定域?yàn)?|1| h由于這里的由于這里的= -100 0= -100 0，則不論，則不論h h取何值，取何值，計(jì)算過程都是穩(wěn)定的。計(jì)算過程都是穩(wěn)定的。7.5 線性多步法線性多步法初值問題初值問題y = f (x, y) y (x

33、0) = y0 與下面的積分方程等價(jià)：與下面的積分方程等價(jià)：, 2, 1, 0)(,()()(11 ndttytfxyxynnxxnn 1)(1nnkxnndttPyy（1）求出開頭幾個(gè)點(diǎn)上的近似值，即計(jì)算）求出開頭幾個(gè)點(diǎn)上的近似值，即計(jì)算“表頭表頭”；線性多步法：線性多步法：（2）利用利用 逐步求后面點(diǎn)逐步求后面點(diǎn)xk上的值上的值yk。 1)(1nnkxnndttPyy1. 阿當(dāng)姆斯外推公式阿當(dāng)姆斯外推公式 以以xn-2，xn-1，xn為節(jié)點(diǎn)作牛頓向后插值多項(xiàng)式為節(jié)點(diǎn)作牛頓向后插值多項(xiàng)式P2(x)。nnnfttftfxP22! 2) 1()(其中其中 hxxtnnnnnnnnnnfffhyd

34、tfttftfhyy2102112521! 2) 1(插值公式的余項(xiàng)為插值公式的余項(xiàng)為 nnxxyttthxr 232)(! 3)2)(1()(則積分公式的截?cái)嗾`差為則積分公式的截?cái)嗾`差為 )(0)(! 3)2)(1(31031hdtyttthRn k = 3時(shí)的外推公式為時(shí)的外推公式為nnnnnnffffhyy3218312521余項(xiàng)為余項(xiàng)為 ：)(0)(7202515)5(51hyhRn將差分表示成函數(shù)值的和的形式：將差分表示成函數(shù)值的和的形式：10) 1(nijjijnjfcf二階阿當(dāng)姆斯外推公式可改寫為：二階阿當(dāng)姆斯外推公式可改寫為： )51623(12211nnnnnfffhyy三階

35、阿當(dāng)姆斯外推公式可改寫為：三階阿當(dāng)姆斯外推公式可改寫為： )9375955(243211nnnnnnffffhyy1!) 11() 1(0jijcijjjc例例 求y = x + y, y (0) = 1當(dāng)x = 0.1到0.5，步長(zhǎng)h = 0.1時(shí)的數(shù)值解。解解 先用前面講過的方法計(jì)算出表頭y0 = 1； y1 = 1.11034； y2 = 1.24281； y3 = 1.39972將上述值代入公式，計(jì)算得：y4 = 1.58364； y5 = 1.79742)9375955(243211nnnnnnffffhyy2.阿當(dāng)姆斯內(nèi)插公式阿當(dāng)姆斯內(nèi)插公式將被積函數(shù)將被積函數(shù)用用以以xn-1，x

36、n，xn+1為插值節(jié)點(diǎn)的內(nèi)插多項(xiàng)式為插值節(jié)點(diǎn)的內(nèi)插多項(xiàng)式 12112! 2) 1()(nnnfttftfxP得到：得到：13111112112121nnnnnnffffhyy)(041hRn)85(12111nnnnnfffhyyk = 2 時(shí)時(shí)131211124112121nnnnnnffffhyy)5199(242111nnnnnnffffhyy)(051hRnk=1為了提高精度，經(jīng)常把為了提高精度，經(jīng)常把與與聯(lián)合起來交替使用。例如聯(lián)合起來交替使用。例如 )51623(1221)0(1 nnnnnfffhyy21) 0 (111519,924nnnnnnnfffyxfhyy 第一個(gè)方程的精

37、度為第一個(gè)方程的精度為0(h4)，用第二個(gè)方程迭代一，用第二個(gè)方程迭代一次精度仍達(dá)次精度仍達(dá)0(h5)。這樣兩個(gè)方程交替使用可達(dá)較好。這樣兩個(gè)方程交替使用可達(dá)較好的效果。第一個(gè)方程是的效果。第一個(gè)方程是，第二個(gè)方程是，第二個(gè)方程是。解解 先用前面講過的方法計(jì)算出表頭y0 = 1； y1 = 1.11034； y2 = 1.24281再用下面的第一個(gè)式子計(jì)算出 ，最后用第二個(gè)式子進(jìn)行迭代，得y3 = 1.39972同樣，可算出y4及y5。 )51623(1221)0(1 nnnnnfffhyy21) 0 (111519,924nnnnnnnfffyxfhyy39964. 1)0(3y解二階常微分方程邊值問題的差分法解二階常微分方程邊值問題的差分法考慮常微分方程的邊值問題：考慮常微分方程的邊值問題： )(,)();()()(byaybxaxfyxqyxpy其中其中p(x)，q(x)和和f (x)均為均為a, b上給定的函數(shù)上給定的函數(shù)， ， 為已