常微分方程及其求解实例

第4章常微分方程及其求解,常微分方程（组）数值解,1初值问题,一阶常微分方程,一阶常微分方程,高阶常微分方程,(1)Euler法,functiontout,yout=euler(ypfun,t0,tfinal,y0,tol,trace)pow=1/3;ifnargincomet3(x(:,1),x(:,2),x(:,3),描述微分方程是常微分方程初值问题数值求解的关键。f1=inline(-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);.-x(1)*x(2)+28*x(2)-x(3),t,x);t_final=100;x0=0;0;1e-10;t,x=ode45(f1,0,t_final,x0);plot(t,x),figure;plot3(x(:,1),x(:,2),x(:,3);axis(1042-2020-2025);得出完全一致的结果。,sol:结构体数据，sol.x:时间向量t,sol.y:状态向量。,（3）延迟微分方程求解,例：,编写函数：functiondx=c7exdde(t,x,z)xlag1=z(:,1);%第一列表示提取xlag2=z(:,2);dx=1-3*x(1)-xlag1(2)-0.2*xlag2(1)3-xlag2(1);x(3);4*x(1)-2*x(2)-3*x(3);历史数据函数：functionS=c7exhist(t)S=zeros(3,1);,求解：lags=10.5;tx=dde23(c7exdde,lags,zeros(3,1),0,10);plot(tx.x,tx.y(2,:)与ode45()等返回的x矩阵不一样，它是按行排列的。,（4）刚性问题,A的特征值为,一类特殊的常微分方程，其中一些解变化缓慢，另一些变化快，且相差悬殊，这类方程常常称为刚性方程。,A=-0.0100-99.99000-100.0000,eig(a)=-0.0100-100.0000,刚性问题的积分函数t,y=ode15s(fun,tspan,y0)t,y=ode15s(fun,tspan,y0,options,p1,p2)t,y=ode23s(fun,tspan,y0)t,y=ode23s(fun,tspan,y0,options,p1,p2),functionStiffODEsclearall;clcy0=0;1;%x1,y1=ode45(f,0100,y0)x2,y2=ode23s(f,0100,y0)%-functiondydx=f(x,y)%DefinesimultaneousODEequationsdy1dx=-0.01*y(1)-99.99*y(2);dy2dx=-100*y(2);dydx=dy1dx;dy2dx;,functiondy=c7exstf2(t,y)dy=0.04*(1-y(1)-(1-y(2)*y(1)+0.0001*(1-y(2)2;-104*y(1)+3000*(1-y(2)2;tict2,y2=ode45(c7exstf2,0,100,0;1);toclength(t2),plot(t2,y2),elapsed_time=229.4700ans=356941,用ode15s()代替ode45()opt=odeset;opt.RelTol=1e-6;tict1,y1=ode15s(c7exstf2,0,100,0;1,opt);toclength(t1),plot(t1,y1),（5）隐式微分方程求解,隐式微分方程:不能转化为显式常微分方程组的方程,functiondx=c7ximp(t,y)A=sin(y(1)cos(y(2);-cos(y(2)sin(y(1);B=1-y(1);-y(2);dx=inv(A)*B;opt=odeset;opt.RelTol=1e-6;t,y=ode45(c7ximp,0,10,0;0,opt);plot(t,y),数值解法:（1）、打靶法（2）、有限差分法（3）、配置法,2边值问题(BVP),（1）打靶法,先将边值问题看作一个初值问题，通过一系列的试验来估计第二个初始条件，判断第二个边界条件是否满足，如果不满足，换一个初值再试，重复该过程直到末端点的条件满足为止。,（2）有限差分法,用有限差分格式近似代替ODEs中的导数项,从而将ODEs方程离散为线性或非线性代数方程组,然后求解此代数方程组.,functiont,y=shooting(f1,f2,tspan,x0f,varargin)t0=tspan(1);tfinal=tspan(2);ga=x0f(1);gb=x0f(2);t,y1=ode45(f1,tspan,1;0,varargin);t,y2=ode45(f1,tspan,0;1,varargin);t,yp=ode45(f2,tspan,0;0,varargin);m=(gb-ga*y1(end,1)-yp(end,1)/y2(end,1);t,y=ode45(f2,tspan,ga;m,varargin);,functionxdot=c7fun1(x,y)xdot=y(2);-(1-x)*y(1)-(1+x)*y(2);functionxdot=c7fun2(x,y)xdot=y(2);-(1-x)*y(1)-(1+x)*y(2)+(1+x2);t,y=shooting(c7fun1,c7fun2,0,1,1;4);plot(t,y),有限差分法,functionx,y=fdiff(funcs,tspan,x0f,n)t0=tspan(1);tfinal=tspan(2);ga=x0f(1);gb=x0f(2);h=(tfinal-t0)/n;fori=1:n,x(i)=t0+h*(i-1);end,x0=x(1:n-1);t=-2+h2*feval(funcs,x0,2);tmp=feval(funcs,x0,1);v=1+h*tmp/2;w=1-h*tmp/2;b=h2*feval(funcs,x0,3);b(1)=b(1)-w(1)*ga;b(n-1)=b(n-1)-v(n-1)*gb;b=b;A=diag(t);fori=1:n-2,A(i,i+1)=v(i);A(i+1,i)=w(i+1);endy=inv(A)*b;x=xtfinal;y=ga;y;gb;,functiony=c7fun5(x,key)switchkeycase1,y=1+x;case2,y=1-x;otherwise,y=1+x.2;endt,y=fdiff(c7fun5,0,1,1,4,50);plot(t,y),（3）配置法MATLAB中提供的ODE-BVPs求解器，采用配置法。调用函数bvpinit生成在区间内的给定初始网格x上的初始猜测解调用函数bvp4c利用配置法得到解sol调用函数deval计算在区间内任意点xint处的解,函数bvpinitsolinit=bvpinit(x,yinit)%x为初始网格，yinit为猜测解,函数bvp4csol=bvp4c(ODEfun,Bcfun,solinit)sol=bvp4c(ODEfun,Bcfun,solinit,options,p1,p2),函数devalyint=deval(sol,xint),例1:,functionxBVP2clearallclca=0;b=1;%solutionisobtainedusinganinitialguessofy1(x)=0,y2(x)=0solinit=bvpinit(linspace(a,b,10),00);sol=bvp4c(ODEfun,BCfun,solinit);x=0:0.1:1;y=deval(sol,x);%-functiondydx=ODEfun(x,y)dydx=y(2);(1+x2)*y(1)+1;%-functionbc=BCfun(ya,yb)bc=ya(1)-1;yb(1)-3;,例2:,functionxBVP2clearallclca=0;b=1;%solutionisobtainedusinganinitialguessofy1(x)=0,y2(x)=0solinit=bvpinit(linspace(a,b,10),00);sol=bvp4c(ODEfun,BCfun,solinit);x=0:0.1:1;y=deval(sol,x);%-functiondydx=ODEfun(x,y)dydx=y(2);1+x2-(1+x)*y(2)-(1-x)*y(1);%-functionbc=BCfun(ya,yb)bc=ya(1)-1;yb(1)-4;,tica=0;b=1;solinit=bvpinit(linspace(a,b,10),00);sol=bvp4c(ODEfun,BCfun,solinit);x=0:0.02:1;y2=deval(sol,x);t0=toctict0,y1=shooting(c7fun1,c7fun2,0,1,1;4);t1=toctict1,y2=fdiff(c7fun5,0,1,1,4,50);t2=tocplot(x,y2(1,:),*,t0,y1(:,1),o,t1,y2,r),IVP应用1:间歇反应器连串-平行复杂反应系统,在间歇反应器中进行液相反应制备产物B,反应网络如下,反应可在180260度温度范围内进行,反应物X过剩,而C、D和E为副产物。各反应均为一级动力学关系：r=-kC，式中k=k0e-Ea/RT。已知参数：k01=5.780521010,k02=3.923171012,k03=1.64254104,Ea1=2.1667104,Ea2=150386,Ea3=77954,Ea4=111528.初始浓度为:t=0,CA=1kmol/m3,CC=CB=CD=CE=0kmol/m3.用最优化方法求得产物B收率最大的最优反应温度为224.6度.试计算(1)在此最优反应温度下各组份浓度随时间的动态变化;(2)最优反应时间.,clearall;clcT=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=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:themaxiconcenofByBmax=CBmax/C0(1)%yBmax:themaxiyieldofBindex=find(C(:,2)=CBmax);topt=t(index)%t_opt:theoptimumbatchtime,sfunctiondCdt=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;,IVP应用2:三釜串联连续搅拌等温反应器,在三个串联的CSTR等温反应器中，进行简单一级反应：AB已知条件为:CA0=1.8,CA10=0.4,CA20=0.2,CA30=0.1,(kmol/m3),k=0.5/min,=2min.试求解三个反应器中组分A浓度随时间的变化规律.,k,clearall;clcCA0=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(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;,IVP应用3：连续搅拌釜反应器例4.2,functionCSTRclearall；clcglobalERk0FCA0VT0UATJrhoCpHrHGHLtauSTOPTIME=5000;F=1.0e-8;%液体体积流量,m3/sCA0=5.0;%进料浓度,kmol/m3T0=300;%进料温度,KTJ=305;%夹套温度,KV=2.0e-6;%反应器体积,m3rho=1000.0;%液体密度,kg/m3Cp=4.187;%液体热容,kJ/kgKUA=5.68e-6;%传热系数与传热面积的乘积,kJ/KsHr=-4.19e4;%反应热,kJ/kmolk0=8.03E+12;%指前因子,1/sE=9.42e+4;%活化能,kJ/kmolR=8.317;%气体常数,kJ/kmolKtau=V/F;,%绘制Qg(Qr)T的关系图以便观察稳态点%-T=(290:380);k=RateConstant(T)CA=CA0./(1+tau*k);%STEADYSTATEMASSBALANCEQg,Qr=Heat(k,CA,T);figureplot(T,Qg,b-,T,Qr,k-);xlabel(T(K);ylabel(Q_g,Q_r(kJ/s);legend(Q_g,Q_r),%-functionk=RateConstant(T)globalk0ER%以下两式由式k=k0exp(-E/R/T)变换而成，通过此变换以避免溢出。temp=k0*exp(-E/(R*273);%Rateconstantat273K,1/sk=temp*exp(-(E/R)*(1./T-1/273);%Thisformavoidsmathoverflow%-functionQg,Qr=Heat(k,CA,T)globalFT0VUATJrhoCpHrQr=F*rho*Cp*(T-T0)+UA*(T-TJ);%移走的热量,kJ/sQg=k.*CA*V*(-Hr);%产生的热量,kJ/s,%动态计算tspan=0STOPTIME;CAI=linspace(0,5,2)%Initialtemp,KTI=linspace(300,340,10)%Initialconc.,kmol/m3%绘制CSTR反应器的相平面图figureholdonfori=1:length(CAI)forj=1:length(TI)t,y=ode45(DynamicModel,tspan,CAI(i),TI(j);CA=y(:,1);T=y(:,2);xA=1-CA/CA0;plot(T,xA);endendxlabel(T(K);ylabel(x_A);holdoff,%-functiondydt=DynamicModel(t,y)%DynamicModelglobalFCA0VT0UATJrhoCpHrtauCA=y(1);T=y(2);k=RateConstant(T);%Dynamicmassbalance,firstorderkineticsdCAdt=(CA0-CA)/tau-k*CA;%DynamicenergybalanceQg,Qr=Heat(k,CA,T);dTdt=(Qg-Qr)/(V*rho*Cp);%热量平衡,K/sdydt=dCAdt;dTdt;,functionNonIsothermTRclearallclcglobalCtAcrhoBCpH1H3UdtTJCpH1H3UTJrhogusE1E2E3RyO2L=1;%反应管长,mG=4684;%表观质量流速,kg/m2hrrhoB=1300;%催化剂堆积密度,kg/m3Mm=29.48;%气体的平均分子量,kg/kmolyA0=0.00924;%入口处邻二甲苯的摩尔分率yO2=0.208;%氧的摩尔分率(恒为常数)Cp=1.047;%比热,kJ/kmolKU=345.686;%传热系数,kJ/m2hrKP=101.325;%假设总压为定值=1atm=101.325kJdt=0.0254;%管径,mTJ=400;%冷却温度,KT0=700;%物料进口温度(初始温度),K,IVP应用4：管式反应器（等温、非等温）例4.5固定床一维稳态拟均相反应器,H1=-1.285e+6;%反应AB的反应热,kJ/kmolH3=-4.564e+6;%反应AC的反应热,kJ/kmol%活化能,kJ/kmolE1=1.1304e5;E2=1.315e5;E3=1.197e5;R=8.314;%理想气体常数,kJ/kmolKAc=pi*(dt/2)2;%反应管的横截面积,m2Ft=G*Ac/Mm;%总摩尔流率,moles/hrus=3600;%线速度，m/hrrhog=G/us;Ct=Ft/(Ac*us);FA0=yA0*Ft;%A的进料摩尔流率,kmol/hrCA0=FA0/(Ac*us);CB0=0;%FB0=0CC0=0;%FC0=0,z,y=ode45(Equations,0L,CA0CB0CC0T0)CA=y(:,1);CB=y(:,2);CC=y(:,3);xA=(CA0-CA)./CA0;%A的转化率xB=CB(2:end)./(CA0-CA(2:end);%生成的B/反应的AxB=0;xBxC=CC(2:end)./(CA0-CA(2:end);%生成的C/反应的AxC=0;xC,%图形输出plot(z,y(:,4)%温度分布xlabel(z)ylabel(T(K)figureplot(z,xA,r-)%转化率分布xlabel(z)ylabel(x_A)figureplot(z,CA,r-,z,CB,k-,z,CC,b:)%浓度分布xlabel(z)ylabel(C_A,C_B,C_C)legend(C_A,C_B,C_C),functiondydz=Equations(z,y)%模型方程组globalyO2CtAcrhoBCpH1H3UdtTJCpH1H3UTJrhogusCA=y(1);CB=y(2);CC=y(3);T=y(4);%摩尔分率yA=CA/Ct;yB=CB/Ct;yC=CC/Ct;%反应速度rA,rB,rC,k1,k2,k3=Rates(yA,yB,yC,T);%物料平衡dCAdz=rhoB*rA/us;dCBdz=rhoB*rB/us;dCCdz=rhoB*rC/us;%热量衡算dTdz=(rhoB*(-H1*k1-H3*k3)*yA*yO2-4*U*(T-TJ)/dt)/(us*rhog*Cp);dydz=dCAdz;dCBdz;dCCdz;dTdz;,%-functionrA,rB,rC,k1,k2,k3=Rates(yA,yB,yC,T)%反应动力学globalE1E2E3RyO2%速度常数,kmol/kgcatalysthrk1=exp(-E1/(R*T)+19.837);k2=exp(-E2/(R*T)+20.86);k3=exp(-E3/(R*T)+18.97);%反应速度,kmol/kgcatalysthrrA=-(k1+k3)*yA*yO2;%A的总反应速度rB=k1*yA*yO2-k2*yB*yO2;%B的净生成速率rC=k2*yB*yO2+k3*yA*yO2;%C的总生成速率,functionMembraneRe