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

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

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

3、程與解為為n n階常微分方程階常微分方程。0 ), , ,()(nyyyyxF如果函數(shù)如果函數(shù) 在區(qū)間在區(qū)間a,ba,b內(nèi)內(nèi)n n階可導(dǎo)，稱方程階可導(dǎo)，稱方程)(xyy )(xyy 滿足方程的函數(shù)滿足方程的函數(shù)稱為微分方程的稱為微分方程的解解。則則如如為任意常數(shù)）為任意常數(shù)）xy2 CCxy(, 2一般稱為方程的一般稱為方程的通解通解。為方程的解。為方程的解。12 xy如果如果則有則有10 )(y為方程滿足定解條件的解。為方程滿足定解條件的解。一一、初值問題的數(shù)值解法初值問題的數(shù)值解法微分方程數(shù)值解法 102)(yxy CCxy1212 xy方程的通解方程的通解滿足定解條件的解滿足定解條件的解

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

5、、數(shù)值解的思想（1 1）將連續(xù)變量）將連續(xù)變量 離散為離散為,bax bxxxxank 10nkxyykk,)(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ù)求函數(shù) y(x) 在一系列節(jié)點(diǎn)在一系列節(jié)點(diǎn) a = x0 x1 xn= b 處的近似值

6、處的近似值 的方法稱為微分方程的數(shù)值解法。的方法稱為微分方程的數(shù)值解法。() (1,., )iiyy xin稱節(jié)點(diǎn)間距稱節(jié)點(diǎn)間距 為步長，為步長，通常采用通常采用等距節(jié)點(diǎn)等距節(jié)點(diǎn)，即取，即取 hi = h (常數(shù)常數(shù))。) 1,., 0(1 nixxhiii1,nyy稱為微分方程的數(shù)值解。稱為微分方程的數(shù)值解。所謂數(shù)值解法：所謂數(shù)值解法：微分方程數(shù)值解法稱稱 在區(qū)域在區(qū)域D上對(duì)上對(duì) 滿足滿足Lipschitz條件條件是指是指:1212120. .( ,)( ,), , , ( ),( )Ls tf x yf x yL yyxa byyy xy x ( , )f x yy( , ), ( )(

7、)Dx y axb y xyy x記記3 3、相關(guān)定義、相關(guān)定義微分方程數(shù)值解法(2) 一般構(gòu)造方法：一般構(gòu)造方法：4、 迭代格式的構(gòu)造迭代格式的構(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ù)插值）等等。 微分方程數(shù)值解法(3) 如何保證迭代公式的穩(wěn)

8、定性與收斂性如何保證迭代公式的穩(wěn)定性與收斂性?5、微分方程的數(shù)值解法需要解決的主要問題、微分方程的數(shù)值解法需要解決的主要問題(1) 如何將微分方程離散化，并建立求其如何將微分方程離散化，并建立求其數(shù)值解的迭代公式？數(shù)值解的迭代公式？(2) 如何估計(jì)迭代公式的局部截?cái)嗾`差與整體誤差？如何估計(jì)迭代公式的局部截?cái)嗾`差與整體誤差？微分方程數(shù)值解法二、初值問題解的存在唯一性二、初值問題解的存在唯一性 考慮一階常微分方程的初值問題考慮一階常微分方程的初值問題 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy| ),(),(|2121yyLyxfyxf 則上

9、述則上述IVP存在唯一解。存在唯一解。只要只要 在在 上連續(xù)上連續(xù), 且關(guān)于且關(guān)于 y 滿足滿足 Lipschitz 條件，條件，( , )f x y1, a bR即存在與即存在與 無關(guān)的常數(shù)無關(guān)的常數(shù) L 使使, x y對(duì)任意定義在對(duì)任意定義在 上的上的 都成立，都成立，, a b 12,yxyx微分方程數(shù)值解法三、初值問題的離散化方法三、初值問題的離散化方法 離散化方法的基本特點(diǎn)是依照某一遞推公式，離散化方法的基本特點(diǎn)是依照某一遞推公式，值值 ,取取 。按節(jié)點(diǎn)從左至右的順序依次求出按節(jié)點(diǎn)從左至右的順序依次求出 的近似的近似( )iy x(1,., )iyin0y 如果計(jì)算如果計(jì)算 ，只用到

10、前一步的值，只用到前一步的值 ，則稱這則稱這類方法為類方法為單步方法單步方法。1iyiy如果計(jì)算如果計(jì)算 需用到前需用到前r步的值步的值 , ,則稱這類方法為則稱這類方法為r步方法步方法。1iy11,ii ryy iy微分方程數(shù)值解法6.2 Euler6.2 Euler方法方法kp0p1p1npnpkx0 x1x1nxnx),(111212yxfxxyy 第一步：連續(xù)變量離散化第一步：連續(xù)變量離散化,nkxxxxx10第二步：用直線步進(jìn)第二步：用直線步進(jìn)),(000101yxfxxyy ),(),(nnnnnnnnnnyxhfyyyxfxxyy 111EulerEuler格式格式1 1、Eul

11、erEuler格式格式 00yxyyxfy)(),(微分方程數(shù)值解法l 18 18世紀(jì)最杰出的數(shù)學(xué)家之一，世紀(jì)最杰出的數(shù)學(xué)家之一，1313歲歲時(shí)入讀巴塞爾大學(xué)，時(shí)入讀巴塞爾大學(xué)，1515歲大學(xué)畢業(yè)，歲大學(xué)畢業(yè)，1616歲獲得碩士學(xué)位。歲獲得碩士學(xué)位。l 17271727年年-1741-1741年（年（2020歲歲-34-34歲）在彼歲）在彼得堡科學(xué)院從事研究工作，在分析學(xué)、得堡科學(xué)院從事研究工作，在分析學(xué)、數(shù)論、力學(xué)方面均有出色成就，并應(yīng)數(shù)論、力學(xué)方面均有出色成就，并應(yīng)俄國政府要求，解決了不少地圖學(xué)、俄國政府要求，解決了不少地圖學(xué)、造船業(yè)等實(shí)際問題。造船業(yè)等實(shí)際問題。l 2424歲晉升物理學(xué)教

12、授。歲晉升物理學(xué)教授。l 17351735年（年（2828歲）右眼失明。歲）右眼失明。微分方程數(shù)值解法l 1741 1741年年 - 1766- 1766（3434歲歲-59-59歲）任德國科學(xué)院物理數(shù)學(xué)所所歲）任德國科學(xué)院物理數(shù)學(xué)所所長，任職長，任職2525年。在行星運(yùn)動(dòng)、剛體運(yùn)動(dòng)、熱力學(xué)、彈道學(xué)、人年。在行星運(yùn)動(dòng)、剛體運(yùn)動(dòng)、熱力學(xué)、彈道學(xué)、人口學(xué)、微分方程、曲面微分幾何等研究領(lǐng)域均有開創(chuàng)性的工作?？趯W(xué)、微分方程、曲面微分幾何等研究領(lǐng)域均有開創(chuàng)性的工作。l 17661766年應(yīng)沙皇禮聘重回彼得堡，在年應(yīng)沙皇禮聘重回彼得堡，在17711771年（年（6464歲）左眼失歲）左眼失明。明。l Eu

13、lerEuler是數(shù)學(xué)史上最多產(chǎn)的數(shù)學(xué)家，平均以每年是數(shù)學(xué)史上最多產(chǎn)的數(shù)學(xué)家，平均以每年800800頁的速頁的速度寫出創(chuàng)造性論文。他去世后，人們用度寫出創(chuàng)造性論文。他去世后，人們用3535年整理出他的研究成年整理出他的研究成果果7474卷。卷。 微分方程數(shù)值解法在假設(shè)在假設(shè) yi = y(xi)，即第，即第 i 步計(jì)算是精確的前提步計(jì)算是精確的前提下，考慮的截?cái)嗾`差下，考慮的截?cái)嗾`差 Ri = y(xi+1) yi+1 稱為稱為局部截?cái)嗑植拷財(cái)嗾`差誤差 /* local truncation error */。定義定義2.2 若某算法的局部截?cái)嗾`差為若某算法的局部截?cái)嗾`差為O(hp+1)，則稱

14、該，則稱該 算法有算法有p 階精度。階精度。定義定義2.12、歐拉法的局部截?cái)嗾`差、歐拉法的局部截?cái)嗾`差微分方程數(shù)值解法 歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：11()iiiRy xy23()()2ihyxO hRi 的的主項(xiàng)主項(xiàng)/* leading term */歐拉法具有歐拉法具有 1 階精階精度度。232 ( )( )( )() ( ,)hiiiiiiy xhy xy xO hyhf x y( )iiyy x( )( , ( )iiiy xf x y x2()O h微分方程數(shù)值解法例例1:1: 用歐拉公式求解初值問題用歐拉公式求解初值問題 2201.201yxyxy （）取步長取步

15、長 。 0.1h 解解: : 應(yīng)用應(yīng)用EulerEuler公式于題給初值問題的具體形式為：公式于題給初值問題的具體形式為： 2120,1,.,1101iiiiyyhx yiy 其中其中 。0.1ixi計(jì)算結(jié)果列于下表：計(jì)算結(jié)果列于下表： 微分方程數(shù)值解法iixiy iy xiiy xy 1234567891011120.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

16、.9615380.9174310.8630690.8000000.7352940.6711410.6097560.5524860.5000000.4524890.4098360.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.002053微分方程數(shù)值解法可用來檢驗(yàn)近似解的準(zhǔn)確程度?？捎脕頇z驗(yàn)近似解的準(zhǔn)確程度。 進(jìn)行計(jì)算，數(shù)值解已達(dá)到了一定的精度。進(jìn)行計(jì)算，數(shù)值解已達(dá)到了一定的精度。這個(gè)初值問題的準(zhǔn)確解為這個(gè)初值問題的準(zhǔn)確解為 ， 21 1y xx從上表最后一列，我們看

17、到取步長從上表最后一列，我們看到取步長0.1h 微分方程數(shù)值解法3、 歐拉公式的改進(jìn)：歐拉公式的改進(jìn)： 隱式歐拉法隱式歐拉法 /* implicit Euler method */向后差商近似導(dǎo)數(shù)向后差商近似導(dǎo)數(shù)hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy 微分方程數(shù)值解法由于未知數(shù)由于未知數(shù) yi+1 同時(shí)出現(xiàn)在等式的兩邊，不能直接同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱為得到，故稱為隱式隱式 /* implicit */ 歐拉公式，而前者歐拉公式，而前者稱為稱為顯式顯式 /* explicit */ 歐拉公式。歐拉公式。111,0,1iiiiyhfxiy

18、yn微分方程數(shù)值解法一般先用顯式計(jì)算一個(gè)初值，再一般先用顯式計(jì)算一個(gè)初值，再迭代迭代求解。求解。隱式隱式歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：11()iiiRy xy23( )()2ihy xO h即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。微分方程數(shù)值解法 梯形公式梯形公式 / /* *trapezoid formula trapezoid formula * */ / 顯、隱式兩種算法的顯、隱式兩種算法的平均平均) 1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：梯形公式的局部截?cái)嗾`差梯形公式的局部截?cái)嗾`差 ，311iiiRy xyO h即

19、梯形公式即梯形公式具有具有2 階精度階精度，比歐拉方法有了進(jìn)步。，比歐拉方法有了進(jìn)步。但注意到該公式是但注意到該公式是隱式公式隱式公式，計(jì)算時(shí)不得不用到，計(jì)算時(shí)不得不用到迭代法，其迭代收斂性與歐拉公式相似。迭代法，其迭代收斂性與歐拉公式相似。微分方程數(shù)值解法中點(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 階精度。階精度。)(),(

20、11iiiixyyxyy )()(311hOyxyRiii 微分方程數(shù)值解法方方 法法 顯式歐拉顯式歐拉隱式歐拉隱式歐拉梯形公式梯形公式中點(diǎn)公式中點(diǎn)公式簡單簡單精度低精度低穩(wěn)定性最好穩(wěn)定性最好精度低精度低, 計(jì)算量大計(jì)算量大精度提高精度提高計(jì)算量大計(jì)算量大精度提高精度提高, 顯式顯式多一個(gè)初值多一個(gè)初值, 可能影響精度可能影響精度微分方程數(shù)值解法 改進(jìn)歐拉法改進(jìn)歐拉法 /* modified Eulers method */Step 1: 先用顯式歐拉公式作預(yù)測，算出先用顯式歐拉公式作預(yù)測，算出Step 2: 再將再將 代入代入隱式隱式梯形公式的右邊作梯形公式的右邊作校正校正，得到，得到1 n

21、y1,nnnnyyhf xy111,2nnnnnnyyyyhf xf x微分方程數(shù)值解法注注: :此法亦稱為此法亦稱為預(yù)測預(yù)測- -校正法校正法 / /* * predictor-corrector method predictor-corrector method * */ /可以證明該算法可以證明該算法具有具有 2 階精度階精度，同時(shí)可以看到它，同時(shí)可以看到它是個(gè)是個(gè)單步單步遞推格式，比隱式公式的迭代求解過程遞推格式，比隱式公式的迭代求解過程簡單簡單。后面將看到，它的。后面將看到，它的穩(wěn)定性高于穩(wěn)定性高于顯式歐拉法。顯式歐拉法。改進(jìn)的歐拉法改進(jìn)的歐拉法111(,)(,)2nnnnnnhyy

22、f xyf xy微分方程數(shù)值解法在實(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)歐拉法, ,如果序列如果序列 收斂收斂, ,)1(1)0(1 nnyy它的極限便滿足方程它的極限便滿足方程111(,)(,)2nnnnnnhyyf xyf xy微分方程數(shù)值解法改進(jìn)歐拉法的截?cái)嗾`差改進(jìn)歐拉法的截?cái)嗾`差)(0)(311hyxynn 因此，改進(jìn)歐拉法公式具有因此，改進(jìn)歐拉法公式具有 2 2 階精度階精度微分方

23、程數(shù)值解法例例2:2: 用改進(jìn)用改進(jìn)Euler公式求解例公式求解例1中的初值問題，中的初值問題， 取步長取步長 。0.1h 解：解：對(duì)此初值問題采用改進(jìn)對(duì)此初值問題采用改進(jìn)EulerEuler公式，公式， 其具體形式為其具體形式為 21( )211111111( ,)2(,)2()12piiiiiiicppiiiiiiipciiiyyhf x yyhx yyyhf xyyhxyyyy 計(jì)算結(jié)果列于下表：計(jì)算結(jié)果列于下表：例例1:1: 用歐拉公式求解初值問題用歐拉公式求解初值問題 2201.201yxyxy （）0,1,.,11i 01y 微分方程數(shù)值解法iixiy 1piy 1ciyiiy x

24、y改進(jìn)的改進(jìn)的Euler法法iiy xyEuler法法01234567891 01 11 20.00.10.20.30.40.50.60.70.80.91.01.11.21.0000000.9900000.9613660.9172460.8619540.8000340.7355270.6175870.6103990.5532890.5009190.4534790.4108591.0000000.9703890.9243970.8667650.8025170.7360290.6706070.6084430.5507850.4981860.4507350.4082370.9800000.95233

25、30.9100950.8571430.7975510.7350250.6725670.6123550.5557930.5036510.4562230.4134810.0000000.0000990.0001730.0001850.0001150.0000340.0002330.0004460.0006430.0008030.0009190.0009900.0010230.0000000.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.002053微分方程數(shù)值解法通過計(jì)

26、算結(jié)果的比較可以看出，改進(jìn)的通過計(jì)算結(jié)果的比較可以看出，改進(jìn)的Euler方法方法的計(jì)算精度比的計(jì)算精度比Euler方法要高。方法要高。微分方程數(shù)值解法歐拉法誤差概述)()(2)O()(y2Euler式對(duì)于)O(h)(2Euler對(duì)于顯式則精度相對(duì)的越高。越大，誤差階，一個(gè)算法，局部截?cái)嘁话?31321321hOxyhThxhTxyhTnnnnnn ：對(duì)于梯形公式：法隱：法說來微分方程數(shù)值解法6.3 龍格庫塔方法 對(duì)許多實(shí)際問題來說，歐拉公式與改進(jìn)歐拉對(duì)許多實(shí)際問題來說，歐拉公式與改進(jìn)歐拉公式精度還不能滿足要求，為此從另一個(gè)角度來分公式精度還不能滿足要求，為此從另一個(gè)角度來分析這兩個(gè)公式的特點(diǎn)，

27、從而探索一條構(gòu)造高精度方析這兩個(gè)公式的特點(diǎn)，從而探索一條構(gòu)造高精度方法的途徑法的途徑. 微分方程數(shù)值解法受改進(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 yk微分方程數(shù)值解法 這是改進(jìn)的Euler方法。1211(1),1,22取可得此時(shí)算式為1121211()2(,) (,)

28、nnnnnnyykkkhf xykhf xh yk微分方程數(shù)值解法 R-K這也是二階方法。12132(3),443取可得又有算式11212113)4(,)22(,)33nnnnnnyykkkhf xykhf xh yk（微分方程數(shù)值解法三階龍格-庫塔方法三階龍格-庫塔方法是用三個(gè)值 k1, k2, k3 的線性組合112312123221331332112(,)(,)(,)nnnnnnnnyykkkkhf xykhf xhcccababykkhf xh ykb k 要使三階龍格-庫塔方法具有三階精度，必須使其局部截?cái)嗾`差為 O(h4)將 k1, k2, k3 代入 yn+1 的表達(dá)式中，在 (xn, yn) 處用二元泰勒公式展開，與 y(xn+1) 在 xn 處的泰勒展開式比較微分方程數(shù)值解法 類似二階龍格-庫塔方法的推導(dǎo)過程，8 個(gè)待定系數(shù) c1, c2, c3, a2, a3, b