资源预览内容
第1页 / 共12页
第2页 / 共12页
第3页 / 共12页
第4页 / 共12页
第5页 / 共12页
第6页 / 共12页
第7页 / 共12页
第8页 / 共12页
第9页 / 共12页
第10页 / 共12页
亲,该文档总共12页,到这儿已超出免费预览范围,如果喜欢就下载吧!
资源描述
MATLAB 解微分方程(2011-07-15 17:35:25)转 载标签: 教育分类: matlab 学习用 matlab 时间也不短了,可是一直没有接触过微分方程。这次看看书,学习学习,记点儿 笔记。 1.可以解析求解的微分方程。dsolve() 调用格式为:y=dsolve(f1,f2,.,fmO;y=dsolve(f1,f2,.,fm,x);如下面的例子,求解了微分方程syms t;u=exp(-5*t)*cos(2*t-1)+5;uu=5*diff(u,t,2)+4*diff(u,t)+2*u;syms t y;y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=87*exp(-5*t)*cos(2*t-1)+92*exp(-5*t)*sin(2*t-1)+10)yc=latex(y)将 yc 的内容 copy 到 latex 中编译,得到结果。关于 Matlab 的微分方程,直到今天才更新第 2 篇,实在是很惭愧的事因为原因都在于 太懒惰,而不是其他的什么。在上一篇中,我们使用 dsolve 可以解决一部分能够解析求解的微分方程、微分方程组,但是对于大多数微分方程(组)而言不能得到解析解,这时数值求解也就是没有办法的办法 了,好在数值解也有很多的用处。数值分析方法中讲解了一些 Eular 法、 Runge-Kutta 法等一些方法,在 matlab 中内置的 ode 求解器可以实现不同求解方法的相同格式的调用,而不必太关心 matlab 究竟是用什么 算法完成的。这一回我们来说明 ode45 求解器的使用方法。1.ode45 求解的上手例子:求解方程组Dx=y+x(1-x2-y2);Dy=-x+y*(1-x2-y2)初值 x=0.1;y=0.2;先说明一下最常用的 ode45 调用方式,和相应的函数文件定义格式。t,x=ode45(odefun,tspan,x0);其中,Fun 就是导函数,tspan 为求解的时间区间(或时间序列,如果采用时间序列,则必 须单调),x0 为初值。这时,函数文件可以采用如下方式定义function dx=odefun(t,x)对于上面的小例子,可以用如下的程序求解。function jixianhuanclear;clcx0=0.1;0.2;t,x=ode45(jxhdot,0,100,x0);plot(x(:,1),x(:,2)function dx=jxhdot(t,x)dx=x(2)+x(1).*(1-x(1).2-x(2).2);-x(1)+x(2).*(1-x(1).2-x(2).2);2.终值问题tspan 可以是递增序列,也可以为递减序列,若为递减则可求解终值问题。t,x=ode45(zhongzhiode,3,0,1;0;2);plot(t,x)function dx=zhongzhiode(t,x) dx=2*x(2)2-2; -x(1)+2*x(2)*x(3)-1;-2*x(2)+2*x(3)2-4;结果如下3.odesetoptions = odeset(name1,value1,name2,value2,.)t,x=solver(fun,tspan,x0,options)通过 odeset 设置 options第一,通过求解选项的设置可以改善求解精度,使得原本可能不收敛的问题收敛。options=odeset(RelTol,1e-10);第二,求解形如 M(t,x)x=f(t,x)的方程。例如,方程x=-0.2x+yz+0.3xyy=2xy-5yz-2y2x+y+z-2=0可以变形为1 0 0x -0.2x+yz+0.3xy0 1 0y=2xy-5yz-2y2 0 0 1z x+y+z-2 这样就可以用如下的代码求解该方程function mydae M=1 0 0;0 1 0;0 0 0; options=odeset(Mass,M); x0=1.6,0.3,0.1; t,x=ode15s(daedot,0,1.5,x0,options);plot(t,x)function dx=daedot(t,x) dx=-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);x(1)+x(2)+x(3)-2;4.带附加参数的 ode45有时我们需要研究微分方程组中的参数对于解的影响,这时采用带有参数的 ode45 求解会 使求解、配合循环使用,可以使得求解的过程更加简捷。使用方法:只需将附加参数放在 options 的后面就可以传递给 odefun 了。看下面的例子。function Rosslerclear;clc a=0.2,0.2; b=0.2,0.5; c=5.7,10;x0=0 0 0; for jj=1:2t,x=ode45(myRossler,0,100,x0,a(jj),b(jj),c(jj);figure;plot3(x(:,1),x(:,2),x(:,3);grid on; endfunction dx=myRossler(t,x,a,b,c)dx=-x(2)-x(3);x(1)+a*x(2);b+(x(1)-c)*x(3);5. 刚性方程的求解刚性方程就是指各个自变量的变化率差异很大,会造成通常的求解方法失效。这是 matlab 中自带的一个例子,使用 ode15s 求解,如果用 ode45 求解就会出现错误。function myode15studyt,Y = ode15s(vdp1000,0 3000,2 0); plot(T,Y(:,1),-o)figure;plot(Y(:,1),Y(:,2)function dy = vdp1000(t,y) dy = zeros(2,1); dy(1) = y(2); dy(2) = 1000*(1 - y(1)2)*y(2) - y(1);6.高阶微分方程的求解通常的方法是进行变量替换,将原方程降阶,转换成更多变量的一阶方程组进行求解。在这个例子里我们求解一个动力学系统里最常见的一个运动方程,其中 f=sin(t)function myhighoder clear;clc x0=zeros(6,1); t,x=ode45(myhigh,0,100,x0); plot(t,x(:,1)function dx=myhigh(t,x) f=sin(t);0;0; M=eye(3); C=eye(3)*0.1; K=eye(3)-0.5*diag(ones(2,1),1)-0.5*diag(ones(2,1),-1); dx=x(4:6);inv(M)*(f-C*x(4:6)-K*x(1:3);7.延迟微分方程matlab 提供了 dde23 求解非中性微分方程。dde23 的调用格式如下:sol = dde23(ddefun,lags,history,tspan)lags 是延迟量,比如方程中包含 y1(t-0.2)和 y2(t-0.3)则可以使用 lags=0.2,0.3。这里的 ddefun 必须采用如下的定义方式:dydt = ddefun(t,y,Z)其中的 Z(:,1)就是 y(t-lags(1),Z(:,2)就是 y(t-lags(2)).下面是个使用 dde23 求解延迟微分方程的例子。function mydde23study % The differential equations % y_1(t) = y_1(t-1) % y_2(t) = y_1(t-1)+y_2(t-0.2) % y_3(t) = y_2(t) % are solved on 0, 5 with history y_1(t) = 1, y_2(t) = 1, y_3(t) = 1 for % t = 0. clear;clc lags=1,0.2; history=1;1;1; tspan=0,5; sol = dde23(myddefun,lags,history,tspan)plot(sol.x,sol.y)function dy = myddefun(t,y,Z) dy=Z(1,1);Z(1)+Z(2,2);y(2) ;8.ode15i 求解隐式微分方程T,Y = ode15i(odefun,tspan,y0,yp0)yp0 为 y的初值。odefun 的格式如下 dy = odefun(t,y,yp),yp 表示 y,而方程中应该使得 f(t,y,y)=0function myodeIMP % The problem is % % y(1) = -0.04*y(1) + 1e4*y(2)*y(3)% y(2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)2 % y(3) = 3e7*y(2)2 % % It is to be solved with initial conditions y(1) = 1, y(2) = 0, y(3) = 0 % to steady state. clear;clc y0=1;0;0;fixed_y0=1;1;1; yp0=0 0 0; fixed_yp0=;y0mod,yp0mod=decic(myodefunimp,0,y0,fixed_y0,yp0,fixed_yp0); tspan=0, logspace(-6,6); t,y = ode15i(myodefunimp,tspan,y0mod,yp0mod); y(:,2)=1e4*y(:,2); semilogx(t,y) function res=myodefunimp(t,y,yp) res=-yp(1)-0.04*y(1)+1e4*y(2)*y(3);-yp(2)+0.04*y(1)-1e4*y(2)*y(3)-3e7*y(2)2;-yp(3)+3e7*y(2)2;这次要接触一个新的求解 ode 的方法,就是使用 simulink 的积分器求解。1.还是做我们研究过的一个例子(在初识 matlab 微分方程(2)中采用的)。Dx=y+x(1-x2-y2);Dy=-x+y*(1-x2-y2)初值 x=0.1;y=0.2;积分器中设置初始条件;f(u)中指定 Dx,Dy 的计算公式。运行这个仿真,scope 中可以看到两个变量的时程如下:在 WorkSpace 里可以得到 tout 和 yout,执行 plot(yout(:,1),yout(:,1)得到与 ode45 求 解相似的结果如下2.这部分解决一个使用 ode 求解器 dde23 没法求解的一类延迟微分方程(中性微分方程)。形如 x(t)=f(x(t-t1),x(t),x(t-t2),x(t-t3)这类方程。dde23 是无法求解的,但是可 以借助 simulink 仿真求解。看下面的这个例子。x(t)=A1*x(t-t1)+A2*x(t-t2)+B*u(t)t1=0.15;t2=0.5A1=-12 3 -3 A2=0.02 0 0 B=0106 -116 62 0 0.03 0 1207 -207 113 0
网站客服QQ:2055934822
金锄头文库版权所有
经营许可证:蜀ICP备13022795号 | 川公网安备 51140202000112号