微分方程數值解法

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

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

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

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

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

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

7、)Dx y axb y xyy x記記3 3、相關定義、相關定義微分方程數值解法(2) 一般構造方法：一般構造方法：4、 迭代格式的構造迭代格式的構造(1) 構造思想：構造思想：將連續(xù)的微分方程及初值條件離散為線性方程將連續(xù)的微分方程及初值條件離散為線性方程組加以求解。由于離散化的出發(fā)點不同，產生出各種不同的數組加以求解。由于離散化的出發(fā)點不同，產生出各種不同的數值方法?；痉椒ㄓ校河邢薏罘址ǎ〝抵滴⒎郑?、有限體積法值方法?；痉椒ㄓ校河邢薏罘址ǎ〝抵滴⒎郑⒂邢摅w積法（數值積分）、有限元法（函數插值）等等。（數值積分）、有限元法（函數插值）等等。 微分方程數值解法(3) 如何保證迭代公式的穩(wěn)

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

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

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

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

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

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

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

15、長 。 0.1h 解解: : 應用應用EulerEuler公式于題給初值問題的具體形式為：公式于題給初值問題的具體形式為： 2120,1,.,1101iiiiyyhx yiy 其中其中 。0.1ixi計算結果列于下表：計算結果列于下表： 微分方程數值解法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微分方程數值解法可用來檢驗近似解的準確程度?？捎脕頇z驗近似解的準確程度。 進行計算，數值解已達到了一定的精度。進行計算，數值解已達到了一定的精度。這個初值問題的準確解為這個初值問題的準確解為 ， 21 1y xx從上表最后一列，我們看

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

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

19、梯形公式即梯形公式具有具有2 階精度階精度，比歐拉方法有了進步。，比歐拉方法有了進步。但注意到該公式是但注意到該公式是隱式公式隱式公式，計算時不得不用到，計算時不得不用到迭代法，其迭代收斂性與歐拉公式相似。迭代法，其迭代收斂性與歐拉公式相似。微分方程數值解法中點歐拉公式中點歐拉公式 /* midpoint formula */中心差商近似導數中心差商近似導數hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假設假設 , 則可以導出則可以導出即中點公式具有即中點公式具有 2 階精度。階精度。)(),(

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

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

22、f xyf xy微分方程數值解法在實際計算時，可將歐拉法與梯形法則相結合，在實際計算時，可將歐拉法與梯形法則相結合，計算公式為計算公式為 ,.2, 1, 0),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyyknnnnnknnnnn應用改進歐拉法應用改進歐拉法, ,如果序列如果序列 收斂收斂, ,)1(1)0(1 nnyy它的極限便滿足方程它的極限便滿足方程111(,)(,)2nnnnnnhyyf xyf xy微分方程數值解法改進歐拉法的截斷誤差改進歐拉法的截斷誤差)(0)(311hyxynn 因此，改進歐拉法公式具有因此，改進歐拉法公式具有 2 2 階精度階精度微分方

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

24、y改進的改進的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微分方程數值解法通過計

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

27、從而探索一條構造高精度方析這兩個公式的特點，從而探索一條構造高精度方法的途徑法的途徑. 微分方程數值解法受改進的Euler方法啟發(fā)，更一般算式可設為1121211()2(,)(,)(0,1,2,.)nnnnnnyykkkhf xykhf xh ykn123111()-()()nnnnnTy xyO hyy x適當選擇參數 ， ， ，使局部截斷誤差,這里仍假定。改進歐拉法改進歐拉法11 122121(,)(0,1,2,.)(,)nnnnnnyykkkhf xynkhf xh yk微分方程數值解法 這是改進的Euler方法。1211(1),1,22取可得此時算式為1121211()2(,) (,)

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