matlab_dsolve.doc

用matlab求解常微分方程在MATLAB中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下：r = dsolve(eq1,eq2,., cond1,cond2,., v)eq1,eq2,.为微分方程或微分方程组，cond1,cond2,.,是初始条件或边界条件，v是独立变量，默认的独立变量是t。函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。例1：求解常微分方程的MATLAB程序为：dsolve(Dy=1/(x+y),x) ，注意，系统缺省的自变量为t，因此这里要把自变量写明。其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X。例2：求解常微分方程的MATLAB程序为：Y2=dsolve(y*D2y-Dy2=0,x) Y2=dsolve(D2y*y-Dy2=0,x) 我们看到有两个解，其中一个是常数0。例3：求常微分方程组通解的MATLAB程序为：X,Y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t),t) 例4：求常微分方程组通解的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) 以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。但是，我们知道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB具有丰富的函数，我们将其统称为solver，其一般格式为：T,Y=solver(odefun,tspan,y0) 该函数表示在区间tspan=t0,tf上，用初始条件y0求解显式常微分方程。solver为命令ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb之一，这些命令各有特点。我们列表说明如下：求解器特点说明ode45一步算法，4,5阶Runge-Kutta方法累积截断误差大部分场合的首选算法ode23一步算法，2,3阶Runge-Kutta方法累积截断误差使用于精度较低的情形ode113多步法，Adams算法，高低精度均可达到计算时间比ode45短ode23t采用梯形算法适度刚性情形ode15s多步法，Gears反向数值积分，精度中等若ode45失效时，可尝试使用ode23s一步法，2阶Rosebrock算法，低精度。当精度较低时，计算时间比ode15s短odefun为显式常微分方程中的tspan为求解区间，要获得问题在其他指定点上的解，则令（要求单调递增或递减），y0初始条件。例5：求解常微分方程，的MATLAB程序如下：y=dsolve(Dy=-2*y+2*x2+2*x,y(0)=1,x)x=0:0.01:0.5;yy=subs(y,x); fun=inline(-2*y+2*x*x+2*x);x,y=ode15s(fun,0:0.01:0.5,1);ys=x.*x+exp(-2*x);plot(x,y,r,x,ys,b) 例6：求解常微分方程的解，并画出解的图形。分析：这是一个二阶非线性方程(函数以及所有偏导数军委一次幂的是现性方程，高于一次的为非线性方程)，用现成的方法均不能求解，但我们可以通过下面的变换，将二阶方程化为一阶方程组，即可求解。令：，则得到：解：function dfy=mytt(t,fy)%f1=y;f2=dy/dt%求二阶非线性微分方程时，把一阶、二阶直到(n-1)阶导数用另外一个函数代替%用ode45命令时，必须表示成Y=f(t,Y)的形式%Y=y1;y2;y3,Y=y1;y2;y3=y2;y3;f(y1,y2,y3),%其中y1=y,y2=y,y3=y%更高阶时类似dfy=fy(2);7*(1-fy(1)2)*fy(2)-fy(1);clear;clct,yy=ode45(mytt,0 40,1;0);plot(t,yy)legend(y,dy) 【例4.14.2.1-1】采用ODE解算指令研究围绕地球旋转的卫星轨道。（1）问题的形成轨道上的卫星，在牛顿第二定律，和万有引力定律作用下有 ，引力常数G=6.672*10-11(N.m2/kg2) ,ME=5.97*1024(kg)是地球的质量。假定卫星以初速度vy(0)=4000m/s在x(0)=-4.2*107(m)处进入轨道。（2）构成一阶微分方程组令Y=y1 y2 y3 y4T=x y vx vyT=x y x yT（3）根据上式dYdt.mfunction Yd=DYdt(t,Y) % t% Yglobal G ME%xy=Y(1:2);Vxy=Y(3:4);%r=sqrt(sum(xy.2);Yd=Vxy;-G*ME*xy/r3;%（4）global G ME%G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;tspan=t0,tf;%Y0=x0;0;0;vy0;%t,YY=ode45(DYdt,tspan,Y0);%X=YY(:,1);%Y=YY(:,2);%plot(X,Y,b,Linewidth,2); hold on%axis(image)%XE,YE,ZE = sphere(10);%RE=0.64e7;%XE=RE*XE;YE=RE*YE;ZE=0*ZE;%mesh(XE,YE,ZE),hold off% 练习：1利用MATLAB求常微分方程的初值问题，的解。r=dsolve(Dy+3*y=0,y(0)=2,x) 2利用MATLAB求常微分方程的初值问题，的解。r=dsolve(D2y*(1+x2)-2*x*Dy=0,y(0)=1,Dy(0)=3,x) 3利用MATLAB求常微分方程的解。解：y=dsolve(D4y-2*D3y+D2y,x) 4利用MATLAB求常微分方程组的特解。X,Y=dsolve(Dx*2+4*x+Dy-y=exp(t),Dx+3*x+y=0,x(0)=1.5,y(0)=0,t) 5求解常微分方程，的特解，并作出解函数的曲线图。r=dsolve(D2y-2*(1-y2)*Dy+y=0,y(0)=0,Dy(0)=0,x) function 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) 求解特殊函数方程勒让德方程的求解r=dsolve(1-x2)*D2y-2*x*Dy+y*l*(l+1)=0,x) 连带勒让德方程的求解：r=dsolve(1-x2)*D2y-2*x*Dy+y*(l*(l+1)-m2/(1-x2)=0,x) 贝塞尔方程r=dsolve(x2*D2y+x*Dy+(x2-v2)*y=0,x) 利用maplemaple dsolve(1-x2)*diff(y(x),x$2)-2*x*diff(y(x),x)+n*(n+1)*y(x) = 0, y(x) 利用MAPLE的深层符号计算资源经典特殊函数的调用MAPLE库函数在线帮助的检索树发挥MAPLE的计算潜力调用MAPLE函数【例5.7.3.1-1】求递推方程的通解。（1）gs1=maple(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.7.3.1-2】求的Hessian矩阵。（1）FH1=maple(hessian(x*y*z,x,y,z);) （2）FH2=maple(hessian,x*y*z,x,y,z) （3）FH=sym(FH2) 【例5.7.3.1-3】求在处展开的截断8阶小量的台劳近似式。（1）TL1=maple(mtaylor(sin(x2+y2),x=0,y=0,8) （2）maple(readlib(mtaylor););TL2=maple(mtaylor(sin(x2+y2),x=0,y=0,8);pretty(sym(TL2) 运行MAPLE程序【例5.7.3.2-1】目标：设计求取一般隐函数的导数解析解的程序，并要求：该程序能象MAPLE原有函数一样可被永久调用。（1）DYDZZY.srcDYDZZY:=proc(f)# DYDZZY(f) is used to get the derivate of# an implicit functionlocal Eq,deq,imderiv;Eq:=Eq;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