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

1、MATLAB解常微分方程及其應(yīng)用授課教師：授課教師：唐躍龍?zhí)栖S龍TelephoneQ: 794218136Email:1 微分方程的基本概念微分方程的基本概念定義1 含有未知函數(shù)的導(dǎo)數(shù)(或微分)，同時(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ù)(或微分)的項(xiàng)稱為微分方程自由項(xiàng).定義6 自由項(xiàng)為零的微分方程稱為齊次方程.定義7 自由項(xiàng)不為零的微分方程稱為非齊次方程.例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 若微分方程的解中所含相互獨(dú)立的任意常數(shù)的個(gè)數(shù)與微分方程的階

3、數(shù)相同，這樣的解叫做微分方程的通解.例2. 驗(yàn)證函數(shù) 是微分方程 的解.255xye210yy 例3. 求方程 的通解.6yx 例4. 求方程 滿足初始條件：當(dāng) 時(shí) 的特解.23yx 1x 2y 注：這里所求的微分方程的解、通解、特解都是解析解（精確解）.2 MATLAB求微分方程的解析解求微分方程的解析解格式一：y=dsolve(f1, f2, , fm)格式二：y=dsolve(f1, f2, , fm ,x) %指明自變量 fi既可以描述微分方程，又可描述初始條件或邊界條件。如： 描述微分方程時(shí) 描述條件時(shí) . 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. 因此，不是所有方程都有解析解，只有極少部分的微分方程有解析解.當(dāng)找不到原方程的解析解時(shí)，我們只能考慮其數(shù)值解.1)(1 (*)()(2txtxtx以下面的一階常微分方程以下面的一階常微分方程(

5、ODE )初值問(wèn)題為例初值問(wèn)題為例 ：000( , ) ()ndyf x yxxxdxy xy()nnyy x 數(shù)值解法就是求數(shù)值解法就是求y(x)在某些分立的節(jié)點(diǎn)在某些分立的節(jié)點(diǎn) xn 上的近似值上的近似值yn，用以近似用以近似y(xn)3.1 簡(jiǎn)單歐拉簡(jiǎn)單歐拉(L.Euler, 1707-1783)方法。方法。00( , )()dyf y xdxy xy012nyyyy 歐拉數(shù)值算法就是由初值通過(guò)歐拉數(shù)值算法就是由初值通過(guò)遞推求解遞推求解，遞推求解，遞推求解就是從初值開(kāi)始，后一個(gè)函數(shù)值由前一個(gè)函數(shù)值得到。關(guān)就是從初值開(kāi)始，后一個(gè)函數(shù)值由前一個(gè)函數(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對(duì)于等間隔節(jié)點(diǎn)對(duì)于等間隔節(jié)點(diǎn)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é)點(diǎn)上，微分方程可以寫為節(jié)點(diǎn)上，微分方程

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

8、量的初值和終值表示自變量的初值和終值%y0表示函數(shù)在表示函數(shù)在x0處的值處的值,其可以為向量形式其可以為向量形式%PointNum表示自變量在表示自變量在x0,xt上取的點(diǎn)數(shù)上取的點(diǎn)數(shù)if nargin5 | PointNum=0 %如果函數(shù)僅輸入如果函數(shù)僅輸入4個(gè)參數(shù)值，則個(gè)參數(shù)值，則PointNum默認(rèn)值為默認(rèn)值為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的符號(hào)解法的符號(hào)解法,求解常微分方程的特解：求解常微分方程的特解：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求解描述某剛性問(wèn)題的微分方程：求解描述某剛性問(wèn)題的微分方程：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)熱計(jì)算傳熱過(guò)程涉及面廣，數(shù)學(xué)模型復(fù)雜。計(jì)算過(guò)程中涉及到傳熱過(guò)程涉及面廣，數(shù)學(xué)模型復(fù)雜。計(jì)算過(guò)程中涉及到許多運(yùn)算方法。導(dǎo)熱雖是容易做數(shù)學(xué)處理的一種熱量傳遞許多運(yùn)算方法。導(dǎo)熱雖是容易做數(shù)學(xué)處理的一種熱量傳遞方式，但其過(guò)程往往涉及常微分、偏微分方程、線性（非方式，但其過(guò)程往往涉及常微分、偏微分方程、線性（非線性）方程組的求解。線性）方程組的求解。有一外徑為

14、有一外徑為4cm，內(nèi)徑為，內(nèi)徑為1.5cm，載有電流密度，載有電流密度I為為5000A/cm2的內(nèi)冷鋼質(zhì)導(dǎo)體。導(dǎo)體單位時(shí)間發(fā)出的熱量等于流體同時(shí)帶走的內(nèi)冷鋼質(zhì)導(dǎo)體。導(dǎo)體單位時(shí)間發(fā)出的熱量等于流體同時(shí)帶走的熱量，導(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。解：這是圓柱坐標(biāo)中常物性一維穩(wěn)態(tài)導(dǎo)熱問(wèn)題，結(jié)合本題具解：這是圓柱坐標(biāo)中常物性一維穩(wěn)態(tài)導(dǎo)熱問(wèn)題，結(jié)合本題具體條件導(dǎo)熱微分方程式可簡(jiǎn)化為：體條件導(dǎo)熱微分方程

15、式可簡(jiǎn)化為：結(jié)合邊界條件，可得這一導(dǎo)熱問(wèn)題的數(shù)學(xué)描述為：結(jié)合邊界條件，可得這一導(dǎo)熱問(wèn)題的數(shù)學(xué)描述為：此常微分方程的分析解，可調(diào)用此常微分方程的分析解，可調(diào)用MATLAB符號(hào)工具箱中的符號(hào)工具箱中的dsolve函數(shù)來(lái)實(shí)現(xiàn)。在命令窗口執(zhí)行下面的代碼：函數(shù)來(lái)實(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)熱問(wèn)題，往往由于物體的幾何形狀復(fù)雜或邊界條件難以描工程上遇到的導(dǎo)熱問(wèn)題，往往由于物體的幾何形狀復(fù)雜或邊界條件難以描述，無(wú)法求出分析解，此時(shí)可采用數(shù)值方法進(jìn)行求解。常微分方程（述，無(wú)法求出分析解，此時(shí)可采用數(shù)值方法進(jìn)行求解。常微分方程（ODE）包括初值問(wèn)題和邊值問(wèn)題兩種，初值問(wèn)題包括初值問(wèn)題和邊值問(wèn)題兩種，初值問(wèn)題ODE的數(shù)值解法常調(diào)用的數(shù)值解法常調(diào)用ode45（）（）和和ode23（）函數(shù)實(shí)現(xiàn)。（）函數(shù)實(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）；）；%用初始值對(duì)解初始化用初始值