工程計算基礎(chǔ)第9章-1

1、第10章 常微分方程的數(shù)值解法微分方程 微分方程是精確表示自然科學(xué)中各種基本定律和各種問題的基本工具之一?，F(xiàn)代建立起來的自然科學(xué)和社會科學(xué)中的數(shù)學(xué)模型大多都是微分方程。 在許多物理、力學(xué)、生物等現(xiàn)象中，不能直接找到聯(lián)系所研究的那些量的規(guī)律，但卻容易建立起這些量與它們的導(dǎo)數(shù)或微分間的關(guān)系。含有未知函數(shù)的導(dǎo)數(shù)（或微分）的方程，稱為微分方程。未知函數(shù)可以不出現(xiàn)，但其導(dǎo)數(shù)一定要出現(xiàn)。未知函數(shù)可以不出現(xiàn)，但其導(dǎo)數(shù)一定要出現(xiàn)。車輛縱向動力學(xué)模型mFa直流電機模型牛頓定律和基爾霍夫定律電機特性其中： T為電磁轉(zhuǎn)矩； i為電樞電流； Kt為轉(zhuǎn)矩系數(shù)； e為反電動勢； Ke為反電動勢系數(shù)。二自由度轉(zhuǎn)向模型（1）

2、忽略轉(zhuǎn)向系統(tǒng)的影響，直接以前輪轉(zhuǎn)角作為輸入；（2）忽略懸架的作用；車身只作平行于地面的平面運動，沿z 軸的位移、繞 y軸的俯仰角和繞 x 軸的側(cè)傾角均為零，且lrZZFF（3）汽車前進速度u視為不變；（4）側(cè)向加速度限定在0.4g一下，確保輪胎側(cè)偏特性處于線性范圍； （5）驅(qū)動力不大，不考慮地面切向力對輪胎側(cè)偏特性的影響，沒有空氣動力的作用。常微分方程 什么是常微分方程 ODE，Ordinary Differential Equations(1) y= kx, k 為常數(shù)；(2) ( y - 2xy) dx + x2 dy = 0；(3) mv(t) = mg - kv(t)；;112yay

3、(4).,(0sind22為常數(shù)lglgtd(5)常微分方程、偏微分方程、微分方程的階 未知函數(shù)是一元函數(shù)的微分方程稱做常微分方程， 未知函數(shù)是多元函數(shù)的微分方程稱做偏微分方程 微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)，稱為微分方程的階 n 階微分方程的一般形式為：F(x, y, y, , y(n) = 0，其中 x 是自變量， y 是未知函數(shù)，F(xiàn)(x, y, y, , y(n)是已知函數(shù)，222222222y=4xy -(y) +12xy=000ydydytydtdtTTTxyz 一階線性二階非線性一階非線性二階偏微分方程例如例如關(guān)于微分方程的解 任何代入微分方程后使其成為恒等式的函數(shù)，都叫

4、做該方程的解. 若微分方程的解中含有任意常數(shù)的個數(shù)與方程的階數(shù)相同，且任意常數(shù)之間不能合并，則稱此解為該方程的通解(或一般解). 當(dāng)通解中的各任意常數(shù)都取特定值時所得到的解，稱為方程的特解. 為了確定特解中的任意常數(shù)，對n階方程，相應(yīng)必須提出n個條件，這樣的條件稱為定解條件。微分方程與定解條件聯(lián)立構(gòu)成定解問題。微分方程的初值問題和邊值問題 給出方程定義區(qū)間左端點處的某些值（初始值），這樣的條件稱為初始條件。方程和初始條件構(gòu)成的問題稱為初值問題。一階微分方程的初始條件，如 二階微分方程的初始條件，如 給出方程解在定義區(qū)間兩端點處的某種值，稱為邊界條件，方程和邊界條件構(gòu)成邊值問題。00|yyxx,

5、|0000yyyyxxxx及bxayyxfy ),(;)(,)(byay結(jié)合下述三種邊界條件之一：;)(,)(byay;)()(,)()(1010bybyayay(3)式中。0, 0, 00000它們分別為第一、第二、第三邊界問題。方程：(1)(2)(3)n階初值問題的轉(zhuǎn)化高階微分方程高階微分方程( (或方程組或方程組) )的初值問題的初值問題, ,原則上都可以歸結(jié)為一階原則上都可以歸結(jié)為一階方程組來求解。例如方程組來求解。例如, ,有二階微分方程的初值問題有二階微分方程的初值問題 0000)(,)(),(yxyyxyyyxfy在引入新的變量在引入新的變量 后后, ,即化為一階方程組初值問題即

6、化為一階方程組初值問題: :yz0000)(,)(,),(yxzyxyzyzyxfz微分方程的解析解與數(shù)值解0(0)21yxyyxtxtxy0dee)(22解：但在許多工程實踐和科學(xué)研究中，除了少數(shù)特殊類型的微分方程能用解析法求得其精確解外，大多數(shù)情況要得出解的解析表達式是極其困難的，甚至是不可能的。常微分方程的數(shù)值問題：由于高階常微分方程可以轉(zhuǎn)化為一階常微分方程組，因此，為了不失一般性，本章主要介紹一類一階常微分方程初值問題：00)()(yxyyx,fy10.110.1 歐拉（歐拉（EulerEuler）法法 10.110.1.1 .1 歐拉法的基本思想歐拉法的基本思想 雖然歐拉法在實際應(yīng)用

7、時較少采用，但由于它的結(jié)構(gòu)簡單，在理論上具有非常重要的意義。歐拉法的基本思想就是用差分方程初值問題 (10.3) 的解來近似微分方程初值問題(10.2)的解，其中 ，式(10.3)也稱為歐拉公式。 )()0,1,()(01ayyiy,xfhyyiiii2)-(/abh 歐拉法的幾何意義是用一條自點 出發(fā)的折線去逼近積分曲線 ，如圖10.1所示。因此，這種方法又稱為折線法。顯然，歐拉法簡單地取折線的端點作為數(shù)值解，精度非常差。 )(00y,x)(xfy 10.110.1.2 .2 實現(xiàn)歐拉法的基本步驟實現(xiàn)歐拉法的基本步驟（1 1）輸入）輸入 的區(qū)間的區(qū)間 ，區(qū)間的等分數(shù)，區(qū)間的等分數(shù) 以及以及

8、在在 處處 的初始值的初始值 ；（2 2）對）對 ，計算：，計算： （3 3）輸出）輸出 ；（4 4）結(jié)束。）結(jié)束。xnx,x0Ny0 x0y10,1,2,Ni)(iiiy,xff iiifhyyhi*xxiiN/xxhn)-(0Ny例例10.10.1 1 用歐拉法求解初值問題用歐拉法求解初值問題 的數(shù)值解（保留的數(shù)值解（保留4 4位小數(shù)）。位小數(shù)）。02(0)y10.y,xxyy# #include include float float funcfunc(float (float x,floatx,float y) y) return(y-x); return(y-x); float fl

9、oat eulereuler(float x0,float (float x0,float xn,floatxn,float y0,int N) y0,int N) float float x,y,hx,y,h; ; intint i i; ; x=x0; x=x0; y=y0; y=y0; h=(xn-x0)/(float)N; / h=(xn-x0)/(float)N; /* * 計算步長計算步長 * */ / for(for(i i=1;i=1;i=N;iN;i+) /+) /* * 歐拉公式歐拉公式 * */ / y=y=y+hy+h* *funcfunc( (x,yx,y);); x

10、=x0+i x=x0+i* *h;h; return(y); return(y); main()main() float x0,xn,y0,e; float x0,xn,y0,e; intint N; N; clrscrclrscr();(); printfprintf(“(“ninputninput n:n ”); / n:n ”); / 輸入?yún)^(qū)間等分數(shù)輸入?yún)^(qū)間等分數(shù) scanfscanf(%(%d,&Nd,&N);); printfprintf(“input x0,xn:n ”); / (“input x0,xn:n ”); / 輸入輸入x x的區(qū)間的區(qū)間 scanfsca

11、nf(%f%f,&x0,&xn);(%f%f,&x0,&xn); printfprintf(input y0:n ); / (input y0:n ); / 輸入輸入x x0 0處的處的y y的值的值 scanfscanf(%f,&y0);(%f,&y0); e= e=eulereuler(x0,xn,y0,N);(x0,xn,y0,N); printfprintf(y(%f)=%6.4f,y0,e);(y(%f)=%6.4f,y0,e); 程序運行結(jié)果：程序運行結(jié)果：input n:input n:5 5input x0,xn:input x0

12、,xn: 0.0 1.0 0.0 1.0input y0:input y0: 2.0 2.0y(2.000000)=4.4883y(2.000000)=4.4883input n:input n:1010input x0,xn:input x0,xn: 0.0 1.0 0.0 1.0input y0:input y0: 2.0 2.0y(2.000000)=4.5937y(2.000000)=4.5937input n:input n:2020input x0,xn:input x0,xn: 0.0 1.0 0.0 1.0input y0:input y0: 2.0 2.0y(2.000000

13、)=4.6533y(2.000000)=4.65331e)(xxyx此問題的精確解為此問題的精確解為 ，相應(yīng)計算值為相應(yīng)計算值為4.71834.7183。從計算結(jié)果。從計算結(jié)果可以看出，雖然等分數(shù)越大，數(shù)值解可以看出，雖然等分數(shù)越大，數(shù)值解越精確，但精度并不理想。越精確，但精度并不理想。10.1.10.1.3 3 改進的歐拉法改進的歐拉法 10.1.10.1.3 3 改進的歐拉法的基本思想改進的歐拉法的基本思想 改進的歐拉法以隱函數(shù)的方式給出 ，計算時常使用迭代法，一般可先由歐拉公式計算出 的初始值 ，再按下面的格式進行 的迭代計算 (10.4) 改進的歐拉法由于添加了迭代過程，計算量較歐拉法

14、增加了一倍，付出這種代價的目的是為了提高精度。1iy1iy(0)1iy1iy)0,1,2,()()(2)()(111)(1(0)1ky,xfy,xfhyyy,xfhyykiiiiikiiiii10.1.410.1.4 實現(xiàn)改進的歐拉法的基本步驟實現(xiàn)改進的歐拉法的基本步驟（1 1）輸入）輸入 的區(qū)間的區(qū)間 ，區(qū)間的等分數(shù)，區(qū)間的等分數(shù) 以及以及 在在 處處 的初始值的初始值 ；（2 2） , , ， ；（3 3） ；（4 4）計算：）計算：（5 5）若）若 ， ，轉(zhuǎn)（，轉(zhuǎn)（4 4）；否則，轉(zhuǎn)（）；否則，轉(zhuǎn)（6 6）；）；（6 6）輸出）輸出 ；（7 7）結(jié)束。）結(jié)束。xnx,x0Ny0 x0y/

15、N)x-x(h0n0 xx 0yy 1i)yx,(fhyyphi*xx0)(pcyx,fh*yy2(/ )yyycpNi 1 iiy# #include include float float funcfunc(float (float x,floatx,float y) y) return(y-x); return(y-x); float float eulereuler(float x0,float (float x0,float xn,floatxn,float y0,int N) y0,int N) float float x,y,yp,yc,hx,y,yp,yc,h; ; intin

16、t i i; ; x=x0; x=x0; y=y0; y=y0; h=(xn-x0)/(float)N; / h=(xn-x0)/(float)N; /* * 計算步長計算步長 * */ / for(for(i i=1;i=1;i=N;iN;i+) /+) /* * 改進的歐拉公式改進的歐拉公式 * */ / ypyp= =y+hy+h* *funcfunc( (x,yx,y);); x=x0+i x=x0+i* *h;h; ycyc= =y+hy+h* *funcfunc( (x,ypx,yp);); y=( y=(yp+ycyp+yc)/2.0;)/2.0; return(y); retu

17、rn(y); 例例10.10.1 1 用改進的歐拉法求解初值問題用改進的歐拉法求解初值問題 的數(shù)值解（保留的數(shù)值解（保留4 4位小數(shù)）。位小數(shù)）。02(0)y10.y,xxyy程序運行結(jié)果：程序運行結(jié)果：input n:input n:5 5input x0,xn:input x0,xn: 0.0 1.0 0.0 1.0input y0:input y0: 2.0 2.0y(2.000000)=4.7027y(2.000000)=4.7027input n:input n:1010input x0,xn:input x0,xn: 0.0 1.0 0.0 1.0input y0:input y0

18、: 2.0 2.0y(2.000000)=4.7141y(2.000000)=4.7141input n:input n:2020input x0,xn:input x0,xn: 0.0 1.0 0.0 1.0input y0:input y0: 2.0 2.0y(2.000000)=4.7172y(2.000000)=4.717210.1.5 歐拉方法的收斂性 局部截斷誤差由歐拉公式計算一步得到的誤差記為： 11212 () Taylor, (,),()()()()()2 ()(, ()() (1)2nnnnnnnnnnnnnny xxxxhy xy xhy xhy xyhy xhf xy

19、xy將在點展開1(),(8.2) () (, () (2)nnnnny xyy xh f xy x用準(zhǔn)確值 利用歐拉公式計算得到的稱為“局部截斷誤差”。令 定義 若給定方法的局部截斷誤差滿足則稱該方法是 p 階的，或稱為具有 p 階精度。11112() ()()(, () (). 2nnnnnnnnTy xyy xy xhf xy xhy2221 max |( )|,|()|(). 22a x bnnMy xhhTyMO h 則11|(),pnTO h 整體截斷誤差 記111().nnney xy112121112 ,(), (), () ,(),.nnnnnnyy yyy xy xy xy

20、xey yy因為計算 時 用到的, , 是 的近似值 每步產(chǎn)生的誤差會累積到計算的誤差中因此 與 , , 都有關(guān) 稱為整體截斷誤差1111111111| | ()| | ()| |. (3)nnnnnnnnnney xyy xyyyTyy1111,(8.2)| | () (, () (,) | | () |(, ()(,)| | () | () | () (1)nnnnnnnnnnnnnnnnnnnnyyyyy xh f xy xyh f xyy xyh f xy xf xyy xyhL y xyhL對應(yīng)用歐拉公式得李普希茲條件|. (4)ne由此知，當(dāng) max |,(4)(3)kkTT記 將

21、代入得1121210112|(1)| (1)(1)| (1)(1) | (1)(1)(1) (1)|(1)1(1)1 ()(1) 1 (nnnnnnnneThLeThL ThLeThL ThLeThL ThL ThL ThLehLhLTO hhLhLO).h10,|0, nhe有歐拉公式是一階收斂的。 注 整體截斷誤差與局部截斷誤差的關(guān)系： 類似于向前歐拉公式收斂性的分析，可以得到： 向后歐拉公式： 局部截斷誤差： 整體截斷誤差： 收斂的條件： 1(1)() 0,(1)e , (1)ee.nnznnhLL b azzhLconst111|(|).nneO hT21|(),nTO h1|( ).

22、neO h01,8.4hL因為由（）知，.(1)( )( )(1)( )(1)11111111|(,)(,)|.Lkkkkkknnnnnnnnyyh f xyf xyhL yy 梯形公式：（二階方法） 局部截斷誤差： 整體截斷誤差： 收斂的條件： 31|();nTO h21|().neO h101,8.92hL因為由（） 知，(1)( )( )(1)111111.( )(1)11|(,)(,) |.2kkkknnnnnnLkknnyyh f xyf xyhLyy10.2 龍格龍格- -庫塔庫塔(Runge-Kutta)方法方法 回顧： 向前Euler: 向后Euler: 梯形公式: 研究提高收

23、斂階的方法 用Taylor展式|( );neO h|( );neO h2|().neO h21( )1()()()()()2! ()(),!nnnnnpppnhy xy xhy xhy xy xhyxO hp采用此式，需要計算高階導(dǎo)數(shù)，如112( )1 (), .2!nnppnnnnny xyhhyyhyyyp22( , ),( , ),()2, xyxyxyxxxyxyyyyyf x yyfx yffyfffyfffffff ffff f 10.2.1 二階二階Runge-Kutta)方法方法 基本思想：分別計算不同點上的函數(shù)值，然后對這些函數(shù)值進行線性組合，構(gòu)造出近似計算公式，然后把此近似

24、公式與解 的Taylor展開式進行比較，使他們前面的若干項保持一致，從而使近似公式達到需要的收斂階。 考察梯形公式： 引入?yún)?shù) k1， k2 ，上式改寫為1()ny x111 (,)(,) .2nnnnnnhyyf xyf xy 為便于后面的推廣，將上式寫成更一般的形式 任務(wù)：適當(dāng)選取參數(shù) 使公式（2）的收斂階盡可能的高。11212111 2(,) (1) (,)(,) nnnnnnnnhyykkkf xykf xyf xh yhk11 122122211 (,) (2) (,) nnnnnnyyh c kc kkf xykf xa h yb hk12221,c c a b 先將k2在(xn,

25、 yn)處Taylor展開 將k1 ，k2 代入（2)中的第一式，得23123()()()()()2! ()(, ()(, ()2 (, ()(, ()() (4)nnnnnnnxnnynnnnhy xy xhy xy xO hhy xhf xy xfxy xfxy xf xy xO h而 22211222112221(,) (,)(,)(,)() (,)(,)(,)(,)()nnnnxnnynnnnxnnnnynnkf xa h yb hkf xya hfxyb hk fxyO hf xya hfxyb hf xyfxyO h2112222321 2(,)(,) (,)(,)() (3)nn

26、nnxnnnnynnyyh ccf xya c h fxyb c h f xyfxyO h要使(3)逼近(4)具有O(h3) 的截斷誤差,比較兩式的系數(shù)得 特別地， 取12222211 2121cca ca b122211, 1,(2)2ccab 由得梯形公式：112121() 2(,) (8.15) (,) nnnnnnhyykkkf xykf xh yhk 取12121 (,) (8.16) (,)22nnnnnnyyhkkf xyhhkf xyk12221101, ,(2)2ccab，由即得中點公式：(5),(6)-也稱為二階龍格 庫塔公式.2（計算 個點函數(shù)值）3121|(), |()

27、.nnTO heO h 10.2.2 三階三階Runge-Kutta方法方法 增加一個參與運算的函數(shù)值，設(shè) 將k2 , k3在(xn, yn)處Taylor展開到三階，然后將k2 , k2 , k3代入（5）式，合并同類項，然后與 y(xn+1) 在 xn 處Taylor展開對應(yīng)項相比較，得參數(shù)滿足的方程：11 1223312221133311322 (, ) (5)(, ) (, ) nnnnnnnnyyh c kc kc kkf xykf xa h yb hkkf xa h yb hkb hk12322133132223322223332 321 1 21 31 6cccababbc ac

28、 ac ac ac a b (6) 取滿足（6）式，代入（5）式，得三階R-K公式（10.17）：12322133132121, ,6361, 1, 1, 2,2cccababb ，11231213124 6(, ) (8.17)11(, )22 (, 2) nnnnnnnnhyykkkkf xykf xh yhkkf xh yhkhk 取滿足（6）式，代入（5）式，得三階R-K公式（10.18）：12322133132130, ,44122, , 0, ,333cccababb，113121323 4(, ) (8.18)11(, )3322 (, ) 33nnnnnnnnhyykkkf x

29、ykf xh yhkkf xh yhk4131|()|()nnTO heO h 取滿足（6）式，代入（5）式，得三階R-K公式（10.19）：12322133132214, ,939133, , 0, ,244cccababb，112312132234 9(, ) (8.19)11(, )2233 (, ) 44nnnnnnnnhyykkkkf xykf xh yhkkf xh yhk10.2.3 四階四階Runge-Kutta方法方法11234121324322 6(, ) 11(, ) 2211(, ) 22 (, ) nnnnnnnnnnhyykkkkkf xykf xh yhkkf x

30、h yhkkf xh yhk(8.20) 局部截斷誤和整體截斷誤差：11234121312412333 8(, ) 11(, ) 3321(, ) 33(, ) nnnnnnnnnnhyykkkkkf xykf xh yhkkf xh yhkhkkf xh yhkhkhk(8.21)5141|()|()nnTO heO h10.10.4.2 4.2 實現(xiàn)標(biāo)準(zhǔn)四階龍格實現(xiàn)標(biāo)準(zhǔn)四階龍格庫塔法的基本步驟庫塔法的基本步驟（1 1）輸入）輸入 的區(qū)間的區(qū)間 ，區(qū)間的等分數(shù)，區(qū)間的等分數(shù) 以及以及 在在 處的初始值處的初始值 ；（2 2） , , ， ；（3 3） ；（4 4）計算：）計算：（5 5）若）

31、若 ， ，轉(zhuǎn)（，轉(zhuǎn)（4 4）；否則，轉(zhuǎn)（）；否則，轉(zhuǎn)（6 6）；）；（6 6）輸出）輸出 ；（7 7）結(jié)束。）結(jié)束。xnx,x0Ny0 x0y/Nxxhn)-(00 xx 0yy 1i2/hxxh)(1yx,fd )21(2/dh*y,xfdh)22(3/dh*y,xfdh)3(4dh*y,xfdh4)/63*22*21(ddddh*yyhi*xx0Ni 1 iiy例例10.10.3 3 用標(biāo)準(zhǔn)四階龍格用標(biāo)準(zhǔn)四階龍格庫塔法求解初值問題的數(shù)值庫塔法求解初值問題的數(shù)值 解（保留解（保留6 6位小數(shù)）。位小數(shù)）。0(0)10yxyxy#include #include float float fun

32、cfunc(float (float x,floatx,float y) y) return(x-y); return(x-y); float float runge_kuttarunge_kutta(float x0,float (float x0,float xn,floatxn,float y0,int N) y0,int N) float x,y,y1,y2,h,xh; float x,y,y1,y2,h,xh; float d1,d2,d3,d4; float d1,d2,d3,d4; intint i i; ; x=x0; x=x0; y=y0; y=y0; h=(xn-x0)/(

33、float)N; / h=(xn-x0)/(float)N; /* * 計算步長計算步長 * */ / for(for(i i=1;i=1;i=N;iN;i+) /+) /* * 標(biāo)準(zhǔn)四階龍格標(biāo)準(zhǔn)四階龍格- -庫塔公式庫塔公式 * */ / xhxh= =x+hx+h/2;/2; d1= d1=funcfunc( (x,yx,y);); d2= d2=funcfunc( (xh,y+hxh,y+h* *d1/2.0);d1/2.0); d3= d3=funcfunc( (xh,y+hxh,y+h* *d2/2.0);d2/2.0); d4= d4=funcfunc( (xh,y+hxh,y+h* *d3);d3); y= y=y+hy+h* *(d1+2(d1+2* *d2+2d2+2* *d3+d4)/6.0;d3+d4)/6.0; x=x0+i x=x0+i* *h; h; return(y); return(y); main()main() float x0,xn,y0,e; float x0,xn,y0,e; intint N; N; clrscrclrscr();(); printfprintf(ninputninput n:n ); / n:n ); /輸入?yún)^(qū)間等分數(shù)輸入?yún)^(qū)間等分數(shù) scanfscanf(%(%