工程数学作业第十次方健(124)

姓名：方健学号：652081701073专业：化学工程班级：化工1503问题：在MATLAB中使用ode函数求解常微分方程，并且改变初始值以及步长的变化观察结果的变化，并且总结有关于ode函数的调用格式，以及使用方法和特点，以及求解不同的方程适用范围。下面首先使用ode45函数进行求解微分函数，在matlab中使用helpode45学习该函数的基本格式，以及使用方法，ode45是适用于非刚性普通4-5阶法解微分方程，例如求解如下的微分方程：'2(0)1xyyyy创建函数fjfun1,将此函数保存在m文件中functionf=fjfun1(x,y)f=y-2*x/y;然后在matlab’中编写m文件如下：首先编写函数my111.mformatlongy0=1;[x1,y1]=ode45(@fjfun1,[0,1],y0);[x2,y2]=ode23(@fjfun1,[0,1],y0);plot(x1,y1,'b-',x2,y2,'m-.'),xlabel('x'),ylabel('y'),legend('ODE45','ODE23','location','Northwest')y1(end)y2(end)在commandwindow中运行结果如下：my111my111ans=1.732050816766531ans=1.732154879417795从上面的结果可看出，对于两个求解微分方程函数结果不同，oode45四五阶Runge-Kutta法结果的计算较高，ode23二三阶Runge-Kutta法低。下面研究关于微分方程组的解法对于这样的一个微分方程组：'112'213'3121230.78(0)0,(0)1,(0)1yyyyyyyyyyyy对于方程组而言，编写方程组的函数不同于一个微分方程，编写fjfun2.m文件如下：functiondy=fjfun2(t,y)dy=zeros(3,1);%%方程组为三个微分方程，生成三行一列的0向量dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.78*y(1)*y(2);end然后再编写求解函数my112.m文件如下：options=odeset('RelTol',1e-4,'AbsTol',[1e-41e-41e-5]);[t,y]=ode45(@fjfun2,[015],[011],options);plot(t,y(:,1),'-',t,y(:,2),'-.',t,y(:,3),'.')xlabel('timet'),00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8xyODE45ODE23ylabel('solutiony'),legend('y1','y2','y3')然后再在commandwindow中运行结果如下：my112下面改变微分方程的区间，观察变化options=odeset('RelTol',1e-4,'AbsTol',[1e-41e-41e-5]);[t,y]=ode45(@fjfun2,[0100],[011],options);plot(t,y(:,1),'-',t,y(:,2),'-.',t,y(:,3),'.')xlabel('timet'),ylabel('solutiony'),legend('y1','y2','y3')my112051015-1-0.8-0.6-0.4-0.200.20.40.60.81timetsolutionyy1y2y3当微分方程的刚性比很大的时候，可以得出如下结果'123'213'312123100*0.01*(0)0,(0)1,(0)1yyyyyyyyyyyyfunctiondy=fjfun2(t,y)dy=zeros(3,1);%·½³Ì×éΪÈý¸ö΢·Ö·½³Ì£¬Éú³ÉÈýÐÐÒ»ÁеÄ0ÏòÁ¿dy(1)=100*y(2)*y(3);dy(2)=-0.01*y(1)*y(3);dy(3)=-y(1)*y(2);endoptions=odeset('RelTol',1e-4,'AbsTol',[1e-41e-41e-5]);[t,y]=ode45(@fjfun2,[015],[011],options);plot(t,y(:,1),'-',t,y(:,2),'-.',t,y(:,3),'.')xlabel('timet'),ylabel('solutiony'),legend('y1','y2','y3')运行结果如下所示：0102030405060708090100-1.5-1-0.500.511.5timetsolutionyy1y2y3下面采用ode15s函数求解：options=odeset('RelTol',1e-4,'AbsTol',[1e-41e-41e-5]);[t,y]=ode15s(@fjfun2,[015],[011],options);plot(t,y(:,1),'-',t,y(:,2),'-.',t,y(:,3),'.')xlabel('timet'),ylabel('solutiony'),legend('y1','y2','y3')051015-10-8-6-4-20246810timetsolutionyy1y2y3结果分析与讨论：从上面的结果可以得出，明显对于ode15s的函数求解微分方程的刚性问题优势更加的明显。ode45基于显式4-5阶龙格库塔公式，其算法属于单步法；ode15s是一个变阶求解器，用的是多步法。对于很多问题，这些求解器都是可以使用的，尽管可能存在一些效率和精度方面的差异。但是，这些求解器并不是可以互相取代的，它们分别适用于不同的精度要求和问题的类型。Ode15s对于求解刚性微分方程十分的适用，但是ode45对于求解非刚性的问题十分不理想。因此对于不同的微分方程组应该根据每一个求解器的特点取选取相对应的求解器。051015-1.5-1-0.500.511.5timetsolutionyy1y2y3