1、2021/8/141 第6章 常微分方程數值解法2021/8/142 微分方程數值解一般可分為：微分方程數值解一般可分為：常微分方程數值解常微分方程數值解和和偏微分偏微分方程數值解方程數值解。自然界與工程技術中的許多現象，其數學表達式。自然界與工程技術中的許多現象，其數學表達式可歸結為常微分方程（組）的定解問題。一些偏微分方程問題可歸結為常微分方程（組）的定解問題。一些偏微分方程問題也可以轉化為常微分方程問題來（近似）求解。也可以轉化為常微分方程問題來（近似）求解。Newton最早采最早采用數學方法研究二體問題，其中需要求解的運動方程就是常微用數學方法研究二體問題，其中需要求解的運動方程就是常

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

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

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

5、462 2、數值解的思想、數值解的思想（1 1）將連續變量）將連續變量 離散為離散為,bax bxxxxank 10nkxyykk,)(21 （2 2）用代數的方法求出解函數）用代數的方法求出解函數 在在 點的近似值點的近似值)(xyy kx)(kxy* *ky)(xyy 數學界關注數學界關注工程師關注工程師關注如果找不到解函數如果找不到解函數數學界還關注：數學界還關注：解的存在性解的存在性解的唯一性解的唯一性解的光滑性解的光滑性解的振動性解的振動性解的周期性解的周期性解的穩定性解的穩定性解的混沌性解的混沌性2021/8/147 求函數求函數 y(x) 在一系列節點在一系列節點 a = x0

6、x1 xn= b 處的近似值處的近似值 的方法稱為微分方程的數值解法。的方法稱為微分方程的數值解法。() (1,., )iiyy xin稱節點間距稱節點間距 為步長，為步長，通常采用通常采用等距節點等距節點，即取，即取 hi = h (常數常數)。) 1,., 0(1 nixxhiii1,nyy稱為微分方程的數值解。稱為微分方程的數值解。所謂數值解法：所謂數值解法：2021/8/148稱稱 在區域在區域D上對上對 滿足滿足Lipschitz條件條件是指是指:1212120. .( ,)( ,), , , ( ),( )Ls tf x yf x yL yyxa byyy xy x ( , )f

7、x yy( , ), ( )( )Dx y axb y xyy x記記3 3、相關定義、相關定義2021/8/149(2) 一般構造方法：一般構造方法：4、 迭代格式的構造迭代格式的構造(1) 構造思想：構造思想：將連續的微分方程及初值條件離散為線性方程將連續的微分方程及初值條件離散為線性方程組加以求解。由于離散化的出發點不同，產生出各種不同的數組加以求解。由于離散化的出發點不同，產生出各種不同的數值方法。基本方法有：有限差分法（數值微分）、有限體積法值方法。基本方法有：有限差分法（數值微分）、有限體積法（數值積分）、有限元法（函數插值）等等。（數值積分）、有限元法（函數插值）等等。 2021

8、/8/1410(3) 如何保證迭代公式的穩定性與收斂性如何保證迭代公式的穩定性與收斂性?5、微分方程的數值解法需要解決的主要問題、微分方程的數值解法需要解決的主要問題(1) 如何將微分方程離散化，并建立求其如何將微分方程離散化，并建立求其數值解的迭代公式？數值解的迭代公式？(2) 如何估計迭代公式的局部截斷誤差與整體誤差？如何估計迭代公式的局部截斷誤差與整體誤差？2021/8/1411二、初值問題解的存在唯一性二、初值問題解的存在唯一性 考慮一階常微分方程的初值問題考慮一階常微分方程的初值問題 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy|

9、 ),(),(|2121yyLyxfyxf 則上述則上述IVP存在唯一解。存在唯一解。只要只要 在在 上連續上連續, 且關于且關于 y 滿足滿足 Lipschitz 條件，條件，( , )f x y1, a bR即存在與即存在與 無關的常數無關的常數 L 使使, x y對任意定義在對任意定義在 上的上的 都成立，都成立，, a b 12,yxyx2021/8/1412三、初值問題的離散化方法三、初值問題的離散化方法 離散化方法的基本特點是依照某一遞推公式，離散化方法的基本特點是依照某一遞推公式，值值 ,取取 。按節點從左至右的順序依次求出按節點從左至右的順序依次求出 的近似的近似( )iy x

10、(1,., )iyin0y 如果計算如果計算 ，只用到前一步的值，只用到前一步的值 ，則稱這則稱這類方法為類方法為單步方法單步方法。1iyiy如果計算如果計算 需用到前需用到前r步的值步的值 , ,則稱這類方法為則稱這類方法為r步方法步方法。1iy11,ii ryy iy2021/8/14136.2 Euler6.2 Euler方法方法kp0p1p1npnpkx0 x1x1nxnx),(111212yxfxxyy 第一步：連續變量離散化第一步：連續變量離散化,nkxxxxx10第二步：用直線步進第二步：用直線步進),(000101yxfxxyy ),(),(nnnnnnnnnnyxhfyyyx

11、fxxyy 111EulerEuler格式格式1 1、EulerEuler格式格式 00yxyyxfy)(),(2021/8/1414l 18 18世紀最杰出的數學家之一，世紀最杰出的數學家之一，1313歲歲時入讀巴塞爾大學，時入讀巴塞爾大學，1515歲大學畢業，歲大學畢業，1616歲獲得碩士學位。歲獲得碩士學位。l 17271727年年-1741-1741年（年（2020歲歲-34-34歲）在彼歲）在彼得堡科學院從事研究工作，在分析學、得堡科學院從事研究工作，在分析學、數論、力學方面均有出色成就，并應數論、力學方面均有出色成就，并應俄國政府要求，解決了不少地圖學、俄國政府要求，解決了不少地圖

12、學、造船業等實際問題。造船業等實際問題。l 2424歲晉升物理學教授。歲晉升物理學教授。l 17351735年（年（2828歲）右眼失明。歲）右眼失明。2021/8/1415l 1741 1741年年 - 1766- 1766（3434歲歲-59-59歲）任德國科學院物理數學所所歲）任德國科學院物理數學所所長，任職長，任職2525年。在行星運動、剛體運動、熱力學、彈道學、人年。在行星運動、剛體運動、熱力學、彈道學、人口學、微分方程、曲面微分幾何等研究領域均有開創性的工作。口學、微分方程、曲面微分幾何等研究領域均有開創性的工作。l 17661766年應沙皇禮聘重回彼得堡，在年應沙皇禮聘重回彼得堡

13、，在17711771年（年（6464歲）左眼失歲）左眼失明。明。l EulerEuler是數學史上最多產的數學家，平均以每年是數學史上最多產的數學家，平均以每年800800頁的速頁的速度寫出創造性論文。他去世后，人們用度寫出創造性論文。他去世后，人們用3535年整理出他的研究成年整理出他的研究成果果7474卷。卷。 2021/8/1416在假設在假設 yi = y(xi)，即第，即第 i 步計算是精確的前提步計算是精確的前提下，考慮的截斷誤差下，考慮的截斷誤差 Ri = y(xi+1) yi+1 稱為稱為局部截斷局部截斷誤差誤差 /* local truncation error */。定義定

14、義2.2 若某算法的局部截斷誤差為若某算法的局部截斷誤差為O(hp+1)，則稱該，則稱該 算法有算法有p 階精度。階精度。定義定義2.12、歐拉法的局部截斷誤差、歐拉法的局部截斷誤差2021/8/1417 歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：11()iiiRy xy23()()2ihyxO hRi 的的主項主項/* leading term */歐拉法具有歐拉法具有 1 階精階精度度。232 ( )( )( )() ( ,)hiiiiiiy xhy xy xO hyhf x y( )iiyy x( )( , ( )iiiy xf x y x2()O h2021/8/1418例例1:1:

15、 用歐拉公式求解初值問題用歐拉公式求解初值問題 2201.201yxyxy （）取步長取步長 。 0.1h 解解: : 應用應用EulerEuler公式于題給初值問題的具體形式為：公式于題給初值問題的具體形式為： 2120,1,.,1101iiiiyyhx yiy 其中其中 。0.1ixi計算結果列于下表：計算結果列于下表： 2021/8/1419iixiy iy xiiy xy 1234567891011120.10.20.30.40.50.60.70.80.91.01.11.21.0000000.9800000.9415840.8883890.8252500.7571470.6883540

16、.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.0076260.0036420.0004220.0020532021/8/1420可用來檢驗近似解的準確程度。可用來檢驗近似解的準確程度。 進行計算，數值解已達到了一定的精度。進行計算，數值解已達到

17、了一定的精度。這個初值問題的準確解為這個初值問題的準確解為 ， 21 1y xx從上表最后一列，我們看到取步長從上表最后一列，我們看到取步長0.1h 2021/8/14213、 歐拉公式的改進：歐拉公式的改進： 隱式歐拉法隱式歐拉法 /* implicit Euler method */向后差商近似導數向后差商近似導數hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy 2021/8/1422由于未知數由于未知數 yi+1 同時出現在等式的兩邊，不能直接同時出現在等式的兩邊，不能直接得到，故稱為得到，故稱為隱式隱式 /* implicit */ 歐拉公式，而前者歐拉

18、公式，而前者稱為稱為顯式顯式 /* explicit */ 歐拉公式。歐拉公式。111,0,1iiiiyhfxiyyn2021/8/1423一般先用顯式計算一個初值，再一般先用顯式計算一個初值，再迭代迭代求解。求解。隱式隱式歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：11()iiiRy xy23( )()2ihy xO h即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。2021/8/1424 梯形公式梯形公式 / /* *trapezoid formula trapezoid formula * */ / 顯、隱式兩種算法的顯、隱式兩種算法的平均平均) 1,., 0(),(),(21

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

20、),(211 niyxfhyyiiii假設假設 , 則可以導出則可以導出即中點公式具有即中點公式具有 2 階精度。階精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 2021/8/1426方方 法法 顯式歐拉顯式歐拉隱式歐拉隱式歐拉梯形公式梯形公式中點公式中點公式簡單簡單精度低精度低穩定性最好穩定性最好精度低精度低, 計算量大計算量大精度提高精度提高計算量大計算量大精度提高精度提高, 顯式顯式多一個初值多一個初值, 可能影響精度可能影響精度2021/8/1427 改進歐拉法改進歐拉法 /* modified Eulers method */Step 1: 先用顯式歐拉

21、公式作預測，算出先用顯式歐拉公式作預測，算出Step 2: 再將再將 代入代入隱式隱式梯形公式的右邊作梯形公式的右邊作校正校正，得到，得到1 ny1,nnnnyyhf xy111,2nnnnnnyyyyhf xf x2021/8uler,1()2(,)(,)(0,1,2,.(,)(,(,).)2nnnnnnnnnnnnnnyykkkhf xykhf xh yhyyf xyf xh yhf xykn上式還常寫成該式稱為改進方法 亦可寫成2021/8/1429注注: :此法亦稱為此法亦稱為預測預測- -校正法校正法 / /* * predictor-corrector m

22、ethod predictor-corrector method * */ /可以證明該算法可以證明該算法具有具有 2 階精度階精度，同時可以看到它，同時可以看到它是個是個單步單步遞推格式，比隱式公式的迭代求解過程遞推格式，比隱式公式的迭代求解過程簡單簡單。后面將看到，它的。后面將看到，它的穩定性高于穩定性高于顯式歐拉法。顯式歐拉法。改進的歐拉法改進的歐拉法111(,)(,)2nnnnnnhyyf xyf xy2021/8/1430在實際計算時，可將歐拉法與梯形法則相結合，在實際計算時，可將歐拉法與梯形法則相結合，計算公式為計算公式為 ,.2, 1, 0),(),(2),()(11)1(1)0

23、(1kyxfyxfhyyyxhfyyknnnnnknnnnn應用改進歐拉法應用改進歐拉法, ,如果序列如果序列 收斂收斂, ,)1(1)0(1 nnyy它的極限便滿足方程它的極限便滿足方程111(,)(,)2nnnnnnhyyf xyf xy2021/8/1431改進歐拉法的截斷誤差改進歐拉法的截斷誤差)(0)(311hyxynn 因此，改進歐拉法公式具有因此，改進歐拉法公式具有 2 2 階精度階精度2021/8/1432例例2:2: 用改進用改進Euler公式求解例公式求解例1中的初值問題，中的初值問題， 取步長取步長 。0.1h 解：解：對此初值問題采用改進對此初值問題采用改進EulerE

24、uler公式，公式， 其具體形式為其具體形式為 21( )211111111( ,)2(,)2()12piiiiiiicppiiiiiiipciiiyyhf x yyhx yyyhf xyyhxyyyy 計算結果列于下表：計算結果列于下表：例例1:1: 用歐拉公式求解初值問題用歐拉公式求解初值問題 2201.201yxyxy （）0,1,.,11i 01y 2021/8/1433iixiy 1piy 1ciyiiy xy改進的改進的Euler法法iiy xyEuler法法01234567891 01 11 20.00.10.20.30.40.50.60.70.80.91.01.11.21.00

25、00000.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.9523330.9100950.8571430.7975510.7350250.6725670.6123550.5557930.5036510.4562230.4134810.000000

26、0.0000990.0001730.0001850.0001150.0000340.0002330.0004460.0006430.0008030.0009190.0009900.0010230.0000000.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.0020532021/8/1434通過計算結果的比較可以看出，改進的通過計算結果的比較可以看出，改進的Euler方法方法的計算精度比的計算精度比Euler方法要高。方法要高。2021/8/1435歐拉法誤差概述)

27、()(2)O()(y2Euler式對于)O(h)(2Euler對于顯式則精度相對的越高。越大，誤差階，一個算法，局部截斷一般431321321hOxyhThxhTxyhTnnnnnn ：對于梯形公式：法隱：法說來2021/8/14366.3 龍格庫塔方法 對許多實際問題來說，歐拉公式與改進歐拉對許多實際問題來說，歐拉公式與改進歐拉公式精度還不能滿足要求，為此從另一個角度來分公式精度還不能滿足要求，為此從另一個角度來分析這兩個公式的特點，從而探索一條構造高精度方析這兩個公式的特點，從而探索一條構造高精度方法的途徑法的途徑. 2021/8/1437受改進的Euler方法啟發，更一般算式可設為112

28、1211()2(,)(,)(0,1,2,.)nnnnnnyykkkhf xykhf xh ykn123111()-()()nnnnnTy xyO hyy x適當選擇參數 ， ， ，使局部截斷誤差,這里仍假定。改進歐拉法改進歐拉法11 122121(,)(0,1,2,.)(,)nnnnnnyykkkhf xynkhf xh yk2021/8)()()(,)(,)(,)()nnnxnnnnynnyy xhy xfxyf xyfxyhO h)(),(),(),()()(),(),(),(:Taylor323122hOyxfyxfyxfhxyhhOyxfhkyxfhyxhfk

29、nnynnnnxnnnynnxnn展開式由二元函數11 122121 (,) (,)(0,1,2,.)nnnnnnnyykkkhf xyhykhf xh ykn2021/8/1439由于四個參數，三個方程，因此有一個自由參數，即解答不唯一。1122232()()()(,)(,)(,)()nnnxnnnnynnyy xhy xfxyf xyfxyhO hTaylor:與展式相比較得1212122122312()()()()()hnnnny xy xhy xyxO h2021/8/1440 這是改進的Euler方法。1211(1),1,22取可得此時算式為1121211()2(,) (,)nnnn

30、nnyykkkhf xykhf xh yk2021/8/1441-R K這是二階方法.121(2)0,1,2取可得此時算式為12121(,)11(,) 22nnnnnnyykkhf xykhf xh yk2021/8/1442 R-K這也是二階方法。12132(3),443取可得又有算式11212113)4(,)22(,)33nnnnnnyykkkhf xykhf xh yk（2021/8/1443三階龍格-庫塔方法三階龍格-庫塔方法是用三個值 k1, k2, k3 的線性組合112312123221331332112(,)(,)(,)nnnnnnnnyykkkkhf xykhf xhcccababykkhf xh ykb k 要使三階龍格-庫塔方法具有三階精度，必須使其局部截斷誤差為 O(h4)將 k1, k2, k3 代入 yn+1 的表達式中，在 (xn, yn) 處用二元泰勒公式展開，與 y(xn+1) 在 xn 處的泰勒展開式比較2021/8/1444 類似二階