Matlab与化工数值计算-第5讲常微分方程数值解

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第五讲常微分方程数值解,化工学院软件应用教科组,2006-10,本章知识要点,数值计算,常微分方程初值问题,常微分方程边值问题,MATLAB,微分方程求解常微分方程的相关函数,ode45 ode23,bvp4c,微分方程在化工模型中的应用,间歇反应器的计算,活塞流反应器的计算,全混流反应器的动态模拟,定态一维热传导问题,逆流壁冷式固定床反应器一维模型,固定床反应器的分散模型,Matlab,常微分方程求解问题分类,初值问题：,定解附加条件在自变量的一端,一般形式为：,初值问题的数值解法一般采用步进法，如,Runge-Kutta,法,边值问题：,在自变量两端均给定附加条件,一般形式：,边值问题可能有解、也可能无解，可能有唯一解、也可能有无数解,边值问题有,3,种基本解法,迭加法,打靶法,松弛法,Matlab,求解常微分方程初值问题方法,将待求解转化为标准形式，并“翻译”成,Matlab,可以理解的语言，即编写,odefile,文件,选择合适的解算指令求解问题,根据求解问题的要求，设置解算指令的调用格式,Matlab求解初值问题,函数,指令,含义,指令,含义,解,算,ode23,普通,2-3,阶法解,ODE,odefile,ODE,文件格式,ode45,普通,4-5,阶法解,ODE,选,项,odeset,创建、更改,ODE,选项的设置,ode113,普通变阶法解,ODE,odeget,读取,ODE,选项的设置,ode23t,解适度刚性,ODE,输,出,odeplot,ODE,的输出时间序列图,ode15s,变阶法解刚性,ODE,odephas2,ODE,的二维相平面图,ode23s,低阶法解刚性,ODE,odephas3,ODE,的三维相空间图,ode23tb,低阶法解刚性,ODE,odeprint,在,Matlab,指令窗显示结果,odefile,所谓的,odefile,实际上是一个,Matlab,函数文件，一般作为整个求解程序的一个子函数，表示,ode,求解问题,Matlab,提供了,odefile,的模板，采用,type odefile,命令显示其详细内容，然后将其复制到脚本编辑窗口，在合适的位置填入所需内容,一般而言，对于程序通用性要求不高的场合，只需将原有模型写成标准形式，然后“翻译”成,Matlab,语言即可,odefile,的编写规定,ode,文件的最简单格式必须有一个自变量,t,和函数,y,作为输入变量，一个,y,的导函数作为输出变量。其中自变量,t,不论在,ode,文件中是否使用都必须作为第一输入变量，,y,则必须作为第二输入变量，位置不能颠倒。,可以向,ode,文件中传递参数，数目不受限制,odefile,的编写,function f=fun(x,y),f=y-2*x/y;,求解初值问题：,自变量在前，因变量在后,ode,输入函数,输出变量为因变量导数的表达式,初值问题：,function f=fun(x,y),f=y,y2;,常微分方程组,odefile,的编写,常微分方程组与单个常微分方程求解方法相同，只需在编写,odefile,时将整个方程组作为一个向量输出。,function f=fun(x,y),dy1dx=0.04*(1-y(1)-(1-y(2).*y(1)+0.0001*(1-y(2).2;,dy2dx=-1e4*dy1dx+3000*(1-y(2).2;,f=dy1dx;dy2dx;,高阶微分方程,odefile,的编写,本例的难度：,求解：,y(0)=0,，,y(0)=1,，,方程系数非线性,可在,odefile,中定义,方程高阶，非标准形式,方程变形：令,y1,y,；,y2,y,则原方程等价于：,function f=fun(t,y),a=-exp(-t)+cos(2*pi*t)*exp(-2*t);,b=cos(2*pi*t);,f=y(2),-a*y(2)2-b*y(1)+exp(t)*b;,解算指令的使用方法,调用格式：,T,Y=ode45(fun,TSPAN,Y0),T,Y=ode45(fun,TSPAN,Y0,options),T,Y=ode45(fun,TSPAN,Y0,options,P1,P2,),T,Y,TE,YE,IE=ode45(fun,TSPAN,Y0,options,P1,P2,),说明：,输出变量,T,为返回时间列向量；解矩阵,Y,的每一行对应于,T,的一个元素，列数与求解变量数相等。,fun,为函数句柄，为根据待求解的,ODE,方程所编写的,ode,文件（,odefile,）；,TSPAN,T0 TFINAL,是微分系统,y,F(t,y),的积分区间；,Y0,为初始条件,options,用于设置一些可选的参数值，缺省时，相对于第一种调用格式。,options,中可以设置的参数参见,odeset,P1,，,P2,，,的作用是传递附加参数,P1,，,P2,，,到,ode,文件。当,options,缺省时，应在相应位置保留,，以便正确传递参数。,常微分方程初值问题解算指令比较,解算指令,算法,精度,ode45,四五阶,Runge-Kutta,法,较高,ode23,二三阶,Runge-Kutta,法,低,ode113,可变阶,Adams-Bashforth-Moulton,法,ode15s,基于数值差分的可变阶方法（,BDFs,，,Gear,）,低中,ode23s,二阶改进的,Rosenbrock,法,低,ode23t,使用梯形规则,适中,ode23tb,TR-BDF2,（隐式,Runge-Kutta,法）,低,ode,解算指令的选择,(1),求解初值问题：,比较,ode45,和,ode23,的求解精度和速度,1.,根据常微分方程要求的求解精度与速度要求,ode45,和,ode23,的比较,function Cha5demo1,%Comparison of results obtained from ode45 and ode23 solver,format long,y0=1;,tic,x1,y1=ode45(fun,0,1,y0);t_ode45=toc,tic,x2,y2=ode23(fun,0,1,y0);t_ode23=toc,plot(x1,y1,b-,x2,y2,m-.),xlabel(x),ylabel(y),legend(ODE45,ODE23,location,Northwest),disp(Comparative Results at x=1:);,fprintf(nODE45ttt y=%.8fnODE23ttt y=%.8fnPrecisive Result.=%.8fn,y1(end),y2(end),1.7320508),%-,function f=fun(x,y),f=y-2*x/y;,ode,解算指令的选择,(2),2.,根据常微分方程组是否为刚性方程,如果系数矩阵,A,的特征值连乘积小于零，且绝对值最大和最小的特征值之比（刚性比）很大，则称此类方程为刚性方程,刚性方程在化学工程和自动控制领域的模型中比较常见。,刚性比：,100/0.01,10000,ode,解算指令的选择,(2),刚性方程的物理意义：常微分方程组所描述的物理化学变化过程中包含了多个子过程，其变化速度相差非常大的数量级，于是常微分方程组含有快变和慢变分量。,常微分方程组数值积分的稳定步长受模值最大的特征值控制，即受快变量分量约束，特征值大则允许步长小；而过程趋于稳定的时间又由模值最小的特征值控制，特征值小则积分到稳定的时间则长。,Matlab,提供了不同种类的刚性方程求解指令：,ode15s ode23s ode23t ode23tb,，可根据实际情况选用,function Cha5demo3,figure,ode23s(fun,0,100,0;1),hold on,ode45(fun,0,100,0;1),%-,function f=fun(x,y),dy1dx=0.04*(1-y(1)-(1-y(2).*y(1)+0.0001*(1-y(2).2;,dy2dx=-1e4*dy1dx+3000*(1-y(2).2;,f=dy1dx;dy2dx;,刚性常微分方程组求解,解算指令的输出控制,T,Y=ode45(fun,TSPAN,Y0,options,P1,P2,),sol=ode45(fun,TSPAN,Y0),将解输出指结构体,sol,中；,SXINT=deval(SOL,XINT),用于计算解在,XINT,处的值，,XINT,必须位于区间,SOL.x(1)SOL.x(end),内。,在无输出变量时，将调用默认的,odeplot,输出解的图形。,除了以,odeplot,形式输出外，还可以以,odephas2,，和,odephas3,的形式输出解向量的二维和三维相平面图。,采用以下语句,options=odeset(outputfcn,odephas2),可以将输出方法改变为平面输出,odeprint,输出求解过程每一步的解,解算指令的,options,选项,RelTol,相对误差，它应用于解向量的所有分量。在每一步积分过程中，第,i,个分量误差,e(i),满足：,e(i)=max(RelTol*abs(y(i),AbsTol(i),。,AbsTol,绝对误差，若是实数，则应用于解向量的所有分量，若是向量，则它的每一个元素应用于对应位置解向量元素。,OutputFcn,可调用的输出函数名。每一步计算完后，这个函数将被调用输出结果，可以选择的值为：,odeplot,，,odephas2,，,odephas3,，,odeprint,。,OutputSel,输出序列选择。指定解向量的哪个分量被传递给,OutputFcn,。,MaxSetp,步长上界，缺省值为求解区间的,1/10,。,InitialStep,初始步长，缺省时自动设置。,采用,odeset,改变原有选项的值,在间歇反应器中进行液相反应制备产物,B,，反应网络如图,5-1,所示。反应可在,180260,的温度范围内进行，反应物,X,大量过剩，而,C,D,和,E,为副产物。各反应均为一级动力学关系：,r,kC,，式中,间歇反应器中串联平行复杂反应系统,已知参数：,k,01,5.7805210,10,，,k,02,3.9231710,12,，,k,03,1.6425410,4,，,k,04,6.26410,8,，,Ea,1,=124670,，,Ea,2,=150386,，,Ea,3,=77954,，,Ea,4,=111528,。初始浓度：,C,A,=1kmol/m,3,，其余物质浓度为,0,。已知是产物,B,收率最大的最优反应温度为,224.6,试计算,1),在最优反应温度下各组分浓度随时间的动态变化；,2),最优反应时间；,3),输出产物,D,对反应物浓度,A,的关系图。,间歇反应器中串联平行复杂反应系统,数学模型,间歇反应器中串联平行复杂反应系统,function Cha5demo4,T=224.6+273.15;R=8.31434;k0=5.78052E+10 3.92317E+12 1.64254E+4 6.264E+8;,Ea=124670 150386 77954 111528;,%Initial concentration:C0(i),kmol/m3,C0=1 0 0 0 0;tspan=0 1e4;,opt=odeset(reltol,1e-4,outputfcn,odephas2,outputsel,1;4),t,C=ode45(MassEquations,tspan,C0,opt,k0,Ea,R,T),%plot,plot(t,C(:,1),r-,t,C(:,2),k:,t,C(:,3),b-.,t,C(:,4),k-);,xlabel(Time(s);,ylabel(Concentration(