第六章36微分方程數(shù)值解法

1、 第6章 常微分方程數(shù)值解法 微分方程數(shù)值解一般可分為：微分方程數(shù)值解一般可分為：常微分方程數(shù)值解常微分方程數(shù)值解和和偏微分偏微分方程數(shù)值解方程數(shù)值解。自然界與工程技術(shù)中的許多現(xiàn)象，其數(shù)學(xué)表達(dá)式。自然界與工程技術(shù)中的許多現(xiàn)象，其數(shù)學(xué)表達(dá)式可歸結(jié)為可歸結(jié)為常微分方程（組）的定解問(wèn)題常微分方程（組）的定解問(wèn)題。一些。一些偏微分方程偏微分方程問(wèn)題問(wèn)題也可以也可以轉(zhuǎn)化轉(zhuǎn)化為為常微分方程常微分方程問(wèn)題來(lái)（近似）求解。問(wèn)題來(lái)（近似）求解。Newton最早采最早采用數(shù)學(xué)方法研究二體問(wèn)題，其中需要求解的運(yùn)動(dòng)方程就是常微用數(shù)學(xué)方法研究二體問(wèn)題，其中需要求解的運(yùn)動(dòng)方程就是常微分方程。許多著名的數(shù)學(xué)家，如分方程。許

2、多著名的數(shù)學(xué)家，如 Bernoulli（家族），（家族），Euler、Gauss、Lagrange和和Laplace等，都遵循歷史傳統(tǒng)，研究重要等，都遵循歷史傳統(tǒng)，研究重要的力學(xué)問(wèn)題的數(shù)學(xué)模型，在這些問(wèn)題中，許多是常微分方程的的力學(xué)問(wèn)題的數(shù)學(xué)模型，在這些問(wèn)題中，許多是常微分方程的求解。作為科學(xué)史上的一段佳話，求解。作為科學(xué)史上的一段佳話，海王星海王星的發(fā)現(xiàn)就是通過(guò)對(duì)常的發(fā)現(xiàn)就是通過(guò)對(duì)常微分方程的近似計(jì)算得到微分方程的近似計(jì)算得到的。本章主要介紹常微分方程數(shù)值解的。本章主要介紹常微分方程數(shù)值解的若干方法。的若干方法。1 1、常微分方程與解、常微分方程與解為為n n階常微分方程階常微分方程。0 )

3、, , ,()(nyyyyxF如果函數(shù)如果函數(shù) 在區(qū)間在區(qū)間a,ba,b內(nèi)內(nèi)n n階可導(dǎo)，稱(chēng)方程階可導(dǎo)，稱(chēng)方程)(xyy )(xyy 滿足方程的函數(shù)滿足方程的函數(shù)稱(chēng)為微分方程的稱(chēng)為微分方程的解解。則則如如為任意常數(shù)）為任意常數(shù)）xy2 CCxy(, 2一般稱(chēng)為方程的一般稱(chēng)為方程的通解通解。為方程的解為方程的解, ,12 xy如果如果則有則有10 )(y為方程滿足定解條件的解。為方程滿足定解條件的解。一一、初值問(wèn)題的數(shù)值解法初值問(wèn)題的數(shù)值解法 102)(yxy CCxy1212 xy方程的通解方程的通解滿足定解條件的解滿足定解條件的解微分關(guān)系（方程）微分關(guān)系（方程）解的圖示解的圖示本教材重點(diǎn)討論

4、定解問(wèn)題本教材重點(diǎn)討論定解問(wèn)題( (初值問(wèn)題）初值問(wèn)題）定解條件（初始條件）定解條件（初始條件） 00yxyyxfy)(),(),(yxf是否能夠找到定解問(wèn)題的解取決于是否能夠找到定解問(wèn)題的解取決于僅有極少數(shù)的方程可以通過(guò)僅有極少數(shù)的方程可以通過(guò)“常數(shù)變易法常數(shù)變易法”、“可分可分離變量法離變量法”等特殊方法求得初等函數(shù)形式的解，絕大等特殊方法求得初等函數(shù)形式的解，絕大部分方程至今無(wú)法理論求解。部分方程至今無(wú)法理論求解。如如xyxyeyxyxyy 212,),sin(sin等等等等2 2、數(shù)值解的思想、數(shù)值解的思想（1 1）將連續(xù)變量）將連續(xù)變量 離散為離散為,bax bxxxxank 10n

5、kxyykk,)(21 （2 2）用代數(shù)的方法求出解函數(shù)）用代數(shù)的方法求出解函數(shù) 在在 點(diǎn)的近似值點(diǎn)的近似值)(xyy kx)(kxy* *ky)(xyy 數(shù)學(xué)界關(guān)注數(shù)學(xué)界關(guān)注工程師關(guān)注工程師關(guān)注如果找不到解函數(shù)如果找不到解函數(shù)數(shù)學(xué)界還關(guān)注：數(shù)學(xué)界還關(guān)注：解的存在性解的存在性解的唯一性解的唯一性解的光滑性解的光滑性解的振動(dòng)性解的振動(dòng)性解的周期性解的周期性解的穩(wěn)定性解的穩(wěn)定性解的混沌性解的混沌性 求函數(shù)求函數(shù) y(x) 在一系列節(jié)點(diǎn)在一系列節(jié)點(diǎn) a = x0 x1 xn= b 處的近似值處的近似值 的方法稱(chēng)為微分方程的數(shù)值解法。的方法稱(chēng)為微分方程的數(shù)值解法。() (1,., )iiyy xin稱(chēng)

6、節(jié)點(diǎn)間距稱(chēng)節(jié)點(diǎn)間距 為步長(zhǎng)，為步長(zhǎng)，通常采用通常采用等距節(jié)點(diǎn)等距節(jié)點(diǎn)，即取，即取 hi = h (常數(shù)常數(shù))。) 1,., 0(1 nixxhiii1,nyy稱(chēng)為微分方程的數(shù)值解。稱(chēng)為微分方程的數(shù)值解。所謂數(shù)值解法：所謂數(shù)值解法：稱(chēng)稱(chēng) 在區(qū)域在區(qū)域D上對(duì)上對(duì) 滿足滿足Lipschitz條件條件是指是指:1212120. .( ,)( ,), , , ( ),( )Ls tf x yf x yL yyxa byyy xy x ( , )f x yy( , ), ( )( )Dx y axb y xyy x記記3 3、相關(guān)定義、相關(guān)定義(2) 一般構(gòu)造方法：一般構(gòu)造方法：4、 迭代格式的構(gòu)造迭代格

7、式的構(gòu)造(1) 構(gòu)造思想：構(gòu)造思想：將連續(xù)的微分方程及初值條件離散為線性方程將連續(xù)的微分方程及初值條件離散為線性方程組加以求解。由于離散化的出發(fā)點(diǎn)不同，產(chǎn)生出各種不同的數(shù)組加以求解。由于離散化的出發(fā)點(diǎn)不同，產(chǎn)生出各種不同的數(shù)值方法?；痉椒ㄓ校河邢薏罘址ǎ〝?shù)值微分）、有限體積法值方法?；痉椒ㄓ校河邢薏罘址ǎ〝?shù)值微分）、有限體積法（數(shù)值積分）、有限元法（函數(shù)插值）等等。（數(shù)值積分）、有限元法（函數(shù)插值）等等。 (3) 如何保證迭代公式的穩(wěn)定性與收斂性如何保證迭代公式的穩(wěn)定性與收斂性?5、微分方程的數(shù)值解法需要解決的主要問(wèn)題、微分方程的數(shù)值解法需要解決的主要問(wèn)題(1) 如何將微分方程離散化，并建

8、立求其如何將微分方程離散化，并建立求其數(shù)值解的迭代公式？數(shù)值解的迭代公式？(2) 如何估計(jì)迭代公式的局部截?cái)嗾`差與整體誤差？如何估計(jì)迭代公式的局部截?cái)嗾`差與整體誤差？二、初值問(wèn)題解的存在唯一性二、初值問(wèn)題解的存在唯一性 考慮一階常微分方程的初值問(wèn)題考慮一階常微分方程的初值問(wèn)題 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy| ),(),(|2121yyLyxfyxf 則上述則上述IVP存在唯一解。存在唯一解。只要只要 在在 上連續(xù)上連續(xù), 且關(guān)于且關(guān)于 y 滿足滿足 Lipschitz 條件，條件，( , )f x y1, a bR即存在與即

9、存在與 無(wú)關(guān)的常數(shù)無(wú)關(guān)的常數(shù) L 使使, x y對(duì)任意定義在對(duì)任意定義在 上的上的 都成立，都成立，, a b 12,yxyx三三、初值問(wèn)題的離散化方法初值問(wèn)題的離散化方法 離散化方法的基本特點(diǎn)是依照某一遞推公式，離散化方法的基本特點(diǎn)是依照某一遞推公式，值值 ,取取 。按節(jié)點(diǎn)從左至右的順序依次求出按節(jié)點(diǎn)從左至右的順序依次求出 的近似的近似( )iy x(1,., )iyin0y 如果計(jì)算如果計(jì)算 ，只用到前一步的值，只用到前一步的值 ，則稱(chēng)這則稱(chēng)這類(lèi)方法為類(lèi)方法為單步方法單步方法。1iyiy如果計(jì)算如果計(jì)算 需用到前需用到前r步的值步的值 , ,則稱(chēng)這類(lèi)方法為則稱(chēng)這類(lèi)方法為r步方法步方法。1i

10、y11,ii ryy iy6.2 Euler6.2 Euler方法方法kp0p1p1npnpkx0 x1x1nxnx),(111212yxfxxyy 第一步：連續(xù)變量離散化第一步：連續(xù)變量離散化,nkxxxxx10第二步：用直線步進(jìn)第二步：用直線步進(jìn)),(000101yxfxxyy ),(),(nnnnnnnnnnyxhfyyyxfxxyy 111EulerEuler格式格式1 1、EulerEuler格式格式 00yxyyxfy)(),(l 18 18世紀(jì)最杰出的數(shù)學(xué)家之一，世紀(jì)最杰出的數(shù)學(xué)家之一，1313歲歲時(shí)入讀巴塞爾大學(xué)，時(shí)入讀巴塞爾大學(xué)，1515歲大學(xué)畢業(yè)，歲大學(xué)畢業(yè)，1616歲獲得

11、碩士學(xué)位。歲獲得碩士學(xué)位。l 17271727年年-1741-1741年（年（2020歲歲-34-34歲）在彼歲）在彼得堡科學(xué)院從事研究工作，在分析學(xué)、得堡科學(xué)院從事研究工作，在分析學(xué)、數(shù)論、力學(xué)方面均有出色成就，并應(yīng)數(shù)論、力學(xué)方面均有出色成就，并應(yīng)俄國(guó)政府要求，解決了不少地圖學(xué)、俄國(guó)政府要求，解決了不少地圖學(xué)、造船業(yè)等實(shí)際問(wèn)題。造船業(yè)等實(shí)際問(wèn)題。l 2424歲晉升物理學(xué)教授。歲晉升物理學(xué)教授。l 17351735年（年（2828歲）右眼失明。歲）右眼失明。l 1741 1741年年 - 1766- 1766（3434歲歲-59-59歲）任德國(guó)科學(xué)院物理數(shù)學(xué)所所歲）任德國(guó)科學(xué)院物理數(shù)學(xué)所所長(zhǎng)，

12、任職長(zhǎng)，任職2525年。在行星運(yùn)動(dòng)、剛體運(yùn)動(dòng)、熱力學(xué)、彈道學(xué)、人年。在行星運(yùn)動(dòng)、剛體運(yùn)動(dòng)、熱力學(xué)、彈道學(xué)、人口學(xué)、微分方程、曲面微分幾何等研究領(lǐng)域均有開(kāi)創(chuàng)性的工作?？趯W(xué)、微分方程、曲面微分幾何等研究領(lǐng)域均有開(kāi)創(chuàng)性的工作。l 17661766年應(yīng)沙皇禮聘重回彼得堡，在年應(yīng)沙皇禮聘重回彼得堡，在17711771年（年（6464歲）左眼失歲）左眼失明。明。l EulerEuler是數(shù)學(xué)史上最多產(chǎn)的數(shù)學(xué)家，平均以每年是數(shù)學(xué)史上最多產(chǎn)的數(shù)學(xué)家，平均以每年800800頁(yè)的速頁(yè)的速度寫(xiě)出創(chuàng)造性論文。他去世后，人們用度寫(xiě)出創(chuàng)造性論文。他去世后，人們用3535年整理出他的研究成年整理出他的研究成果果7474卷。

13、卷。 在假設(shè)在假設(shè) yi = y(xi)，即第，即第 i 步計(jì)算是精確的前提步計(jì)算是精確的前提下，考慮的截?cái)嗾`差下，考慮的截?cái)嗾`差 Ri = y(xi+1) yi+1 稱(chēng)為稱(chēng)為局部截?cái)嗑植拷財(cái)嗾`差誤差 /* local truncation error */。定義定義2.2 若某算法的局部截?cái)嗾`差為若某算法的局部截?cái)嗾`差為O(hp+1)，則稱(chēng)該，則稱(chēng)該 算法有算法有p 階精度。階精度。定義定義2.12、歐拉法的局部截?cái)嗾`差、歐拉法的局部截?cái)嗾`差 歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：11()iiiRy xy23()()2ihyxO hRi 的的主項(xiàng)主項(xiàng)/* leading term */

14、歐拉法具有歐拉法具有 1 階精度階精度。232 ( )( )( )() ( ,)hiiiiiiy xhy xy xO hyhf x y( )iiyy x( )( , ( )iiiy xf x y x2()O h例例1:1: 用歐拉公式求解初值問(wèn)題用歐拉公式求解初值問(wèn)題 2201.201yxyxy （）取步長(zhǎng)取步長(zhǎng) 。 0.1h 解解: : 應(yīng)用應(yīng)用EulerEuler公式于題給初值問(wèn)題的具體形式為：公式于題給初值問(wèn)題的具體形式為： 2120,1,.,1101iiiiyyhx yiy 其中其中 。0.1ixi計(jì)算結(jié)果列于下表：計(jì)算結(jié)果列于下表： iixiy iy xiiy xy 12345678

15、91011120.10.20.30.40.50.60.70.80.91.01.11.21.0000000.9800000.9415840.8883890.8252500.7571470.6883540.6220180.5601130.5036420.4529110.4077830.9900990.9615380.9174310.8630690.8000000.7352940.6711410.6097560.5524860.5000000.4524890.4098360.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620

16、.0076260.0036420.0004220.002053可用來(lái)檢驗(yàn)近似解的準(zhǔn)確程度?？捎脕?lái)檢驗(yàn)近似解的準(zhǔn)確程度。 進(jìn)行計(jì)算，數(shù)值解已達(dá)到了一定的精度。進(jìn)行計(jì)算，數(shù)值解已達(dá)到了一定的精度。這個(gè)初值問(wèn)題的準(zhǔn)確解為這個(gè)初值問(wèn)題的準(zhǔn)確解為 ， 21 1y xx從上表最后一列，我們看到取步長(zhǎng)從上表最后一列，我們看到取步長(zhǎng)0.1h 3、 歐拉公式的改進(jìn)：歐拉公式的改進(jìn)： 隱式歐拉法隱式歐拉法 /* implicit Euler method */向后差商近似導(dǎo)數(shù)向后差商近似導(dǎo)數(shù)hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy 由于未知數(shù)由于未知數(shù) yi+1 同時(shí)出現(xiàn)在

17、等式的兩邊，不能直接同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱(chēng)為得到，故稱(chēng)為隱式隱式 /* implicit */ 歐拉公式，而前者歐拉公式，而前者稱(chēng)為稱(chēng)為顯式顯式 /* explicit */ 歐拉公式。歐拉公式。111,0,1iiiiyhf xiyyn 一般先用顯式計(jì)算一個(gè)初值，再一般先用顯式計(jì)算一個(gè)初值，再迭代迭代求解。求解。隱式隱式歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：11()iiiRy xy23( )()2ihy xO h即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。 梯形公式梯形公式 / /* *trapezoid formula trapezoid formula

18、 * */ / 顯、隱式兩種算法的顯、隱式兩種算法的平均平均) 1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：梯形公式的局部截?cái)嗾`差梯形公式的局部截?cái)嗾`差 ，311iiiRy xyO h即梯形公式即梯形公式具有具有2 階精度階精度，比歐拉方法有了進(jìn)步。，比歐拉方法有了進(jìn)步。但注意到該公式是但注意到該公式是隱式公式隱式公式，計(jì)算時(shí)不得不用到，計(jì)算時(shí)不得不用到迭代法。迭代法。中點(diǎn)歐拉公式中點(diǎn)歐拉公式 /* midpoint formula */中心差商近似導(dǎo)數(shù)中心差商近似導(dǎo)數(shù)hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyx

19、y 1,., 1),(211 niyxfhyyiiii假設(shè)假設(shè) , 則可以導(dǎo)出則可以導(dǎo)出即中點(diǎn)公式具有即中點(diǎn)公式具有 2 階精度。階精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 方方 法法 顯式歐拉顯式歐拉隱式歐拉隱式歐拉梯形公式梯形公式中點(diǎn)公式中點(diǎn)公式簡(jiǎn)單簡(jiǎn)單精度低精度低穩(wěn)定性最好穩(wěn)定性最好精度低精度低, 計(jì)算量大計(jì)算量大精度提高精度提高計(jì)算量大計(jì)算量大精度提高精度提高, 顯式顯式多一個(gè)初值多一個(gè)初值, 可能影響精度可能影響精度 改進(jìn)歐拉法改進(jìn)歐拉法 /* modified Eulers method */Step 1: 先用顯式歐拉公式作預(yù)測(cè)，算出先用顯式歐拉

20、公式作預(yù)測(cè)，算出Step 2: 再將再將 代入代入隱式隱式梯形公式的右邊作梯形公式的右邊作校正校正，得到，得到1 ny1,nnnnyyhf xy111,2nnnnnnyyyyhf xf x11121211()2(,) (0,1,2,.)Euler, (,)(,(,)2,()nnnnnnnnnnnnnnyykkkhhyyf xyf xh yhf xf xynkhfyxh yk上式還常寫(xiě)成該式稱(chēng)為改進(jìn)方法 亦可寫(xiě)成注注: :此法亦稱(chēng)為此法亦稱(chēng)為預(yù)測(cè)預(yù)測(cè)- -校正法校正法 / /* * predictor-corrector method predictor-corrector method * *

21、/ /可以證明該算法可以證明該算法具有具有 2 階精度階精度，同時(shí)可以看到它，同時(shí)可以看到它是個(gè)是個(gè)單步單步遞推格式，比隱式公式的迭代求解過(guò)程遞推格式，比隱式公式的迭代求解過(guò)程簡(jiǎn)單簡(jiǎn)單。后面將看到，它的。后面將看到，它的穩(wěn)定性高于穩(wěn)定性高于顯式歐拉法。顯式歐拉法。改進(jìn)的歐拉法改進(jìn)的歐拉法111(,)(,)2nnnnnnhyyf xyf xy在實(shí)際計(jì)算時(shí)，可將歐拉法與梯形法則相結(jié)合，在實(shí)際計(jì)算時(shí)，可將歐拉法與梯形法則相結(jié)合，計(jì)算公式為計(jì)算公式為 ,.2, 1, 0),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyyknnnnnknnnnn應(yīng)用改進(jìn)歐拉法應(yīng)用改進(jìn)歐拉法, ,

22、如果序列如果序列 收斂收斂, ,)1(1)0(1 nnyy它的極限便滿足方程它的極限便滿足方程111(,)(,)2nnnnnnhyyf xyf xy改進(jìn)歐拉法的截?cái)嗾`差改進(jìn)歐拉法的截?cái)嗾`差)(0)(311hyxynn 因此，改進(jìn)歐拉法公式具有因此，改進(jìn)歐拉法公式具有 2 2 階精度階精度例例2:2: 用改進(jìn)用改進(jìn)Euler公式求解例公式求解例1中的初值問(wèn)題，中的初值問(wèn)題， 取步長(zhǎng)取步長(zhǎng) 。0.1h 解：解：對(duì)此初值問(wèn)題采用改進(jìn)對(duì)此初值問(wèn)題采用改進(jìn)EulerEuler公式，公式， 其具體形式為其具體形式為 21( )211111111( ,)2(,)2()12piiiiiiicppiiiiiii

23、pciiiyyhf x yyhx yyyhf xyyhxyyyy 計(jì)算結(jié)果列于下表：計(jì)算結(jié)果列于下表：例例1:1: 用歐拉公式求解初值問(wèn)題用歐拉公式求解初值問(wèn)題 2201.201yxyxy （）0,1,.,11i 01y iixiy 1piy 1ciy iiy xy改進(jìn)的改進(jìn)的Euler法法 iiy xyEuler法法01234567891011120.00.10.20.30.40.50.60.70.80.91.01.11.21.0000000.9900000.9613660.9172460.8619540.8000340.7355270.6175870.6103990.5532890.500

24、9190.4534790.4108591.0000000.9703890.9243970.8667650.8025170.7360290.6706070.6084430.5507850.4981860.4507350.4082370.9800000.9523330.9100950.8571430.7975510.7350250.6725670.6123550.5557930.5036510.4562230.4134810.0000000.0000990.0001730.0001850.0001150.0000340.0002330.0004460.0006430.0008030.0009190

25、.0009900.0010230.0000000.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.0020532(0)101EulerEuler,12dyxydxyyxyx練習(xí) 設(shè)初值問(wèn)題試分別用法和改進(jìn)法求解 并與精確解進(jìn)行比較。 xEuler法y改進(jìn)的Euler法y精確解01.000000 1.0000001.0000000.11.000000 1.0959091.0954450.21.191818 1.1840971.1832160.31.277438 1.26

26、62011.2649110.41.358213 1.3433601.3416410.51.435133 1.4164021.4142140.61.508966 1.4859561.4832400.71.580338 1.5525141.5491930.81.649783 1.6164751.6124520.91.717779 1.6781661.6733201.01.784770 1.7378671.732051通過(guò)計(jì)算結(jié)果的比較可以看出，改進(jìn)的通過(guò)計(jì)算結(jié)果的比較可以看出，改進(jìn)的Euler方法方法的計(jì)算精度比的計(jì)算精度比Euler方法要高。方法要高。Euler方法、隱式方法、隱式Euler方法

27、、梯形方法與單步法計(jì)算公式的方法、梯形方法與單步法計(jì)算公式的顯式單步法顯式單步法對(duì)應(yīng)關(guān)系對(duì)應(yīng)關(guān)系),(1hyxhyykkkk 隱式單步法隱式單步法),(11hyyxhyykkkkk ),(),(kkkkyxfhyx ),(1kkkkkyxfhyy 顯式顯式 Euler方法方法隱式隱式Euler方法方法),(111 kkkkkyxfhyy),(),(2111 kkkkkkkyxfyxfhyy梯形方法梯形方法(隱式隱式),(),(111 kkkkkyxfhyyx),(21),(21),(111 kkkkkkkyxfyxfhyyx單步法小結(jié)單步法小結(jié) 收斂性收斂性 /* Convergency */

28、定義定義 若某算法對(duì)于任意固定的若某算法對(duì)于任意固定的 x = xi = x0 + i h，當(dāng)，當(dāng) h0 ( 同時(shí)同時(shí) i ) 時(shí)有時(shí)有 yi y( xi )，則稱(chēng)該算法是，則稱(chēng)該算法是收斂的。收斂的。 收斂性定理收斂性定理若某單步法滿足以上條件，則該方法是收斂的若某單步法滿足以上條件，則該方法是收斂的則該單步法的整體截?cái)嗾`差為：則該單步法的整體截?cái)嗾`差為：若單步法若單步法 具有具有 p 階精度，階精度，且增量函數(shù)且增量函數(shù) 關(guān)于關(guān)于 y 滿足：滿足：1(, )nnnnyyhxy h ( , , )x y h ()()pnny xyO h( , , )( , , )x y hx y hL yy

29、 Lipschitz 條件：條件： 初值初值 y0 是準(zhǔn)確的是準(zhǔn)確的 穩(wěn)定性穩(wěn)定性 /* Stability */例：例：考察初值問(wèn)題考察初值問(wèn)題 在區(qū)間在區(qū)間0, 0.5上的解上的解.分別用歐拉顯、隱式格式和改進(jìn)的歐拉格分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解。式計(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.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 1

30、0 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 ?! An Engineer complains: Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!定義定

31、義若某算法在計(jì)算過(guò)程中任一步產(chǎn)生的誤差在若某算法在計(jì)算過(guò)程中任一步產(chǎn)生的誤差在以后的計(jì)算中都以后的計(jì)算中都逐步衰減逐步衰減，則稱(chēng)該算法是，則稱(chēng)該算法是絕對(duì)穩(wěn)定絕對(duì)穩(wěn)定的的 /*absolutely stable */。一般分析時(shí)為簡(jiǎn)單起見(jiàn)，只考慮一般分析時(shí)為簡(jiǎn)單起見(jiàn)，只考慮試驗(yàn)方程試驗(yàn)方程 /* test equation */ 0yy 常數(shù)，可以常數(shù)，可以是復(fù)數(shù)是復(fù)數(shù)我們稱(chēng)我們稱(chēng)算法算法A 比算法比算法B 穩(wěn)定穩(wěn)定，就是指，就是指 A 的絕對(duì)穩(wěn)定的絕對(duì)穩(wěn)定區(qū)域比區(qū)域比 B 的的大大。當(dāng)步長(zhǎng)取為當(dāng)步長(zhǎng)取為 h 時(shí)，將某算法應(yīng)用于上式，并假設(shè)只在時(shí)，將某算法應(yīng)用于上式，并假設(shè)只在初值產(chǎn)生誤差初值

32、產(chǎn)生誤差 ，則若此誤差以后逐步衰減，則若此誤差以后逐步衰減,000yy 就稱(chēng)該算法相對(duì)于就稱(chēng)該算法相對(duì)于 絕對(duì)穩(wěn)定，絕對(duì)穩(wěn)定， 的全體構(gòu)成的全體構(gòu)成h h h絕對(duì)穩(wěn)定區(qū)域絕對(duì)穩(wěn)定區(qū)域。yy 例：例：考察顯式歐拉法考察顯式歐拉法110(1)iiiiyyhyhy000yy 011)1 (yhyii 01111)1 ( iiiihyy0-1-2ReImg由此可見(jiàn)，要保證初始誤差由此可見(jiàn)，要保證初始誤差 0 以后逐步衰減，以后逐步衰減，必須滿足：必須滿足：hh 1|1| h例：例：考察隱式歐拉法考察隱式歐拉法11 iiiyhyy iiyhy 11101111 iih可見(jiàn)絕對(duì)穩(wěn)定區(qū)域?yàn)椋嚎梢?jiàn)絕對(duì)穩(wěn)定區(qū)域

33、為：1|1 | h210ReImg注：注：一般來(lái)說(shuō)，隱式歐拉法的絕對(duì)穩(wěn)定性比同階一般來(lái)說(shuō)，隱式歐拉法的絕對(duì)穩(wěn)定性比同階的顯式法的好。的顯式法的好。定理定理 某單步法用于解模型方程某單步法用于解模型方程y =y，若得到的解，若得到的解yn+1=E(h)yn，滿足，滿足|E(h)|1，則稱(chēng)該單步法是，則稱(chēng)該單步法是絕對(duì)絕對(duì)穩(wěn)定穩(wěn)定的的. 在平面上在平面上, 使使|E(h)|1的變量圍成的區(qū)域，的變量圍成的區(qū)域，稱(chēng)為稱(chēng)為絕對(duì)穩(wěn)定區(qū)域絕對(duì)穩(wěn)定區(qū)域，它與實(shí)軸的交稱(chēng)為，它與實(shí)軸的交稱(chēng)為絕對(duì)穩(wěn)定區(qū)絕對(duì)穩(wěn)定區(qū)間間.穩(wěn)定性定理穩(wěn)定性定理6.3 龍格庫(kù)塔方法 對(duì)許多實(shí)際問(wèn)題來(lái)說(shuō)，歐拉公式與改進(jìn)歐拉對(duì)許多實(shí)際問(wèn)題來(lái)

34、說(shuō)，歐拉公式與改進(jìn)歐拉公式精度還不能滿足要求，為此從另一個(gè)角度來(lái)分公式精度還不能滿足要求，為此從另一個(gè)角度來(lái)分析這兩個(gè)公式的特點(diǎn)，從而探索一條構(gòu)造高精度方析這兩個(gè)公式的特點(diǎn)，從而探索一條構(gòu)造高精度方法的途徑法的途徑. 受改進(jìn)的Euler方法啟發(fā)，更一般算式可設(shè)為1121211()2(,)(,)(0,1,2,.)nnnnnnyykkkhf xykhf xh ykn123111()-()()nnnnnTy xyO hyy x適當(dāng)選擇參數(shù) ， ， ，使局部截?cái)嗾`差,這里仍假定。改進(jìn)歐拉法改進(jìn)歐拉法11 122121(,)(0,1,2,.)(,)nnnnnnyykkkhf xynkhf xh yk11

35、22232()()()(,)(,)(,)()nnnxnnnnynnyy xhy xfxyf xyfxyhO h)(),(),(),()()(),(),(),(:Taylor323122hOyxfyxfyxfhxyhhOyxfhkyxfhyxhfknnynnnnxnnnynnxnn展開(kāi)式由二元函數(shù)11 122121 (,) (,)(0,1,2,.)nnnnnnnyykkkhf xyhykhf xh ykn由于四個(gè)參數(shù)，三個(gè)方程，因此有一個(gè)自由參數(shù)，即解答不唯一。1122232()()()(,)(,)(,)()nnnxnnnnynnyy xhy xfxyf xyfxyhO hTaylor:與展式相

36、比較得1212122122312()()()()()hnnnny xy xhy xy xO h 這是改進(jìn)的Euler方法。1211(1),1,22取可得此時(shí)算式為1121211()2(,) (,)nnnnnnyykkkhf xykhf xh yk-R K這是二階方法.121(2)0,1,2取可得此時(shí)算式為12121(,)11(,) 22nnnnnnyykkhf xykhf xh yk R-K這也是二階方法。12132(3),443取可得又有算式11212113)4(,)22(,)33nnnnnnyykkkhf xykhf xh yk（三階龍格-庫(kù)塔方法三階龍格-庫(kù)塔方法是用三個(gè)值 k1, k2, k3 的線性組合112312123221331332112(,)(,)(,)nnnnnnnnyykkkkhf xykhf xhcccababykkhf xh ykb k 要使三階龍格-庫(kù)塔方法具有三階精度，必須使其局部截?cái)嗾`差為 O(h4)將 k1, k2, k3 代入 yn+1 的