如何使用ODE45

1如何使用ODE45ODE45ODE45ODE45在MATLAB中ode23,ode45,ode113,ode15s,ode23s,ode23t,ode23tb等函数都是用来解决常微分方程的初值问题。根据MATLAB的帮助文档，应优先尝试使用ODE45求解器。基本语法：[T,Y]=solver(odefun,tspan,y0)[T,Y]=solver(odefun,tspan,y0,options)其中的求解器为ode45,ode23,ode113,ode15s,ode23s,ode23t,orode23tb之一。相关参数介绍如下：参数名称参数说明odefun用于存放待求解的方程的m文件名，方程必须用y’=f(t,y)的形式存放tspan指定自变量范围的向量，通常用[t0tf]指定y0函数的边界条件，即y0=y(t0)，对于方程组，y0也可以是向量options设置求解的相关选项，可以使用odeset函数创建选项下面看一个简单的例子：求解方程：'1,(1)1yyyt−==在14t≤≤的解事实上这个微分方程可以采用一般的公式求解，其理论解为(ln1)ytt=+，下面用ODE45来求解。首先要准备好一个M文件，我们命名为func.m，内容如下：functiondy=func(t,y)dy=y/t+1;%±ØÐëʹÓÃy'=f(t,y)µÄÐÎʽend然后在MATLAB命令窗口中输入：[t,y]=ode45(@func,[1,4],1);或者：[t,y]=ode45(‘func’,[1,4],1);运行后得到的t和y都是45*1的向量，可以使用函数plot(t,y)得到函数图：可以验证和实际的函数图形是很接近的，方法如下：2新建一个test.m文件，内容如下：[t,y]=ode45(@func,[1,4],1);plot(t,y,['r','*'])x=1:0.01:4;z=x.*(log(x)+1);holdonplot(x,z);grid;蓝色的线条是真实函数的曲线，而红色的点是数值计算的结果，可见两者符合的很好。下面是一个比较复杂一点的方程，所谓的复杂是说用理论方法很难求解的：求微分方程232,(1)2yytytyt′=+++=−在14t≤≤时的数值解首先要将方程改写为：2232yyytttt′=+++然后编写一个M文件，例如yy.m：functiondy=yy(t,y)dy=(y+3*t+y*t+2)/t/t;end然后在命令窗口中输入：[t,y]=ode45(@yy,[1,4],-2);即可求出y在t取不同值的对应值。另外，也可以求微分方程组的问题，例如Matlab帮助文档中提供的一个例子：3求方程组'123'213'3220.51yyyyyyyyy==−=−在满足初值条件123(0)0,(0)1,(0)1yyy===时的数值解和前面一样，写一个名为rigid的M文件：functiondy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);end这里将演示一下odeset的用法。在命令窗口中输入：options=odeset('RelTol',1e-4,'AbsTol',[1e-41e-41e-5]);[t,y]=ode45(@rigid,[012],[011],options);计算，然后绘图有：plot(t,y(:,1),'-',t,y(:,2),'-.',t,y(:,3),'.')grid对于二阶微分方程，基本原理是一样的，但是M文件却有很大的不同。二阶常微分方程的一般形式为：0001''()'()(),(),'()yptyqtygtytyyty++===下面的一个例子将做出说明。4求微分方程''''3sin2,(0)1,(0)1tytyeytyy+−===−在02t≤≤时的数值解。这个时候需要做一点技术上的处理：12'xyxy==这样，原方程可改写为：122123sin(2)''txxxtxext==−++现在变为和一阶常微分方程组了，编写一个M文件：functionxp=order2(t,x)xp=zeros(2,1);%注意这里是方程组xp(1)=x(2);xp(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);end那么在MATLAB中又该如何使用呢ODE命令呢？一般的格式是这样的：[t,x]=ode45('F',[t0,tf],[x10,x20]);F—指向待求解的函数文件t0,tf—t的初值和终值x10—x1的初始值x20—x2的初始值那么上例在MATLAB命令窗口中的代码如下：[t,x]=ode45('order2',[0,2],[1,-1]);执行后可得到t和x的值。如果需要绘制t-y函数图，可以使用命令：plot(t,x(:,1))或者想同时显示y和y’，那么可以使用：plot(t,x)但是这样不好区分各条曲线代表的函数，更好一点，指定颜色和绘图样式：plot(t,x(:,1),['r','+'],t,x(:,2),['g','*'])显示效果如下图，其中红色的是y(t)，绿色的是y’(t)。5