第八章常微分方程的初值問題1ppt課件

1、第八章 常微分方程的初值問題常微分方程常微分方程ODE） ：只含有未知函數(shù)的導(dǎo)數(shù)的方程。：只含有未知函數(shù)的導(dǎo)數(shù)的方程。 ODE的階數(shù)：的階數(shù)： 最高階導(dǎo)數(shù)的階數(shù)最高階導(dǎo)數(shù)的階數(shù)ODE問題問題初值問題：初值問題：邊值問題：邊值問題：已知初始點處的條件已知初始點處的條件已知初始點和末端點處的條件已知初始點和末端點處的條件本章只討論初值問題本章只討論初值問題 一階一階ODE問題的形式問題的形式 00( )( , ( )()yxf x y xy xy 1、假如、假如 f 與與 y 無關(guān)，則計算變?yōu)榈跓o關(guān)，則計算變?yōu)榈?章數(shù)值積分章數(shù)值積分中討論的任一種直接積分方法應(yīng)用中討論的任一種直接積分方法應(yīng)用 初

2、始條件初始條件2、假如、假如 f 是是 y 的函數(shù)的函數(shù) ，積分過程將不同于前者。，積分過程將不同于前者。假設(shè)假設(shè) f 是是 y 的線性函數(shù)，如：的線性函數(shù)，如：f=ay+b 其中其中a,b是常數(shù)或是是常數(shù)或是 t 的函數(shù)，的函數(shù)，此時原方程稱為線性此時原方程稱為線性O(shè)DE假設(shè)假設(shè) f 不是線性函數(shù)，方程就稱為非線性不是線性函數(shù)，方程就稱為非線性O(shè)DE。一、求一、求ODEODE的解析解的解析解dsolvedsolve輸出變量列表輸出變量列表=dsolve(eq1,eq2, . , eqn, cond1,cond2, . , condn, v1,v2,vn)其中其中 eq1、eq2、.、eqn為

3、微分方程，為微分方程，cond1、cond2、.、condn為初值條件，為初值條件，v1,v2,vn 為自變量。為自變量。注注1： 微分方程中用微分方程中用 D 表示對表示對 自變量自變量 的導(dǎo)數(shù)，如：的導(dǎo)數(shù)，如：Dy y； D2y y； D3y y例例 ：求微分方程：求微分方程 的通解，并驗證。的通解，并驗證。22xdyxyxedx y=dsolve(Dy+2 y=dsolve(Dy+2* *x x* *y=xy=x* *exp(-x2),x)exp(-x2),x) 結(jié)果為結(jié)果為 y =(1/2y =(1/2* *x2+C1)x2+C1)* *exp(-x2)exp(-x2) syms x;

4、 diff(y)+2 syms x; diff(y)+2* *x x* *y - xy - x* *exp(-x2)exp(-x2) 結(jié)果為結(jié)果為 ans = 0ans = 0注注2：如果省略初值條件，則表示求通解；：如果省略初值條件，則表示求通解；注注3：如果省略自變量，則默認自變量為：如果省略自變量，則默認自變量為 t 例例 ： dsolve(Dy=2*x) %dy/dt = 2x 結(jié)果為結(jié)果為 ans = 2*x*t+C1例例 ：求微分方程：求微分方程 在初值條件在初值條件 下的下的特解，并畫出解函數(shù)的圖形。特解，并畫出解函數(shù)的圖形。0 xxyye y=dsolve(x*Dy+y-exp

5、(x)=0,y(1)=2*exp(1),x)結(jié)果為結(jié)果為 y = (exp(x)+exp(1)/x ezplot(y) %此時的此時的y已經(jīng)是符號變量，故不用已經(jīng)是符號變量，故不用ezplot(y) 12( )ye 3232(0)(0), (0)1,0,0d ydydyy ydtdtdt dsolve(D3y=-y,y(0)=1,Dy(0)=0,D2y(0)=0)結(jié)果為結(jié)果為ans = 1/3*exp(-t)+2/3*exp(1/2*t)*cos(1/2*3(1/2)*t) 例：例：注注4 4：解微分方程組時，如果所給的輸出個數(shù)與方程個：解微分方程組時，如果所給的輸出個數(shù)與方程個數(shù)相同，則方程

6、組的解按詞典順序輸出；如果只給一個數(shù)相同，則方程組的解按詞典順序輸出；如果只給一個輸出，則輸出的是一個包含解的結(jié)構(gòu)輸出，則輸出的是一個包含解的結(jié)構(gòu)structurestructure類類型的數(shù)據(jù)。型的數(shù)據(jù)。例：求微分方程組例：求微分方程組 在初值條件在初值條件 下的特解下的特解530tdxxyedtdyxydt x,y=dsolve(Dx+5*x+y=exp(t), . Dy-x-3*y=0, x(0)=1, y(0)=0, t)0010|ttxy x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0, . x(0)=1, y(0)=0, t)或或r=dsolve(Dx+

7、5*x+y=exp(t), Dy-x-3*y=0, . x(0)=1, y(0)=0, t)這里返回的這里返回的 r 是一個是一個 結(jié)構(gòu)類型結(jié)構(gòu)類型 的數(shù)據(jù)的數(shù)據(jù)r.x %查看解函數(shù)查看解函數(shù) x(t)r.y %查看解函數(shù)查看解函數(shù) y(t)dsolve的輸出個數(shù)只能為一個的輸出個數(shù)只能為一個 或或 與方程個數(shù)相等。與方程個數(shù)相等。例例 ：用：用dsolve求解著名的求解著名的Van der Pol方程方程 222( )( )( )1)( )0d y tdy tyty tdtdt syms mu; y=dsolve(D2y+mu*(y2-1)*Dy+y)結(jié)果：結(jié)果：Warning, canno

8、t find an explicit solutiony =&where(_a,diff(_b(_a),_a)*_b(_a)+mu*_b(_a)*_a2-mu*_b(_a)+_a = 0, _b(_a) = diff(y(t),t), _a = y(t), y(t) = _a, t = Int(1/_b(_a),_a)+C1)注注5：若找不到解析解，則返回其積分形式。：若找不到解析解，則返回其積分形式。只有很少一部分微分方程組能求出解析解。只有很少一部分微分方程組能求出解析解。大部分微分方程組只能利用數(shù)值方法求數(shù)值解。大部分微分方程組只能利用數(shù)值方法求數(shù)值解。 Euler歐拉方法是求解

9、一階歐拉方法是求解一階ODE的一種簡的一種簡便方法。盡管精度不高，但由于簡單，特別適用于便方法。盡管精度不高，但由于簡單，特別適用于快速編程求解。它有三種格式：快速編程求解。它有三種格式： 二、用二、用Euler Euler 方法求數(shù)值解方法求數(shù)值解（a向前向前Euler 法法（b改進的改進的Euler 法法（c向后向后Euler法法介紹這些方法是為了了解初值問題求解的基本思想。介紹這些方法是為了了解初值問題求解的基本思想。 一階一階ODE00( )( , ( )()yxf x y xy xy 對對 ( )( , ( )yxf x y x 兩邊從兩邊從x0 到到 x 積分得：積分得：00( )

10、( , ( )xxxxyx dxf x y x dx 00( )()( , ( )xxy xy xf x y x dx 00( )()( , ( )xxy xy xf x y x dx ( 積分方程）積分方程）00( )()( , ( )xxy xy xf x y x dx 1 1、向前、向前EulerEuler法法推導(dǎo)推導(dǎo)1：設(shè)節(jié)點為：設(shè)節(jié)點為1(,),0,1,2,nnnnyyhf xyn 向前向前Euler法：法：1()()()(, ()nnnnny xy xy xf xy xh 00( )( , ( ), ()yxf x y xy xy 用向前差分公式代替導(dǎo)數(shù)：用向前差分公式代替導(dǎo)數(shù)：0

11、,(0,1,2,)nxxnh n1()()(, (),0,1,2,nnnny xy xhf xy xn 此處，此處，y (xn)表示表示 xn 處的理論解，處的理論解，yn表示表示y (xn)的近似的近似解解 一階一階ODE00( )( , ( )()yxf x y xy xy 對對 ( )( , ( )yxf x y x 兩邊從兩邊從xn 到到 xn+1 積分得：積分得：11( )( , ( )nnnnxxxxyx dxf x y x dx 11()()( , ( )nnxnnxy xy xf x y x dx 推導(dǎo)推導(dǎo)2：11()()( , ( )nnxnnxy xy xf x y x d

12、x yn 近似代替近似代替y (xn)1(,)nnnnyyhf xy 用矩形代替右邊的積分用矩形代替右邊的積分1()()(, ()nnnny xy xhf xy x 向前向前Euler法法例例 求出求出2=0,2(0)1dyxydxyxy 解析解和數(shù)值解并畫解析解和數(shù)值解并畫圖比較圖比較先用先用dsolve求解析解求解析解y=dsolve(Dy=y+2*x/(y2),y(0)=1,x)結(jié)果為結(jié)果為 y =1/3*(-18-54*x+45*exp(3*x)(1/3)解析解：解析解： 1 3311854453/xyxe function x,y=Euler_bf(fun,x0,y0,xmax,h)

13、 %fun為顯式一階微分方程右端的函數(shù)為顯式一階微分方程右端的函數(shù) %x0,y0為初始條件為初始條件,滿足滿足y(x0)=y0 %xmax為為x的取值最大值的取值最大值 %h為將為將x等分后的步長等分后的步長N=(xmax-x0)/h+1; %N為總的節(jié)點數(shù)為總的節(jié)點數(shù)x(1)=x0;y(1)=y0;for n=1:N-1 x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun,x(n),y(n);end下面再用向前歐拉法數(shù)值求解，為此先編寫向前歐拉法下面再用向前歐拉法數(shù)值求解，為此先編寫向前歐拉法的程序的程序最后通過圖形比較用向前歐拉得到的數(shù)值解和解析解的最后通過圖形比

14、較用向前歐拉得到的數(shù)值解和解析解的誤差誤差 fun=inline(y+2*x/y2,x,y); x1,y1=Euler_bf(fun,0,1,2,0.1); x2,y2=Euler_bf(fun,0,1,2,0.05); x=0:0.01:2; y=1/3*(-18-54*x+45*exp(3*x).(1/3); plot(x1,y1,*,x2,y2,x,x,y) axis(0,2,0,9)向前歐拉方法截斷誤差為向前歐拉方法截斷誤差為211()()nny xyo h00( )( , ( )()yxf x y xy xy 對對 ( )( , ( )yxf x y x 兩邊從兩邊從xn 到到 xn

15、+1 積分得：積分得：2、改進的、改進的Euler法法11()()( , ( )nnxnnxy xy xf x y x dx yn 近似代替近似代替y (xn)111 (,)(,)2nnnnnnhyyf xyf xy 用梯形代替右邊的積分用梯形代替右邊的積分111()() (, ()(, ()2nnnnnnhy xy xf xy xf xy x 梯形法梯形法從從n=0開始計算，每步都要求解一個關(guān)于開始計算，每步都要求解一個關(guān)于yn+1的方程的方程一般是一個非線性方程），可用如下的迭代法計算：一般是一個非線性方程），可用如下的迭代法計算：111 (,)(,)2nnnnnnhyyf xyf xy

16、(0)1(1)( )111(,) (,)(,)2nnnnkknnnnnnyyhf xyhyyf xyf xy (0,1,2,)k 利用此算法，可得：利用此算法，可得： (0)(1)(2)( )(1)11111,kknnnnnyyyyy 利用利用(1)( )11kknnyy （ 為允許誤差來控制是否繼續(xù)為允許誤差來控制是否繼續(xù)迭代迭代 (0)1(1)( )111(,) (,)(,)2nnnnkknnnnnnyyhf xyhyyf xyf xy 迭代法太麻煩，實際上，當?shù)ㄌ闊?，實際上，當h取得很小時，只讓上式中取得很小時，只讓上式中的第二式迭代一次就可以，即的第二式迭代一次就可以，即(0)1

17、(0)111(,) (,)(,)2nnnnnnnnnnyyhf xyhyyf xyf xy (0,1,2,)k 改進的改進的Euler法也叫歐拉預(yù)估法也叫歐拉預(yù)估校正法）校正法）預(yù)估算式預(yù)估算式校正算式校正算式改進的改進的Euler法法=向前歐拉法向前歐拉法+梯形法梯形法向后向后Euler法依賴于向后差分近似，其格式為：法依賴于向后差分近似，其格式為： 111(,)nnnnyyhf yt 3 3、向后、向后EulerEuler法法 精度與向前歐拉法相同。假如精度與向前歐拉法相同。假如 f 為非線性函數(shù)，則為非線性函數(shù)，則與改進的與改進的Euler算法一樣，在每一步中需要采用迭代算法一樣，在每一

18、步中需要采用迭代法。該方法有兩個優(yōu)點：法。該方法有兩個優(yōu)點：（a絕對穩(wěn)定；絕對穩(wěn)定；（b如果解為正，則可保證數(shù)值計算結(jié)果也為正。如果解為正，則可保證數(shù)值計算結(jié)果也為正。三、龍格庫塔三、龍格庫塔Runge-kutta方法方法 Euler法的一個重要缺陷是精度太低。為了獲得法的一個重要缺陷是精度太低。為了獲得高精度，就要減小高精度，就要減小 h ，這不僅會增加計算時間，也，這不僅會增加計算時間，也會產(chǎn)生舍入誤差。會產(chǎn)生舍入誤差。1 1、二階、二階Runge-kuttaRunge-kutta方法方法11 122121(,)(,)nnnnnnyyRkR kkhf xykhf xah ybk121211

19、212RRaRbR其中11 122121(,)(,)nnnnnnyyRkR kkhf xykhf xah ybk121211212RRaRbR其中121=1,2RRab若取，則變?yōu)?1212111()2(,)(,)nnnnnnyykkkhf xykhf xyk或或1111(,)(,)(,)2nnnnnnnnnnyyhf xyhyyf xyf xy其實就是歐拉預(yù)估其實就是歐拉預(yù)估校正方法或改進的歐拉法）校正方法或改進的歐拉法）例例 2=0,2(0)1dyxydxyxy 用改進的歐拉法即二階用改進的歐拉法即二階龍格龍格-庫塔法數(shù)值求解并庫塔法數(shù)值求解并與向前歐拉法、解析解畫與向前歐拉法、解析解畫圖

20、比較圖比較前面已求出解析解：前面已求出解析解： 1 3311854453/xyxe function x,y=Euler_mend(fun,x0,y0,xmax,h) %fun為顯式一階微分方程右端的函數(shù)為顯式一階微分方程右端的函數(shù) %x0,y0為初始條件為初始條件,滿足滿足y(x0)=y0 %xmax為為x的取值最大值的取值最大值 %h為將為將x等分后的步長等分后的步長N=(xmax-x0)/h+1;%N為總的節(jié)點數(shù)為總的節(jié)點數(shù)x(1)=x0;y(1)=y0;for n=1:N-1 x(n+1)=x(n)+h; k1=h*feval(fun,x(n),y(n); k2=h*feval(fun

21、,x(n+1),y(n)+k1); y(n+1)=y(n)+1/2*(k1+k2);end先編寫改進的歐拉法的程序：先編寫改進的歐拉法的程序：再分別調(diào)用再分別調(diào)用Euler_bf.mEuler_bf.m和和Euler_mend.mEuler_mend.m求解：求解： fun=inline(y+2*x/y2,x,y); x1,y1=Euler_bf(fun,0,1,2,0.1); x2,y2=Euler_mend(fun,0,1,2,0.1); x=0:0.01:2; y=1/3*(-18-54*x+45*exp(3*x).(1/3); plot(x1,y1,*,x2,y2,s,x,y) axi

22、s(0,2,0,9)3 3、三階龍格庫塔方法、三階龍格庫塔方法常見的三階龍格庫塔方法的格式為：常見的三階龍格庫塔方法的格式為：1211/231211123(, )(/2,)(2,)1(4)6nnnnnnnnkhf y tkhf yktkhf ykk tyykkk二階龍格庫塔方法截斷誤差二階龍格庫塔方法截斷誤差為為311()()nny xyo h精度不高，希望通過增加計算精度不高，希望通過增加計算f的次數(shù)提高截斷誤差的的次數(shù)提高截斷誤差的階數(shù)，為此引入階數(shù)，為此引入4 4、四階龍格庫塔方法、四階龍格庫塔方法常見的四階龍格常見的四階龍格-庫塔方法有兩種，庫塔方法有兩種，一種基于一種基于1/3辛普森

23、法，格式：辛普森法，格式：1211/2321/243111234(, )(/2,)(/2,)(,)1(22)6nnnnnnnnnnkhf y tkhf yktkhf yktkhf yk tyykkkk另一種基于另一種基于3/8辛普森法，格式：辛普森法，格式：1211/33122/34123111234(, )(/3,)(/3/3,)(,)1(33)8nnnnnnnnnnkhf y tkhf yktkhf ykktkhf ykkk tyykkkk例例 2=0,2(0)1dyxydxyxy 用二階龍格用二階龍格-庫塔法和四階庫塔法和四階龍格龍格-庫塔法數(shù)值求解并與、庫塔法數(shù)值求解并與、解析解畫圖比

24、較解析解畫圖比較前面已求出解析解：前面已求出解析解： 1 3311854453/xyxe 先來編寫四階龍格先來編寫四階龍格- -庫塔法基于庫塔法基于1/31/3辛普森法的程序：辛普森法的程序：function x,y=RK4(fun,x0,y0,xmax,h) %fun為顯式一階微分方程右端的函數(shù)為顯式一階微分方程右端的函數(shù) %x0,y0為初始條件為初始條件,滿足滿足y(x0)=y0 %xmax為為x的取值最大值的取值最大值 %h為將為將x等分后的步長等分后的步長N=(xmax-x0)/h+1;%N為總的節(jié)點數(shù)為總的節(jié)點數(shù)x(1)=x0;y(1)=y0;for n=1:N-1 x(n+1)=x

25、(n)+h; k1=h*feval(fun,x(n),y(n); k2=h*feval(fun,x(n)+1/2*h,y(n)+k1/2); k3=h*feval(fun,x(n)+1/2*h,y(n)+k2/2); k4=h*feval(fun,x(n+1),y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end fun=inline(y+2*x/y2,x,y); x1,y1=Euler_mend(fun,0,1,2,0.1); x2,y2=RK4(fun,0,1,2,0.1); x=0:0.1:2; y=1/3*(-18-54*x+45*exp(3*

26、x).(1/3); plot(x1,y1,*,x2,y2,s,x,y) axis(0,2,0,9) 再分別調(diào)用再分別調(diào)用Euler_mend.mEuler_mend.m和和RK4.mRK4.m求解：求解：00.20.40.60.811.21.41.61.820123456789 二 階 RK四 階 RK解 析 解從圖形上看，似乎二階龍格從圖形上看，似乎二階龍格庫塔法與四階龍格庫塔法與四階龍格- -庫塔法庫塔法同樣接近解析解，故再從數(shù)值結(jié)果比較看看哪種誤差小同樣接近解析解，故再從數(shù)值結(jié)果比較看看哪種誤差小err=abs(y-y1),abs(y-y2) 結(jié)果為結(jié)果為 err = 0 0 0.000

27、9 0.0000 0.0016 0.0000 0.0022 0.0000 0.0026 0.0000 0.0031 0.0000 0.0035 0.0000 0.0040 0.0000 0.0047 0.0000 0.0054 0.0000 0.0062 0.0000 0.0073 0.0000 0.0084 0.0000 0.0098 0.0000 0.0115 0.0000 0.0133 0.0000 0.0155 0.0000 0.0181 0.0000 0.0210 0.0000 0.0243 0.0000 0.0281 0.0000四、二階ODE問題 二階二階ODE的一般形式為：的一

28、般形式為：00( )( )( )( ), (0)(0),u tau tbu ts t uu uu其中其中 是常數(shù)或是是常數(shù)或是 的函數(shù)，后兩個方程為的函數(shù)，后兩個方程為初初始條件。始條件。, ,a b s, ,t u u假如假如 與與u無關(guān)，則該方程稱為線性無關(guān)，則該方程稱為線性O(shè)DE;, ,a b s如果三者之中有一個是如果三者之中有一個是u或或 的函數(shù)，稱為非線性的。的函數(shù)，稱為非線性的。u 下面著重研究用向前下面著重研究用向前Euler法求解二階法求解二階ODE，及，及MATLAB程序。程序。 所以二階所以二階ODE分解為兩個一階分解為兩個一階ODE：設(shè)：設(shè)：uv0( )( )( )(

29、), (0)v tav tbu ts t vu12000( , , )( , , )(0)(0)uf u v tvvf u v tavbusuuvvu 00( )( )( )( )(0)(0),u tau tbu ts tuu uu200( , , )(0)vavbusf u v tvuv 定義定義12,fuvyfvfavbus 則上述兩個一階則上述兩個一階ODE寫為：寫為：( , )yf y t 其向前其向前Euler法的格式為：法的格式為：1(, )nnnnyyhf y t12000( , , )( , , )(0)(0)uf u v tvvf u v tavbusuuvvu 例例 求解二

30、階求解二階ODE5|200, (0)0,(0)1uuuuuu解：設(shè)解：設(shè)uv5200520vv vuvv vu 520uvvv vu 520vuv vuv 令令,5|20uvyfvv vu 則原方程的向量形式為：則原方程的向量形式為：( , )yf y t 0(0)0(0)(0)1uyyv 向前向前Euler法的格式為：法的格式為：1(, )nnnnyyhf y t0,5th = 0.05; t_max=5; n=1;y(:,1) = 0; 1;t(1) = 0;while t(n) t_max y(:,n+1) = y(:,n) + h*f_def(y(:,n),t); t(n+1) = t

31、(n)+h ; n=n+1;endaxis(0 5 -1 1);plot(t, y(1,:), t, y(2,:),:)function f=f_def(y,t)f=y(2);(-5*abs(y(2)*y(2)-20*y(1);先用函數(shù)文件定義先用函數(shù)文件定義f(u,v,t)再用向前再用向前Euler法求解法求解0123456-0.6-0.4-0.200.20.40.60.81 uv五、五、matlabmatlab相關(guān)命令數(shù)值求解相關(guān)命令數(shù)值求解ODEODEX,Y = 求解函數(shù)求解函數(shù)(fun,x0,xf,y0,option，p1,p2,.)其中：（其中：（1fun是用是用inline或函數(shù)文

32、件定義的顯式常微或函數(shù)文件定義的顯式常微分方程的函數(shù)名分方程的函數(shù)名函數(shù)文件格式：函數(shù)文件格式：或或inline格式：格式：function yd=函數(shù)名函數(shù)名(x,y,flag,p1,p2,)yd=.函數(shù)名函數(shù)名=inline(顯式方顯式方程程,x,y,flag,p1,p2,.)x為自變量，為自變量，y為因變量，為因變量，yd為為y的導(dǎo)數(shù)的導(dǎo)數(shù),flag用于控制求解過程，用于控制求解過程，p1，p2為方程中的參數(shù)為方程中的參數(shù)X,Y = 求解函數(shù)求解函數(shù)(fun,x0,xf,y0,option，p1,p2,.)其中：其中：（2）x0,xf為求解區(qū)間，為求解區(qū)間， y0 為初值條件；為初值條件

33、；（3option可省略為控制選項用于設(shè)置誤差限，可省略為控制選項用于設(shè)置誤差限，求解方程最大允許步長等等），求解方程最大允許步長等等），（4p1,p2為微分方程中的附加參數(shù)為微分方程中的附加參數(shù)（5Matlab在數(shù)值求解時自動對求解區(qū)間進行分割，在數(shù)值求解時自動對求解區(qū)間進行分割，X (向量向量) 中返回的是分割點的值中返回的是分割點的值(自變量自變量)，Y (向量向量) 中返回的是解函數(shù)在這些分割點上的函數(shù)值。中返回的是解函數(shù)在這些分割點上的函數(shù)值。（6求解函數(shù)可以是求解函數(shù)可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb求解函數(shù)求解函

34、數(shù) ODE類型類型特點特點說明說明ode45非剛性非剛性單步法；單步法；4，5 階階 R-K 方法；方法；累計截斷誤差為累計截斷誤差為 (x)3大部分場合的首選方法大部分場合的首選方法ode23非剛性非剛性單步法；單步法；2，3 階階 R-K 方法；方法；累計截斷誤差為累計截斷誤差為 (x)3使用于精度較低的情形使用于精度較低的情形ode113非剛性非剛性多步法；多步法；Adams算法；高低精算法；高低精度均可到度均可到 10-310-6計算時間比計算時間比 ode45 短短ode23t適度剛性適度剛性 采用梯形算法采用梯形算法適度剛性情形適度剛性情形ode15s剛性剛性多步法；多步法；Gea

35、rs 反向數(shù)值微反向數(shù)值微分；精度中等分；精度中等若若 ode45 失效時，可失效時，可嘗試使用嘗試使用ode23s剛性剛性單步法；單步法；2 階階Rosebrock 算算法；低精度法；低精度當精度較低時，計算時當精度較低時，計算時間比間比 ode15s 短短ode23tb剛性剛性梯形算法；低精度梯形算法；低精度當精度較低時，計算時當精度較低時，計算時間比間比ode15s短短沒有一種算法可以有效地解決所有的沒有一種算法可以有效地解決所有的 ODE 問題，因此問題，因此MATLAB 提供了多種提供了多種ODE求解函數(shù)，對于不同的求解函數(shù)，對于不同的ODE，可以，可以調(diào)用不同的求解函數(shù)。調(diào)用不同的

36、求解函數(shù)。 求初值問題求初值問題 的數(shù)值解，求解范的數(shù)值解，求解范圍為圍為 0,0.5222201( )dyyxxdxy 例例 先定義需要求解的顯式微分方程為一個函數(shù)先定義需要求解的顯式微分方程為一個函數(shù)function yd=fun1(x,y) yd=-2*y+2*x2+2*x最后在命令窗口調(diào)用函數(shù)求解最后在命令窗口調(diào)用函數(shù)求解x,y=ode23(fun1,0,0.5,1);fun2=inline(-2*y+2*x2+2*x,x,y); 求初值問題求初值問題 的數(shù)值解，求解范的數(shù)值解，求解范圍為圍為 0,0.5222201( )dyyxxdxy 例例 先定義需要求解的顯式微分方程為一個函數(shù)先

37、定義需要求解的顯式微分方程為一個函數(shù)在命令窗口用在命令窗口用inline定義定義最后在命令窗口調(diào)用函數(shù)求解最后在命令窗口調(diào)用函數(shù)求解x,y=ode23(fun2,0,0.5,1);如果需求解的問題是高階常微分方程，則需將其化為如果需求解的問題是高階常微分方程，則需將其化為一階常微分方程組，此時需用函數(shù)文件來定義該常微一階常微分方程組，此時需用函數(shù)文件來定義該常微分方程組。分方程組。122212112100 2 00 7/()( ). ,( ).dxdtxdxdtxxxxx 令令 ，則原方程可化為，則原方程可化為12,dyxy xdt 求解求解 Ver der Pol 初值問題初值問題222+1

38、000 2 00 7()( ). ,( ).d ydyyydtdtyy 例例 ： 前面演示過，該方程沒有解析解，該方程不是一階前面演示過，該方程沒有解析解，該方程不是一階顯式微分方程，故需要進行變換顯式微分方程，故需要進行變換先用函數(shù)文件定義一階顯式微分方程組先用函數(shù)文件定義一階顯式微分方程組function y=vdp_eq1(t,x,flag,mu)y=x(2); -mu*(x(1)2-1)*x(2)-x(1);l 再編寫腳本文件再編寫腳本文件 vdpl.m，在命令窗口直接運行該文件。，在命令窗口直接運行該文件。 x0=-0.2;-0.7;tf=20;mu=1;t1,y1=ode45(vd

39、p_eq1,0,tf,x0, ,mu);mu=2;t2,y2=ode45(vdp_eq1,0,tf,x0, ,mu);plot(t1,y1,t2,y2,:)figure;plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),:)法法1：vdp_eq2=inline(x(2);-mu*(x(1)2-1)*x(2)-x(1), t,x,flag,mu);x0=-0.2;-0.7;tf=20;mu=1;t1,y1=ode45(vdp_eq2,0,tf,x0, ,mu);mu=2;t2,y2=ode45(vdp_eq2,0,tf,x0, ,mu);plot(t1,y1,t2,y2,: