matlab求解常微分方程

1、用matlab求解常微分方程在MATLAB中，由函数dsolve()解决常微分方程(组)的求解问题，其具体格式如下：r = dsolve('eq1,eq2,.', 'cond1,cond2,.', V)'eq1,eq2,为微分方程或微分方程组，con d1,c on d2,.'是初始条件或边界条件，V是 独立变量，默认的独立变量是t'。函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如 果有初始条件，则求出特解。dy 1例1 ：求解常微分方程dxx y的MATLAB程序为：dsolve('Dy=1/(x

2、+y)','x')，注意，系统缺省的自变量为t,因此这里要把自变量写明其中：Y=lambertw(X)表示函数关系 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')我们看到有两个解，其中一个是常数 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、我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB具有丰富的函数，我们将其统称为solver，其一般格式为：T,Y二solver(odefu n,tspa n,yO)该函数表示在区间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'反向 数值积分，精度中等若ode45失效时，可尝试使用ode23s步法，2 阶 Rosebrock 算法，低精度。当精度较低时，计算时间比ode15s 短odefun为显式常微分方程y-f(t,y)中的f(t,y)tspan为求解区间，要获得问题在其他指定点上。,上的解，则令tspan=t0,t1,t2右(要求1单调递增或递减)，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:求解常微分方程栄的解，并画出解的图形。分析：这是一个二阶非线性方程（函数以及所有偏导数军委一次幕的是现性方程，

7、高 于一次的为非线性方程），用现成的方法均不能求解，但我们可以通过下面的变换，将二 阶方程化为一阶方程组，即可求解。=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)阶导数用另外一个函数代替%用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解算指令研究围绕地球旋转的卫星轨道。(1) 问题的形成轨道上的卫星，在牛顿第二定律F二ma二m乓，和万有引力定律F二-G呼：作用下有 dtr32a等一G 警，引力常数 G=6.672*10-11(N.m2/kg2) ,M E=5.97*1024(kg)是地球的质量。假 dtr定卫星以初速度Vy(

9、0)=4000m/s在x(0)=-4.2*107(m)处进入轨道。(2) 构成一阶微分方程组令 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) 根据上式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的特解，并作出解函数

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')求解特殊函数方程 勒让德方程的求解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的深层符号计算资源经典特殊函数的调用MAPLE库函数在线帮助的检索树发挥MAPLE的计算潜力 调用MAPLE函数【例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)调用格式二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】目标：设计求取一般隐函数f(x,y)=0的导数y(x)解析解的程序，并要 求：该程序能象MAPLE原有函数一样可被永久调用。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(