《【2017年整理】matlab_dsolve》由会员分享，可在线阅读，更多相关《【2017年整理】matlab_dsolve（14页珍藏版）》请在金锄头文库上搜索。
1、用 matlab 求解常微分方程在 MATLAB 中，由函数 dsolve()解决常微分方程（组）的求解问题，其具体格式如下：r = dsolve(eq1,eq2,., cond1,cond2,., v)eq1,eq2,.为微分方程或微分方程组， cond1,cond2,.,是初始条件或边界条件，v 是独立变量，默认的独立变量是t。函数 dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。例 1：求解常微分方程1dyx的 MATLAB 程序为：dsolve(Dy=1/(x+y),x) ，注意，系统缺省的自变量为 t，因此这里要把自变量写明。其中
2、：Y=lambertw(X)表示函数关系 Y*exp(Y)=X。例 2：求解常微分方程 20y的 MATLAB 程序为：Y2=dsolve(y*D2y-Dy2=0,x) Y2=dsolve(D2y*y-Dy2=0,x) 我们看到有两个解，其中一个是常数 0。例 3：求常微分方程组253ttdxyet通解的 MATLAB 程序为：X,Y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t),t) 例 4：求常微分方程组0221cos,4,tttdxyxttey通解的 MATLAB 程序为：X,Y=dsolve(Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y
3、=4*exp(-2*t),x(0)=2,y(0)=0,t) 以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。但是，我们知道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为 solver，其一般格式为：T,Y=solver(odefun,tspan,y0) 该函数表示在区间 tspan=t0,tf上，用初始条件 y0 求解显式常微分方程 (,)yft。solver 为命令 ode45， ode23，ode113，ode15s，ode23s，ode23t，o
4、de23tb 之一，这些命令各有特点。我们列表说明如下：求解器 特点 说明ode45一步算法，4,5 阶 Runge-Kutta方法累积截断误差 3()x大部分场合的首选算法ode23一步算法，2,3 阶 Runge-Kutta方法累积截断误差 3()x使用于精度较低的情形ode113多步法，Adams 算法，高低精度均可达到 3610计算时间比 ode45短ode23t 采用梯形算法 适度刚性情形ode15s 多步法，Gears 反向数值积分，精度中等 若 ode45 失效时，可尝试使用ode23s一步法，2 阶 Rosebrock算法，低精度。当精度较低时，计算时间比 ode15s短ode
5、fun 为显式常微分方程 (,)yft中的 (,)ftytspan 为求解区间，要获得问题在其他指定点 012,t 上的解，则令012,ftspantt（要求 it单调递增或递减）， y0 初始条件。例 5：求解常微分方程 2yx， .5x， ()y的 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：求解常
6、微分方程22(1)0,()1,(0)dydyytt的解，并画出解的图形。分析：这是一个二阶非线性方程(函数以及所有偏导数军委一次幂的是现性方程，高于一次的为非线性方程)，用现成的方法均不能求解，但我们可以通过下面的变换，将二阶方程化为一阶方程组，即可求解。令： 1xy， 2dt， 7，则得到：21122,(0)7(),dxtx解：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
7、;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）问题的形成轨道上的卫星，在牛顿第二定律 ，和万有引力定律 作用下有2drFmat 3EmMFGr，引力常数 G=6.672*10-11(N.m2/kg2) ,ME=5.97*1024(kg)是地球的质量。假23EMdraGrt定卫星以初速度 vy(0)
8、=4000m/s 在 x(0)=-4.2*107(m)处进入轨道。（2）构成一阶微分方程组令 Y=y1 y2 y3 y4T=x y vx vyT=x y x yT3142 123/234 /() ()xyExyvyYtGMaxA（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*6
9、0*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 求常微分方程的初值问题 38dyx， 02xy的解。r=dsolve(Dy+3*y=0,y(0)=2,x) 2利用 MAT
10、LAB 求常微分方程的初值问题 2(1)xy， 01xy， 03x的解。r=dsolve(D2y*(1+x2)-2*x*Dy=0,y(0)=1,Dy(0)=3,x) 3利用 MATLAB 求常微分方程 (4)20yy的解。解：y=dsolve(D4y-2*D3y+D2y,x) 4利用 MATLAB 求常微分方程组0324,230,tttdxyextty的特解。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求解常微分方程 2(1)0yy， 30x， ()1y， (0)的特解，并作出解函数的曲线图。r=dsolve(D
11、2y-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) 求解特殊函数方程勒让德方程的求解 22(1)(1)0()()dyxldlyx即 r=dsolve(1-x2)*D2y-2*x*Dy+y*l*(l+1)=0,x) 连带勒让德方程的求解：2 22(1) (1) 0dydymxxl yx r=dsolve(1-x2)*D2y-2*x*Dy+y*(l*(l+1
12、)-m2/(1-x2)=0,x) 贝塞尔方程 2 222()0dydyxxxvyr=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】求递推方程 的通解。fnffn()()()312（1）gs1=maple(rsolve(f(n)=-3*f(n-1)-2
13、*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 矩阵。fxyz（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 阶小量的台劳近似式。sin()xy2xy0,（1）TL1=maple(mtaylor(sin(x2+y2),x=0,y=0,8) （2）maple(readlib(mtaylo
14、r););TL2=maple(mtaylor(sin(x2+y2),x=0,y=0,8);pretty(sym(TL2) 运行 MAPLE 程序【例 5.7.3.2-1】目标：设计求取一般隐函数 的导数 解析解的程序，并要求：0),(yxf )(xy该程序能象 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(DYDZZY(x2*y(x)-exp(2*x)=sin(y(x)s3=maple(DYDZZY,cos(x+sin(y(x)=sin(y(x) （4）clear maplemexprocread(DYDZZY.src);maple(save(DYDZZY.m); （5）maple(read,
