用matlab求解常微分方程在MATLAB中，由函数dsolve()解决常微分方程(组)的求解问题，其具体格式如下：r 二dsolve('eql,eq2,・••字condl,cond2,・.・；V)匕ql,eq2,・・・*为微分方程或微分方程组，，condl,cond2,.・・；是初始条件或边界条件，P是独立变量，默认的独立变量是讥函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。

dy _1例1：求解常微分方程莎一石的MATLAB程序为：dsolve(* Dy=l/(x+y) 1r!x1),注意，系统缺省的自变量为t,因此这里要把自变量写明。

其中：Y=lambertw(X)表示函数关系Y*exp(Y)二X。

例2：求解常微分方程E'-y— 0的MATLAB程序为：Y2=dsolve(1y*D2y-Dy A2=01, 1x f)Y2=dsolve(!D2y*y-Dy A2=0 J )我们看到有两个解，其中一个是常数0。

dx 心 ？ —+ 5x + y = e dt空_兀_3『= g2f 例3：求常微分方程组〔力 ' 通解的MATLAB 程序为： [X,Y]=dsolve(f Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t) 1, 111) [X,Y]=dsolve(1Dx+2 *x-Dy=l0 * cos(t),Dx+Dy+2 *y=4 *exp(-2*t) T ,f x(0)=2,y(0)=0f ,f t T) 以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。

但是，我们知 道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析 解，此吋，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰 富的函数，我们将其统称为solver,其一般格式为：i°cosr, 7=2 例4：求常微分方程组 y = 0 z 通解的MATLAB 程序为:[T,Y]=solver(odefun,tspan,yO)该函数表示在区间tspan=[tO,tf]±,用初始条件yO求解显式常微分方程卩="，刃。

solver 为命令odc45, odc23, ode 113, ode 15s, ode23s, ode23t, ode23tb 之一,这些命令各有特点。

我们列表说明如下:odefun为显式常微分方程）,”，刃中的“丿）tspan为求解区间，要获得问题在其他指定点心，也，…上的解，则令切期二［心,/“，…心］（要求-单调递增或递减），yO初始条件。

例5：求解常微分方程y' = -2y + 2x2+2x, 0<x<0.5, y（0） = l的MATLAB程序如下：y=dsolve（1Dy=-2*y+2*x A2+2*x1, 1y （0）=11, 1x!） x=0:0.01:0.5;yy=subs(y,x);fun=inline (*-2*y+2*x*x+2*x1); [x,y]=odel5s(fun, [0:0.01:0.5],1);ys=x ・ *x+exp(一2*x);Plot(x,y, f r f,x,ys, f b f )例6：求解常微分方程分冷+尸°"3)"'八°2()的解，并画出解的图形。

分析：这是一个二阶非线性方程(函数以及所有偏导数军委一次幕的是现性方程，高于一次的为非线性方程)，用现成的方法均不能求解，但我们可以通过下面的变换，将二阶方程化为一阶方程组，即可求解。

= dy_令：兀产儿氐-亦，“ =7,则得到：~~ =兀2 ，兀| (0) = 1at<= 7(1-%,2 )x2 - Xj, x2 (0) = 0dt解:function [dfy]=mytt(t,fy)%fl=y;f2=dy/dt%求二阶非线性微分方程时，把一阶、二阶直到(ml)阶导数用另外一个函数代替%用ode45命令时，必须表示成Y』f(t,Y)的形式%Y=[yl;y2;y3],Y,=[yl\y2\y3>[y2;y3;f(yl,y2,y3)],%其中yl=y,y2=y',y3=y"%更高阶时类似dfy=[fy(2);7*(l ・fy ⑴八2)*fy(2)・fy ⑴]；clear;clc[t,yy]=ode45Cmytt\[0 40],[1;0]);plot(t,yy)legend('yTdy')【例4.14.2.1-1］采用ODE解算扌旨令研究围绕地球旋转的卫星轨道。

(1)问题的形成轨道上的卫星，在牛顿第二定律F=ma = m^,和万有引力定律F T叫；作用下有dt r2a^L = -G^r ,引力常数G=6.672*10 H(N.m2/kg2) ,ME=5.97*1024(kg)是地球的质量。

假dr r定卫星以初速度v y(0)=4000m/s在x(0)=4.2*l()7(m)处进入轨道。

(2) 构成一阶微分方程组令 Y 二[刃力力为]「二[尤 V V x v y ]T =[x y * y]T儿 (…严 儿 (宀)计(3) 根据上式[dYdt.m]function Ycl=DYdt (t, Y)% t% Yglobal G ME% xy 二Y(l:2) ;Vxy=Y(3:4) ; %r=sqrt(sum(xy.八2)); Yd 二[Vxy ;-G*ME*xy//3] ; %(4)global G ME % <1>G=6 ・ 672e-ll;ME=5 ・ 97e24;vy0=4000;x0=-4 ・ 2e7;t0=0;tf=60*60*24*9; tspan=[tO,tf];%-GM E —GMY0=[xO；0;0;vyO]; %X=YY(:,1);Y=YY(:,2) ; %plot(X,Y, , linewidth1,2); hold on%axis(* image 1) %[XE,YE,ZE] = sphere (10); %RE=0・64e7; %XE=RE*XE;YE=RE*YE;ZE=0*ZE; %mesh(XE,YE,ZE),hold off %练习：dy_1.利用MATLAB求常微分方程的初值问题杰* " , >?L=o = 2的解。

r=dsolve(1Dy+3*y=01, !y (0)=21, 1x *)2.利用MATLAB求常微分方程的初值问题（1 +扌）* = 23,儿=。

二1,儿乜”的解。

r=dsolve(1D2y*(l+x A2)-2*x*Dy=0',1y(0)=1,Dy(0)=31, *)3. 利用MATLAB 求常微分方程y ⑷-2厂+* = 0的解。

解:y=dsolve(,D4y-2*D3y+D2y ,/x ,) [X,Y]=dsolve(1Dx*2+4*x+Dy-y=exp(t),Dx+3*x+y=01, 1x(0)=1.5,y (0)=01, 丫 t 】) 5.求解常微分方程/，-2(l-/)/+y = 0, 0<x<30, y(0) = l, /(0) = 0的特解，并作出解函 数的曲线图。

r=dsolve(1D2y-2*(l-y A 2)*Dy+y=01, 1y(0)=0,Dy(0)=01, * x 1)4.利用MATLAB 求常微分方程组r dx A dy.2— + 4x + -- y = e dt dt 3 X l /=0 ~ 2 y = 0 /=0的特解。

function DY=mytt2(t,Y)DY=[Y(2);2*( 1 -Y(l )A2)*Y(2)+Y( 1)]; clear;clc[t,YY]=ode45(,mytt2\[0 30],[l;0]); plot(t,YY)legendC^Vdy 1) 求解特殊函数方程勒让德方程的求解dy dx + /(/ + l)y = 0dx d , o x—(1一厂)即dy —2— + /(/ + l)y = 0 dxr=dsolve(1(l-x A2)*D2y-2*x*Dy+y*l*(1+1)=01,1x')连带勒让徳方程的求解:r=dsolve(1(l-x A 2)*D2y-2*x*Dy+y*(1*(l+l)-m A 2/(l-x A 2))=0\)贝塞尔方程r=dsolve(!x A 2*D2y+x*Dy+(x A 2-v A 2)*y=01)利用maplemaple dsolve((l-x A 2)*diff(y(x),x$2)-2*x*diff(y(x)r x)+n*(n+1)*y(x) =0, y (x)) 利用MAPLE 的深层符号计算资源d)等—力孚+ (+1)- 2m i 21-x y = 0经典特殊函数的调用MAPLE库函数在线帮助的检索树发挥MAPLE的计算潜力调用MAPLE函数【例 5.7.3.1-1 ]求递推方程f(n) = -3/(n- 1) -2/(n-2)的通解。

(1)gsl=maple(!rsolve(f(n)=-3*f(n-1)-2*f (n-2),f (k)); 1)(2)调用格式二gs2=maple(T rsolve1, T f(n)=-3*f(n-1)-2*f(n-2) T, !f(k) !)【例5.7.3.1-2]求/* =供的Hessian矩阵。

(1)FHl=maple(!hessian(x*y*z f [x,y,z]); *)FH2=maple(* hessian!, * x*y*z1, 1 [x,y,z] 1)FH=sym(FH2)【例5.7.3」・3】求sin(>2 + y2)在兀二Q y = 0处展开的截断8阶小量的台劳近似式。

(1)TLl=maple(^taylor(sin(x A2+y A2),[x=0,y=0],8)1)(2)maple(1 readlib(mtaylor); 1);TL2=maple(1mtaylor(sin(x A2+y A2), [x=0,y=0],8) 1);pretty(sym(TL2))运行MAPLE程序【例5.7.321】目标：设计求取一般隐函数f(x.y) = 0的导数才⑴解析解的程序，并要求：该程序能象MAPLE原有函数一样可被永久调用。

[DYDZZY.src]DYDZZY:=proc (f)#DYDZZY (f) is used to get the derivate of#an implicit functionlocal Eq,deq,imderiv;Eq: = 1 Eq 1 ;Eq:=f;deq:=diff(Eq, x);readlib(isolate);imderiv:=isolate(deq,diff(y(x) , x));end;procread(1DYDZZY•src1)(3)sl=maple(1DYDZZY(x=log(x+y(x)));1)s2=maple(f DYDZZY(x A2*y(x)-exp(2*x)=sin(y(x))) !) s3=maple(1DYDZZY1, 1 cos(x+sin(y(x)))=sin(y(x)) 1)clear maplemexprocread(1DYDZZY ・ src1);maple(1 save(、DYDZZY.m s) 1);(5)maple(* read'J'DYDZZY ・m、T);ss2=maple(1DYDZZY(x A2*y(x)-exp(2*x)=sin(y(x)))1)。