常微分方程(组)的数值解法

2.3常微分方程（组）的数值解法知识要点常微分方程初值问题---ode45,0de23微分方程在化工模型中的应用•间歇反应器的计算•活塞流反应器的计算•全混流反应器的动态模拟•定态一维热传导问题•逆流壁冷式固定床反应器一维模型•固定床反应器的分散模型Matlab常微分方程求解问题分类初值问题：•定解附加条件在自变量的一端•一般形式为：•初值问题的数值解法一般采用步进法，如Runge-Kutta法边值问题：在自变量两端均给定附加条件一般形式：边值问题可能有解、也可能无解，可能有唯一解、也可能有无数解边值问题有3种基本解法•迭加法•打靶法•松弛法0)(),('yayyxfy21)(,)(),('ybyyayyxfyMatlab求解常微分方程初值问题方法将待求解微分方程（组）转化为标准形式，“翻译”成Matlab可以理解的语言编写odefile文件选择合适的解算指令求解问题根据求解问题的要求，设置解算指令的调用格式Matlab求解初值问题函数指令含义指令含义解算ode23普通2-3阶法解ODEodefileODE文件格式ode45普通4-5阶法解ODE选项odeset创建、更改ODE选项的设置ode113普通变阶法解ODEodeget读取ODE选项的设置ode23t解适度刚性ODEode15s变阶法解刚性ODEode23s低阶法解刚性ODEode23tb低阶法解刚性ODEodefileodefile是一个Matlab函数文件，一般作为整个求解程序的一个子函数，表示ode求解问题对于程序通用性要求不高的场合，只需将原有模型写成标准形式，然后“翻译”成Matlab语言即可odefile的编写---常微分方程functionf=fun(x,y)f=y-2*x/y;求解初值问题：1)0(2'yyxyy)10(x自变量在前，因变量在后ode输入函数输出变量为因变量导数的表达式0)(),('yayyxfy初值问题：1)0('2yyyy)10(xfunctionf=fun(x,y)f=y＋y^2;odefile的编写---常微分方程组1000,1)0(,0)0()1(3000'10')1(0001.0)1()1(04.0'2122142221211xyyyyyyyyyy常微分方程组与单个常微分方程求解方法相同，只需在编写odefile时将整个方程组作为一个向量输出。functionf=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的编写本例的难度：teytbytayt2cos)()')((2)2cos()(,2cos)(2ttbteetatt求解：y(0)=0，y'(0)=1，方程系数非线性可在odefile中定义方程高阶，非标准形式方程变形：令y1＝y；y2＝y’则原方程等价于：tebyayyyyt2cos122'22'1functionf=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的每一行对应于T的一个元素，列数与求解变量数相等。@fun为函数句柄，为根据待求解的ODE方程所编写的ode文件（odefile）；TSPAN＝[T0TFINAL]是微分系统y'＝F(t,y)的积分区间；Y0为初始条件常微分方程初值问题解算指令比较解算指令算法精度ode45四五阶Runge-Kutta法较高ode23二三阶Runge-Kutta法低ode113可变阶Adams-Bashforth-Moulton法ode15s基于数值差分的可变阶方法（BDFs，Gear）低～中ode23s二阶改进的Rosenbrock法低ode23t使用梯形规则适中ode23tbTR-BDF2（隐式Runge-Kutta法）低ode解算指令的选择(1)求解初值问题：1)0(2'yyxyy)10(x比较ode45和ode23的求解精度和速度1.根据常微分方程要求的求解精度与速度要求ode45和ode23的比较-1functionxODEclearallclcformatlongy0=1;[x1,y1]=ode45(@f,[0,1],y0);[x2,y2]=ode23(@f,[0,1],y0);plot(x1,y1,'k-',x2,y2,'b--')xlabel('x')ylabel('y')%------------------------------------------------------------------functiondydx=f(x,y)dydx=y-2*x/y;ode45和ode23的比较-2functionxODEsclearallclcy0=[01];[x1,y1]=ode45(@f,[0:0.2:1],y0);[x2,y2]=ode23(@f,[0:0.2:1],y0);disp('Resultsbyusingode45()：')disp('xy(1)y(2)')disp([x1y1]);disp('Resultsbyusingode23()：')disp('xy(1)y(2)')disp([x2y2]);%------------------------------------------------------------------functiondydx=f(x,y)%DefinesimultaneousODEequationsdydx=[3*y(1)+2*y(2);4*y(1)+y(2)];ode解算指令的选择(2)2.根据常微分方程组是否为刚性方程如果系数矩阵A的特征值连乘积小于零，且绝对值最大和最小的特征值之比（刚性比）很大，则称此类方程为刚性方程1)0(,2)0(10099.9901.0212'221'1yyyyyyy00)()('yxyxbAyy•刚性方程在化学工程和自动控制领域的模型中比较常见。刚性比：100/0.01＝10000ode解算指令的选择(2)•刚性方程的物理意义：常微分方程组所描述的物理化学变化过程中包含了多个子过程，其变化速度相差非常大的数量级，于是常微分方程组含有快变和慢变分量。•Matlab提供了不同种类的刚性方程求解指令：ode15sode23sode23tode23tb，可根据实际情况选用functiondemo1figureode23s(@fun,[0,100],[0;1])holdon,pauseode45(@fun,[0,100],[0;1])%--------------------------------------------------------------------------functionf=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];刚性常微分方程组求解间歇反应器中串联－平行复杂反应系统（例4-1）functionBatchReactorclearallclcT=224.6+273.15;%Reactortemperature,KelvinR=8.31434;%Gasconstant,kJ/kmolK%Arrheniusconstant,1/sk0=[5.78052E+103.92317E+121.64254E+46.264E+8];%Activationenergy,kJ/kmolEa=[12467015038677954111528];%初始浓度C0(i),kmol/m^3C0=[10000];tspan=[01e4];[t,C]=ode45(@MassEquations,tspan,C0,[],k0,Ea,R,T);%绘图plot(t,C(:,1),'r-',t,C(:,2),'k:',t,C(:,3),'b-.',t,C(:,4),'k--');xlabel('Time(s)');ylabel('Concentration(kmol/m^3)');legend('A','B','C','D')CBmax=max(C(:,2));%CBmax:themaximumconcentrationofB,kmol/m^3yBmax=CBmax/C0(1)%yBmax:themaximumyieldofBindex=find(C(:,2)==CBmax);t_opt=t(index)%t_opt:theoptimumbatchtime,s%------------------------------------------------------------------functiondCdt=MassEquations(t,C,k0,Ea,R,T)%Reactionrateconstants,1/sk=k0.*exp(-Ea/(R*T));k(5)=2.16667E-04;%Reactionrates,kmoles/m3srA=-(k(1)+k(2))*C(1);rB=k(1)*C(1)-k(3)*C(2);rC=k(2)*C(1)-k(4)*C(3);rD=k(3)*C(2)-k(5)*C(4);rE=k(4)*C(3)+k(5)*C(4);%MassbalancesdCdt=[rA;rB;rC;rD;rE];三个串联的CSTR等温反应器（例4-3）functionIsothermCSTRsclearallclcCA0=1.8;%kmol/m^3CA10=0.4;%kmol/m^3CA20=0.2;%kmol/m^3CA30=0.1;%kmol/m^3k=0.5;%1/mintau=2;stoptime=2.9;%min[t,y]=ode45(@Equations,[0stoptime],[CA10CA20CA30],[],k,CA0,tau);disp('Results:')disp('tCA1CA2CA3')disp([t,y])plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-')legend('CA_1','CA_2','CA_3')xlabel('Time(min)')ylabel('Concentration')%------------------------------------------------------------------functiondydt=Equations(t,y,k,CA0,tau)CA1=y(1);CA2=y(2);CA3=y(3);dCA1dt=(CA0-CA1)/tau-k*CA1;dCA2dt=(CA1-CA2)/tau-k*CA2;dCA3dt=(CA2-CA3)/tau-k*CA3;dydt=[dCA1dt;dCA2dt;dCA3dt];