常微分方程组的数值解法

2.3常微分方程（组）的数值解法,知识要点,常微分方程初值问题-ode45,0de23,微分方程在化工模型中的应用,间歇反应器的计算活塞流反应器的计算全混流反应器的动态模拟定态一维热传导问题逆流壁冷式固定床反应器一维模型固定床反应器的分散模型,Matlab常微分方程求解问题分类,初值问题：定解附加条件在自变量的一端一般形式为：初值问题的数值解法一般采用步进法，如Runge-Kutta法,边值问题：在自变量两端均给定附加条件一般形式：边值问题可能有解、也可能无解，可能有唯一解、也可能有无数解边值问题有3种基本解法迭加法打靶法松弛法,Matlab求解常微分方程初值问题方法,将待求解微分方程（组）转化为标准形式，“翻译”成Matlab可以理解的语言编写odefile文件,选择合适的解算指令求解问题,根据求解问题的要求，设置解算指令的调用格式,Matlab求解初值问题函数,odefile,odefile是一个Matlab函数文件，一般作为整个求解程序的一个子函数，表示ode求解问题对于程序通用性要求不高的场合，只需将原有模型写成标准形式，然后“翻译”成Matlab语言即可,odefile的编写-常微分方程,functionf=fun(x,y)f=y-2*x/y;,求解初值问题：,自变量在前，因变量在后,ode输入函数,输出变量为因变量导数的表达式,初值问题：,functionf=fun(x,y)f=yy2;,odefile的编写-常微分方程组,常微分方程组与单个常微分方程求解方法相同，只需在编写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的编写,本例的难度：,求解：,y(0)=0，y(0)=1，,方程系数非线性,可在odefile中定义,方程高阶，非标准形式,方程变形：令y1y；y2y则原方程等价于：,functionf=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）；TSPANT0TFINAL是微分系统yF(t,y)的积分区间；Y0为初始条件,常微分方程初值问题解算指令比较,ode解算指令的选择(1),求解初值问题：,比较ode45和ode23的求解精度和速度,1.根据常微分方程要求的求解精度与速度要求,ode45和ode23的比较-1,functionxODEclearallclcformatlongy0=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的比较-2,functionxODEsclearallclcy0=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的特征值连乘积小于零，且绝对值最大和最小的特征值之比（刚性比）很大，则称此类方程为刚性方程,刚性方程在化学工程和自动控制领域的模型中比较常见。,刚性比：100/0.0110000,ode解算指令的选择(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/m3C0=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/m3);legend(A,B,C,D)CBmax=max(C(:,2);%CBmax:themaximumconcentrationofB,kmol/m3yBmax=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/m3CA10=0.4;%kmol/m3CA20=0.2;%kmol/m3CA30=0.1;%kmol/m3k=0.5;%1/mintau=2;stoptime=2.9;%mint,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(Concentrat