matlab求解常微分方程

1、用matlab求解常微分方程在MATLAB中，由函數(shù)dsolve()解決常微分方程(組)的求解問題，其具體格式如下：r = dsolve('eq1,eq2,.', 'cond1,cond2,.', V)'eq1,eq2,為微分方程或微分方程組，con d1,c on d2,.'是初始條件或邊界條件，V是 獨立變量，默認的獨立變量是t'。函數(shù)dsolve用來解符號常微分方程、方程組，如果沒有初始條件，則求出通解，如 果有初始條件，則求出特解。dy 1例1 ：求解常微分方程dxx y的MATLAB程序為：dsolve('Dy=1/(x

2、+y)','x')，注意，系統(tǒng)缺省的自變量為t,因此這里要把自變量寫明其中：Y=lambertw(X)表示函數(shù)關(guān)系 Y*exp(Y)=X。例2：求解常微分方程yy''-y'2=0的MATLAB程序為：Y 2=dsolve('y*D2y-DyA2=0','x')Y 2=dsolve('D2y*y-DyA2=0','x')我們看到有兩個解，其中一個是常數(shù) 0dx 5x y = d例3:求常微分方程組dt通解的MATLAB程序為:X,Y =dsolve('Dx+5*x+y二exp(

3、t),Dy-x-3*y=exp(2*t)','t')fd+2d =10cost, xt =2 彳 dt dt7空+吐+2y = 4蘭y =0例4:求常微分方程組Idt dt y '兒弓 通解的MATLAB程序為:X,Y =dsolve('Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t)','x(0)=2,y(0)=0','t')以上這些都是常微分方程的精確解法，也稱為常微分方程的符號解。但是，我們知道，有大量的常微分方程雖然從理論上講，其解是存在的，但我們卻無法求出其解析解，此時，

4、我們需要尋求方程的數(shù)值解，在求常微分方程數(shù)值解方面，MATLAB具有豐富的函數(shù)，我們將其統(tǒng)稱為solver，其一般格式為：T,Y二solver(odefu n,tspa n,yO)該函數(shù)表示在區(qū)間tspan=tO,tf上，用初始條件y0求解顯式常微分方程yJf(t,y)。solver 為命令 ode45, ode23, ode113, ode15s, ode23s ode23t, ode23tb之一，這些命令各有特點。我們列表說明如下:求解 器特點說明ode45步算法，4,5階Runge-Kutta方法累積截斷誤差(")3大部分場合的首選 算法ode23步算法，2,3階Runge-K

5、utta方法累積截斷誤差Qx)3使用于精度較低的 情形ode113多步法，Adams算法， 高低精度均可達到1010工計算時間比ode45短ode23t采用梯形算法適度剛性情形ode15s多步法，Gear'反向 數(shù)值積分，精度中等若ode45失效時，可嘗試使用ode23s步法，2 階 Rosebrock 算法，低精度。當精度較低時，計算時間比ode15s 短odefun為顯式常微分方程y-f(t,y)中的f(t,y)tspan為求解區(qū)間，要獲得問題在其他指定點上。,上的解，則令tspan=t0,t1,t2右(要求1單調(diào)遞增或遞減)，y0初始條件。例5:求解常微分方程y-2y 2x2 2

6、x，0 25, y(0)的MATLAB程序如下： y=dsolve('Dy=-2*y+2*xA2+2*x','y(0)=1','x') x=0:0.01:0.5;yy二subs(y,x);x,y=ode15s(fun,0:0.01:0.5,1);ys二x.*x+exp(-fun二in li ne('- 2*y+2*x*x+2*x')2*x);plot(x,y,'r',x,ys,'b')例6:求解常微分方程栄的解，并畫出解的圖形。分析：這是一個二階非線性方程（函數(shù)以及所有偏導數(shù)軍委一次幕的是現(xiàn)性方程，

7、高 于一次的為非線性方程），用現(xiàn)成的方法均不能求解，但我們可以通過下面的變換，將二 階方程化為一階方程組，即可求解。=dy.令：x1=y，X懇，7，則得到:解:電dtdx2.dt2= 7(1 -為)X2為(0) =1X2(0) =0function dfy二mytt(t,fy) %f1二y;f2二dy/dt%求二階非線性微分方程時，把一階、二階直到(n-1)階導數(shù)用另外一個函數(shù)代替%用ode45命令時，必須表示成 Y'=f(t，丫)的形式%Y 二y1;y2;y3, 丫匕y1'y2'y3'=y2;y3;f(y1,y2,y3),%其中 y仁y,y2=y',y

8、3=y”%更高階時類似dfy=fy (2);7*(1-fy(1)A2)*fy (2)-fy(1);clear;clct,yy=ode45('mytt',0 40,1;0);plot(t,yy)legend('y','dy')【例4.1421-1】采用ODE解算指令研究圍繞地球旋轉(zhuǎn)的衛(wèi)星軌道。(1) 問題的形成軌道上的衛(wèi)星，在牛頓第二定律F二ma二m乓，和萬有引力定律F二-G呼：作用下有 dtr32a等一G 警，引力常數(shù) G=6.672*10-11(N.m2/kg2) ,M E=5.97*1024(kg)是地球的質(zhì)量。假 dtr定衛(wèi)星以初速度Vy(

9、0)=4000m/s在x(0)=-4.2*107(m)處進入軌道。(2) 構(gòu)成一階微分方程組令 Y=yi y2y3y4T=x y Vx VyT=x y x' y'Ty3Y'(t)=y 'i _vx'y'2Xyy'3ax-GM EY4yiy2,22、3/ 2(x y )(3) 根據(jù)上式dYdt mfun ction Y d=D Ydt(t，丫)% t% Yglobal G ME%xy=Y(1:2);Vxy 二Y (3:4);%r=sqrt(sum(xy.A2);Y d=Vxy;-G*ME*xy/rA3;%(4)global G ME%<

10、;1>G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;tspa n=tO,tf;%YO 二x0;0;0;vy0;%<8>t,YY=ode45('D Ydt',tspa n,Y 0);%X二YY (:,1);%丫二丫丫 （：,2）；%plot(X, Y,'b','Li newidth',2); hold on%axis('image')%XE, YE,ZE = sphere(10);%RE=0.64e7;%XE=RE*XE; YE二RE* Y

11、E;ZE=0*ZE;%mesh(XE,YE,ZE),hold off%練習：魚 +3v =81利用MATLAB求常微分方程的初值問題dx 3, 5=2的解r二dsolve('Dy+3*y=0','y(0)=2','x')2.利用MATLAB求常微分方程的初值問題(1 + x2)y'' = 2xy',人弓“ ,y'xd=3的解。r=dsolve('D2y*(1+xA2)-2*x*Dy=0','y(0)=1,Dy(0)=3','x')3.利用MATLAB求常微分方程y-2

12、y'y'' = o的解解：y=dsolve('D4y-2*D3y+D2y','x')4.利用 MATLAB求常微分方程組c dx , dyt2 4xy = e,dtdt生 3x y =0,dt3xy =?y 0的特解X, Y=dsolve('Dx*2+4*x+Dy-y=exp(t),Dx+3*x+y=0','x(0)=1.5,y(0)=0','t')25.求解常微分方程y''-2（1-y ）yy = o, 0沁豈30 , y（o）可,y'（o）=o的特解，并作出解函數(shù)

13、的曲線圖。r=dsolve('D2y-2*(1-yA2)*Dy+y=0','y(0)=0,Dy(0)=0','x')fun ction DY二 mytt2(t, Y)DY二Y (2);2*(1- Y(1)八2)* Y(2)+Y(1);clear;clct,YY=ode45('mytt2',0 30,1;0);plot(t,YY)legend('y','dy')求解特殊函數(shù)方程 勒讓德方程的求解ddX (1l(l 1)y =0(1-x2)1(11)y =0r=dsolve('(1-xA2)*D

14、2y-2*x*Dy+y*l*(l+1)=0','x')連帶勒讓德方程的求解:(1 X2)dXy - -2xdx-l(l 1)-1*廠0r=dsolve('(1-xA2)*D2y-2*x*Dy+y*(l*(l+1)-mA2/(1-xA2)=0','x')貝塞爾方程dx2xd (x2 V2)廠 0dxr=dsolve('xA2*D2y+x*Dy+(xA2-vA2)*y=0','x')利用maplemaple dsolve(1-xA2)*diff(y(x),x$2)-2*x*diff(y(x),x)+n*(n+1

15、)*y(x)=0, y(x)利用MAPLE的深層符號計算資源經(jīng)典特殊函數(shù)的調(diào)用MAPLE庫函數(shù)在線幫助的檢索樹發(fā)揮MAPLE的計算潛力 調(diào)用MAPLE函數(shù)【例5.731-1】求遞推方程f (n) 3f (n 1) 2f( n 一2)的通解(1)gs 仁m aple('rsolve(f( n)=-3*f( n-1)-2*f( n-2),f(k);')(2)調(diào)用格式二gs2=maple('rsolve','f( n)=-3*f( n-1)-2*f( n-2)','f(k)')【例5.731-2】求f =xyz的Hessian矩陣(1)

16、FH1=maple('hessia n(x*y*z,x,y,z);')FH2=maple('hessia n','x*y*z','x,y,z')FH二sym(FH2)【例5.731-3】求sin(x2 y2)在x=O,y=O處展開的截斷8階小量的臺勞近似式(1)TL1二maple('mtaylor(si n( xA2+yA2),x=0,y=0,8)')(2)maple('readlib(mtaylor);');TL2=maple('mtaylor(si n( xA2+yA2),x=0,y=0

17、,8)');pretty(sym(TL2)運行MAPLE程序【例5.732-1】目標：設計求取一般隱函數(shù)f(x,y)=0的導數(shù)y(x)解析解的程序，并要 求：該程序能象MAPLE原有函數(shù)一樣可被永久調(diào)用。1)DYDZZY.srcDYDZZY:=proc(f)# DYDZZY(f) is used to get the derivate of#an implicit functionlocal Eq,deq,imderiv;Eq:=f;deq:=diff(Eq,x);readlib(isolate);imderiv:=isolate(deq,diff(y(x),x);end;(2) procread('DYDZZY.src')(3) s1=maple('DYDZZY(x=log(x+y(x);')s2二maple('D YDZZ Y( x2*y(x)-exp(2*x)二si n(