MATLAB解常微分方程及其應(yīng)用

1、MATLAB解常微分方程及其應(yīng)用授課教師：授課教師：唐躍龍?zhí)栖S龍TelephoneQ: 794218136Email:1 微分方程的基本概念微分方程的基本概念定義1 含有未知函數(shù)的導(dǎo)數(shù)(或微分)，同時也可能含有未知函數(shù)與自變量的方程，稱為微分方程，未知函數(shù)是一元函數(shù)的，稱為常微分方程；未知函數(shù)是多元函數(shù)的，稱為偏微分方程.定義2 微分方程中所出現(xiàn)的未知函數(shù)的最高階導(dǎo)(或微分)的階數(shù)，叫做微分方程的階.定義3 若 是 的一次有理整式， 稱為 階線性常微分方程.( )( , ,)nF x y yy ( ),ny yy ( )( , ,)0nF x y yyn定義4 不是線

2、性的微分方程稱為非線性常微分方程.定義8 若某函數(shù)代入微分方程能使該方程成為恒等式，該函數(shù)就叫做該微分方程的解.定義5 不含有未知函數(shù)及其導(dǎo)數(shù)(或微分)的項稱為微分方程自由項.定義6 自由項為零的微分方程稱為齊次方程.定義7 自由項不為零的微分方程稱為非齊次方程.例1. 指出下列微分方程的階數(shù)并判斷是線性還是非線性、齊次還是非齊次方程：54(1)0;xyy y x y 4212(2)(7 6)0;x yxy 4(3 )sin0;xy y xyy 2(4)60.xxyy 定義5 確定了通解中任意常數(shù)以后得到的解稱為微分方程的特解.定義4 若微分方程的解中所含相互獨立的任意常數(shù)的個數(shù)與微分方程的階

3、數(shù)相同，這樣的解叫做微分方程的通解.例2. 驗證函數(shù) 是微分方程 的解.255xye210yy 例3. 求方程 的通解.6yx 例4. 求方程 滿足初始條件：當 時 的特解.23yx 1x 2y 注：這里所求的微分方程的解、通解、特解都是解析解（精確解）.2 MATLAB求微分方程的解析解求微分方程的解析解格式一：y=dsolve(f1, f2, , fm)格式二：y=dsolve(f1, f2, , fm ,x) %指明自變量 fi既可以描述微分方程，又可描述初始條件或邊界條件。如： 描述微分方程時 描述條件時 . 747)()4(yDty. 3)2(23)2( yDy例5 求方程 的解析解

4、.解： syms x,t;x=dsolve(D1x=x*(1-x2),t)x = 1 0 -1 (-1/(exp(C3 - 2*t) - 1)(1/2)(1 (*)()(2txtxtx若方程變?yōu)椋航猓?syms x,t; x=dsolve(Dx=x*(1-x2)+1)警告: Explicit solution could not be found; implicit solution returned. 因此，不是所有方程都有解析解，只有極少部分的微分方程有解析解.當找不到原方程的解析解時，我們只能考慮其數(shù)值解.1)(1 (*)()(2txtxtx以下面的一階常微分方程以下面的一階常微分方程(

5、ODE )初值問題為例初值問題為例 ：000( , ) ()ndyf x yxxxdxy xy()nnyy x 數(shù)值解法就是求數(shù)值解法就是求y(x)在某些分立的節(jié)點在某些分立的節(jié)點 xn 上的近似值上的近似值yn，用以近似用以近似y(xn)3.1 簡單歐拉簡單歐拉(L.Euler, 1707-1783)方法。方法。00( , )()dyf y xdxy xy012nyyyy 歐拉數(shù)值算法就是由初值通過歐拉數(shù)值算法就是由初值通過遞推求解遞推求解，遞推求解，遞推求解就是從初值開始，后一個函數(shù)值由前一個函數(shù)值得到。關(guān)就是從初值開始，后一個函數(shù)值由前一個函數(shù)值得到。關(guān)鍵是構(gòu)造遞推公式。鍵是構(gòu)造遞推公式

6、。3 歐拉近似方法歐拉近似方法yyxxdd歐拉數(shù)值算法遞推公式構(gòu)造歐拉數(shù)值算法遞推公式構(gòu)造3.1.1 差分法差分法差分法就是用差商近似代替微商，即差分法就是用差商近似代替微商，即代入微分方程得到：代入微分方程得到：()()( )( ( ),)( ( ), )yy xxy xxxy xf y xy x xxxxf對于等間隔節(jié)點對于等間隔節(jié)點11 n=0,1,2nnnnxxxhxxh 可以得到：可以得到：xn x0 x1 x2 . xn .y精確值精確值y(x0) y(x1) y(x2) . y(xn) .y近似值近似值y0 y1 y2 . yn .在在xn節(jié)點上，微分方程可以寫為節(jié)點上，微分方程

7、可以寫為1()()() , nnnny xy xf y xxh作如下近似：作如下近似：( )nnyy t則得到歐拉解法遞推公式的一般形式：則得到歐拉解法遞推公式的一般形式：1( , )nnnnyyf yxh具體求解過程為：具體求解過程為：1000( , x )yyf yh2111( , x )yyf yh3222( , x ) yyf yh簡單歐拉方法程序簡單歐拉方法程序function outx,outy=MyEuler(fun,x0,xt,y0,PointNum)%MyEuler 用前向差分的歐拉方法解微分方程用前向差分的歐拉方法解微分方程%fun 表示表示f（x，y）%x0,xt表示自變

8、量的初值和終值表示自變量的初值和終值%y0表示函數(shù)在表示函數(shù)在x0處的值處的值,其可以為向量形式其可以為向量形式%PointNum表示自變量在表示自變量在x0,xt上取的點數(shù)上取的點數(shù)if nargin5 | PointNum=0 %如果函數(shù)僅輸入如果函數(shù)僅輸入4個參數(shù)值，則個參數(shù)值，則PointNum默認值為默認值為100 PointNum=100;endif nargin dsolve(Dy+3*x*y=x*exp(-x2) ans =(1/3*exp(-x*(x-3*t)+C1)*exp(-3*x*t) dsolve(Dy+3*x*y=x*exp(-x2),x) ans = exp(-x

9、2)+exp(-3/2*x2)*C1例題：用例題：用MATLAB的符號解法的符號解法,求解常微分方程的特解：求解常微分方程的特解：1 20,2xxxyyeye dsolve(x*Dy+2*y-exp(x)=0,y(1)=2*exp(1),x) ans = (exp(x)*x-exp(x)+2*exp(1)/x2例題：采用例題：采用ODE45求解描述某剛性問題的微分方程：求解描述某剛性問題的微分方程：123121323123(0)0,(0)10.51(0)1yy yyyy yyyy yy function dy = rigid(t,y)dy = zeros(3,1); % 行向量行向量dy(1)

10、 = y(2) * y(3);dy(2) = -y(1) * y(3);dy(3) = -0.51 * y(1) * y(2);options = odeset(RelTol,1e-4,AbsTol,1e-4 1e-4 1e-5);T,Y = ode45(rigid,0 12,0 1 1,options);plot(T,Y(:,1),-,T,Y(:,2),-.,T,Y(:,3),.)legend(y1,y2,y3)例題：用例題：用MATLAB函數(shù)函數(shù)ode23,ode45,ode113,求解多階常微分方程求解多階常微分方程：323232233(1)1,(1)10, (1)30,1,10 xd

11、yd ydyxedxdxdxyyyx12232332333233xdyydxdyydxdyeyydxxxx1223333010000103230 xydYdyYdxdxyexxxfunction dy=myfun03(x,y)dy=zeros(3,1) %初始化變量初始化變量dydy(1)=y(2); %dy(1)表示表示y的一階導(dǎo)數(shù)的一階導(dǎo)數(shù),其等于其等于y的第二列值的第二列值dy(2)=y(3); %dy(2)表示表示y的二階導(dǎo)數(shù)的二階導(dǎo)數(shù)dy(3)=2*y(3)/x3+3*y(2)/x3+3*exp(2*x)/x3 %dy(3)表示表示y的三階導(dǎo)數(shù)的三階導(dǎo)數(shù)% 用用ode23 ode45

12、 ode113解多階微分方程解多階微分方程clear,clcx23,y23=ode23(myfun03,1,10,1 10 30);x45,y45=ode45(myfun03,1,10,1 10 30);x113,y113=ode113(myfun03,1,10,1 10 30);figure(1) %第一幅圖第一幅圖plot(x23,y23(:,1),*r,x45,y45(:,1),ob,x113,y113(:,1),+g) %作出各種函數(shù)所得結(jié)果作出各種函數(shù)所得結(jié)果legend(ode23解解,ode45解解,ode113解解)title(ODE函數(shù)求解結(jié)果函數(shù)求解結(jié)果)figure(2)

13、plot(x45,y45) %以以ode45為例作出函數(shù)以及其各階導(dǎo)數(shù)圖為例作出函數(shù)以及其各階導(dǎo)數(shù)圖legend(y,y一階導(dǎo)數(shù)一階導(dǎo)數(shù),y兩階導(dǎo)數(shù)兩階導(dǎo)數(shù))title(y,y一階導(dǎo)數(shù)一階導(dǎo)數(shù),y二階導(dǎo)數(shù)函數(shù)圖二階導(dǎo)數(shù)函數(shù)圖)例題：MATLAB在導(dǎo)熱計算傳熱過程涉及面廣，數(shù)學(xué)模型復(fù)雜。計算過程中涉及到傳熱過程涉及面廣，數(shù)學(xué)模型復(fù)雜。計算過程中涉及到許多運算方法。導(dǎo)熱雖是容易做數(shù)學(xué)處理的一種熱量傳遞許多運算方法。導(dǎo)熱雖是容易做數(shù)學(xué)處理的一種熱量傳遞方式，但其過程往往涉及常微分、偏微分方程、線性（非方式，但其過程往往涉及常微分、偏微分方程、線性（非線性）方程組的求解。線性）方程組的求解。有一外徑為

14、有一外徑為4cm，內(nèi)徑為，內(nèi)徑為1.5cm，載有電流密度，載有電流密度I為為5000A/cm2的內(nèi)冷鋼質(zhì)導(dǎo)體。導(dǎo)體單位時間發(fā)出的熱量等于流體同時帶走的內(nèi)冷鋼質(zhì)導(dǎo)體。導(dǎo)體單位時間發(fā)出的熱量等于流體同時帶走的熱量，導(dǎo)體內(nèi)壁面的溫度為的熱量，導(dǎo)體內(nèi)壁面的溫度為70。假定外壁面完全絕熱。試。假定外壁面完全絕熱。試確定，導(dǎo)體內(nèi)部的溫度分布；已知鋼的導(dǎo)熱系數(shù)確定，導(dǎo)體內(nèi)部的溫度分布；已知鋼的導(dǎo)熱系數(shù)k=0.38Kw/mK，電導(dǎo)率，電導(dǎo)率=21011 km。解：這是圓柱坐標中常物性一維穩(wěn)態(tài)導(dǎo)熱問題，結(jié)合本題具解：這是圓柱坐標中常物性一維穩(wěn)態(tài)導(dǎo)熱問題，結(jié)合本題具體條件導(dǎo)熱微分方程式可簡化為：體條件導(dǎo)熱微分方程

15、式可簡化為：結(jié)合邊界條件，可得這一導(dǎo)熱問題的數(shù)學(xué)描述為：結(jié)合邊界條件，可得這一導(dǎo)熱問題的數(shù)學(xué)描述為：此常微分方程的分析解，可調(diào)用此常微分方程的分析解，可調(diào)用MATLAB符號工具箱中的符號工具箱中的dsolve函數(shù)來實現(xiàn)。在命令窗口執(zhí)行下面的代碼：函數(shù)來實現(xiàn)。在命令窗口執(zhí)行下面的代碼：cleard_equat=D2t+Dt/r+131579=0;condition=（t0.0075）=70,D（t0.02）=0;%邊界條件邊界條件t=dsolve（d_equat,condition,r）程序執(zhí)行結(jié)果：程序執(zhí)行結(jié)果：程序執(zhí)行結(jié)果：程序執(zhí)行結(jié)果：即求出溫度分布方程為即求出溫度分布方程為: :工程上遇

16、到的導(dǎo)熱問題，往往由于物體的幾何形狀復(fù)雜或邊界條件難以描工程上遇到的導(dǎo)熱問題，往往由于物體的幾何形狀復(fù)雜或邊界條件難以描述，無法求出分析解，此時可采用數(shù)值方法進行求解。常微分方程（述，無法求出分析解，此時可采用數(shù)值方法進行求解。常微分方程（ODE）包括初值問題和邊值問題兩種，初值問題包括初值問題和邊值問題兩種，初值問題ODE的數(shù)值解法常調(diào)用的數(shù)值解法常調(diào)用ode45（）（）和和ode23（）函數(shù)實現(xiàn)。（）函數(shù)實現(xiàn)。在在MATLAB編輯器中編寫函數(shù)編輯器中編寫函數(shù)BZ.m，存盤。，存盤。function BZclear allclca=0.0075;b=0.02;%r值的范圍值的范圍solinit=bvpini（tlinspace（a,b,10）,70 80）；）；%用初始值對解初始化用初始值