常微分方程數(shù)值解補充知識

1、2022-3-191 常微分方程數(shù)值解常微分方程數(shù)值解 1 1 Euler折線法折線法 1.1.Euler法法 2.2.改進改進Euler法法 3. Euler法的預估法的預估校正法校正法 2 2 RungeKutta法法 1.1.二級二級RKRK法法 法法 2.2.二級二級RKRK法法 法法 3.3.三級三階法三級三階法 2022-3-192對于常微分方程初值問題對于常微分方程初值問題則（則（1）在區(qū)間）在區(qū)間a,b上存在唯一解上存在唯一解 y=y(x). 如果如果f(x,y) 在在 a,b(-,+)上連續(xù)，且關于上連續(xù)，且關于 y 滿滿足足 Lipschtz 條件：條件：（1）|f (x,

2、y1) f(x,y2)|L|y1-y2| （2） 對于對于(1)在區(qū)間在區(qū)間a,b上的唯一解上的唯一解 y=y(x)，一般情況下很，一般情況下很難求出其解析解，因此只能通過數(shù)值解法求其近似解。難求出其解析解，因此只能通過數(shù)值解法求其近似解。 也就說，構造適當?shù)臄?shù)值方法，利用也就說，構造適當?shù)臄?shù)值方法，利用(1)求出求出 y=y(x) 在在節(jié)點節(jié)點 x1,x2, ,xn 處的近似函數(shù)值處的近似函數(shù)值 y1 , y2 , , yn 。常用方法主要有兩種：常用方法主要有兩種：Euler 折線折線法和法和 Rune-Kutta 法。法。 00)(),(yxyyxfyx a,b2022-3-1931 1

3、 歐拉折線法歐拉折線法一一Euler 法法xi=x0+ih , i=0,1,2, ，n 對于初值問題對于初值問題將區(qū)間將區(qū)間a,b n 等分，步長為等分，步長為 h=(b-a)/n ，得到，得到n+1個分點個分點已知已知 y=y(x) 在在 x0 處的函數(shù)值為處的函數(shù)值為 y0 ，為求出函數(shù)在，為求出函數(shù)在 x1 點點的函數(shù)值的函數(shù)值y(x1) ，先將方程，先將方程(1)進行轉化。進行轉化。 00)(),(yxyyxfyx a,b（1）2022-3-194在區(qū)間在區(qū)間x0 ,x1上將微分方程化為積分方程：上將微分方程化為積分方程：對于右端積分采用左矩形積分公式，得到近似積分：對于右端積分采用左

4、矩形積分公式，得到近似積分：這個近似值我們表示為：這個近似值我們表示為：1000(,)yyhf xy= =+ +即：即：000(,)yhf xy= =+ +x0 x1xyO 00)(),(yxyyxfdxdyx a,b 1010),(xxyydxyxfdy),()()(0001yxhfxyxy ),()()(0001yxhfxyxy 2022-3-1952111(,)yyhf xy= =+ +并稱該計算方法為并稱該計算方法為Euler折線性折線性。（3）依此類推可以求得函數(shù)依此類推可以求得函數(shù)y=y(x) 在所有分點在所有分點 x1,x2, ,xn 處的近似函數(shù)值處的近似函數(shù)值 y1 , y2

5、 , , yn ：采用同樣的方法，可以求出采用同樣的方法，可以求出y=y(x) 在在 x2 處的近似值處的近似值 y2 ： )(1, 2 , 1 , 0),(001xyyniyxhfyyiiii第第n次近似解的整體誤差為：次近似解的整體誤差為：),()()(1nnnnnnnyxhfyxyyxy 2022-3-196二、改進二、改進Euler法法 前面給出的前面給出的Euler折線性，由于采用的左矩形積分折線性，由于采用的左矩形積分公式，精度較低，如果我們采用梯形公式就可以加以改公式，精度較低，如果我們采用梯形公式就可以加以改進，提高計算精度。進，提高計算精度。對于下式的右端積分對于下式的右端積

6、分利用梯形公式得到：利用梯形公式得到：進而得到近似計算式：進而得到近似計算式： 100011(,)(,)2hyyf xyf xy= =+ + +依此類推可以推得一般的計算公式依此類推可以推得一般的計算公式 ： 1010),(xxyydxyxfdy ),(),(2)()(110001yxfyxfhxyxy 2022-3-197并稱其為并稱其為改進改進Euler法法，它是一個隱式計算格式。具體，它是一個隱式計算格式。具體計算時，需要從中解出計算時，需要從中解出 yi+1 來。來。（4）例例1 用用Euler法和改進法和改進Euler法計算初值問題法計算初值問題 1, 2 , 1 , 0),(),(

7、),(200111nixyyyxfyxfhyyiiiiii 1)0(1 . 0 , 0,219 . 0yxxyy2022-3-198解解:以以 h=0.02 為步長進行計算為步長進行計算,這時得區(qū)間這時得區(qū)間0,0.1上的分點上的分點由原方程由原方程xi = 0+ih=0.02i , i=0,1,2,3,4,5及及Euler折線公式折線公式得具體計算公式得具體計算公式 1)0(1 . 0 , 0,219 . 0yxxyy )(1, 2 , 1 , 0),(001xyyniyxhfyyiiii 14 , 3 , 2 , 1 , 0,21018. 001yixyyyiiii2022-3-199再由

8、原方程再由原方程0.9,0,0.112(0)1yyxxy 改進改進Euler折線公式折線公式得到得到這是一個隱式計算公式，但從中很容易解出這是一個隱式計算公式，但從中很容易解出 yi+1來：來： 1, 2 , 1 , 0),(),(),(200111nixyyyxfyxfhyyiiiiii 4 , 3 , 2 , 1 , 0, 12121009. 00111iyxyxyyyiiiiii2022-3-1910110.009112,0,1,2,3,40.009112iiiixyyix+ + +- -+ += = =+ + +y0 =1該初值問題的真解為該初值問題的真解為 y=(1+2x)-0.45

9、。用兩種算法計算出用兩種算法計算出5個點得近似值，再計算出精確解在個點得近似值，再計算出精確解在這些點的值，其結果列表如下：這些點的值，其結果列表如下：2022-3-1911ixiEuler解解 yj改進改進Euler解解 yj精確解精確解 yj001.000001.000001.0000010.020.982000.982500.9825120.040.965000.965950.9659630.060.948920.950260.9502840.080.933670.935370.9353950.100.919180.921200.92123表表5-1:三種解的比較三種解的比較從中可以看出

10、，改進從中可以看出，改進Euler法的結果要更精確一些法的結果要更精確一些。2022-3-1912三、三、Euler法的預估法的預估校正法校正法 在改進在改進Euler法中，有時并不容易解出法中，有時并不容易解出yi+1來，這時來，這時可以通過迭代法求解，得到如下的迭代公式：可以通過迭代法求解，得到如下的迭代公式：(0)1(,)iiiiyyhf xy+ += =+ +其中初值其中初值 通過通過Euler公式計算公式計算(0)1iy+ +合并起來就是如下的形式：合并起來就是如下的形式： ,2, 1 ,0),(),(2)(11)1(1kyxfyxfhyykiiiiiki ),(),(2),()(1

11、1)1(1)0(1kiiiiikiiiiiyxfyxfhyyyxhfyy2022-3-1913 用用Euler法提供初值，往往可以得到較好的結果，法提供初值，往往可以得到較好的結果，只需要迭代一次就可以求得很好的近似，因此上面的公只需要迭代一次就可以求得很好的近似，因此上面的公式可以改為如下的形式：式可以改為如下的形式：并稱其為并稱其為預估一校正法預估一校正法，其中，其中 稱為預估值，稱為預估值， yi+1為為校正值。校正值。1iy+ +如果進行編程計算，則改為下式：如果進行編程計算，則改為下式： ),(),(2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy 1, 1 ,0)

12、,(),()(2121211nihkyhxfkyxfkkkhyyiiiiii2022-3-1914例例2 用預估一校正法求解：用預估一校正法求解：取步長取步長 h=0.1 , xi=ih, i=0,1,2, ,10 。解：由公式解：由公式預估預估-校正校正計算公式計算公式 , 1)0(1 , 0,2yxyxyy 1, 1 , 0),(),()(2121211nihkyhxfkyxfkkkhyyiiiiii112()2iihyykk+ += =+ + +12iiixkyy= =- -2112()0,1,2,9iiixhkyhkiyhk+ += =+ +- -= =+ +L L2022-3-191

13、5010021xkyy= =- -= =100.1(10.9102)1.09592yy= =+ + += =21.1841y = =依此類推可以計算出依此類推可以計算出首先，由首先，由 y0=1，計算出，計算出 0201012()xhkyhkyhk+ += =+ +- -+ +101.7379y= =112()2iihyykk+ += =+ + +12iiixkyy= =- -2112()0,1,2,9iiixhkyhkiyhk+ += =+ +- -= =+ +L L9102. 011 . 011 . 0211 . 01 2022-3-19162 2 、RungeKutta 法法關于預估一校

14、正法，如果將其推廣為關于預估一校正法，如果將其推廣為1122()iiimmyyh k k k+ += =+ + + + +L L12221233311322112211(,)(,(,)(,)iiiiiimimimmmmmkf x ykf x h y hkkf x h y hk hkkf x h y hk hkhk- - - = = = =+ + + = =+ + + + = =+ + + + + L L L LL L則稱其為則稱其為m級級RungeKutta 法，其中法，其中,(1, 2,1, 2,1)iiijimjm= = =- -L LL L為常數(shù)，這些常數(shù)的選取，應該使得誤差盡可能的高。

15、為常數(shù)，這些常數(shù)的選取，應該使得誤差盡可能的高。2022-3-1917一、二級一、二級二階二階RK 法法11122121()(,)(2,)iiiiiiyyh k kkfxykfxhy h k+ + = =+ + + = = = =+ + + 其誤差為其誤差為11()()(,()iiiiiRy xhy x hfxy x+ += =+ +- - -2(2,()(,()iiii hfxh y xhfxy x- -+ + +由由Talor展式得展式得2()()()2!iixyhy xhy xhffff+ += =+ + + +322()3!xxxyyyyxyhff ff ff ff f+ + + +

16、+ + +L L2022-3-1918由二元函數(shù)的由二元函數(shù)的Taylor展式展式(,)( , )( , )( , )ryf ar bsf a brfa bsfa b+ + += =+ + +221( , )2( , )( , )2!xxxyyyr fa brsfa bs fa b+ + + + + L L(,(, ()iiiixyf xh yhf x y xfhfhff+ + += =+ + +222222122xxxyyy h ffh ff h f f+ + + + + L L21122211(1)()()22ixyRhfhaf ff f+ += =- - -+ +- -+ +- -322

17、2222211111()()()62362xxxyyyh ff ff ff +-+-+-+-+-+- 21()6yxyffff + + + L L得到得到帶入帶入 Ri+1 得到：得到：2022-3-1919 要使要使Ri+1 的階數(shù)盡可能的高，應選取的階數(shù)盡可能的高，應選取1 、2、 使使 h、h2 的系數(shù)為零，根據(jù)的系數(shù)為零，根據(jù) f 的任意性，應使的任意性，應使122210102102 - - -= = - -= = - -= = 122112 =-=- = 解得解得這時這時 R Ri+1i+1=O( =O( h h3 3 ) ) 有有 p=2 2 階精度。這樣得到的方法稱階精度。這樣得

18、到的方法稱為為二級二階二級二階R-KR-K法法。 2 可以任意取定可以任意取定.2022-3-19201). 當當2=1/2 時時1 =1/2 ，=1 112121()2(,)(,)iiiiiihyykkkfxykfxhyhk+ + = =+ + + = = = =+ + + 12121(,)(,)22iiiiiiyyh kkfxyhhkfxyk+ + = =+ + = = = =+ + + 則有則有為預估一校正法為預估一校正法2). 當當2=1 時時1 =0 ，=1/2 具體為具體為11(,(,)22iiiiiihyyhfxyhf xy+ += =+ + + +2022-3-1921二、三級三階二、三級三階R-K 法法1123121312()6(,)(,)22(,2)iiiiiiiihyykkkkf x yhhkf xykkf xh yhkhk+ + = =+ + + + = = = =+ + + = =+ +- -+ + 利用和二級二階同樣的推理方法可以得出相應的三級利用和二級二階同樣的推理方法可以得出相應的三級三階三階R-K法的計算公式法的計算公式2022-3-1922三、四級四階標準三、四級四階標準R-K 法法112341213243(22)6(,)(,)22(,)22(,iiiiiiiiiihyykkkkkf x yhhkf xykhhkf xy