常微分方程求解

本章知识要点•数值计算–常微分方程初值问题–常微分方程边值问题•MATLAB–微分方程求解常微分方程的相关函数•ode45ode23•bvp4c微分方程在化工模型中的应用•间歇反应器的计算•活塞流反应器的计算•全混流反应器的动态模拟•定态一维热传导问题•逆流壁冷式固定床反应器一维模型•固定床反应器的分散模型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解适度刚性ODE输出odeplotODE的输出时间序列图ode15s变阶法解刚性ODEodephas2ODE的二维相平面图ode23s低阶法解刚性ODEodephas3ODE的三维相空间图ode23tb低阶法解刚性ODEodeprint在Matlab指令窗显示结果odefile所谓的odefile实际上是一个Matlab函数文件，一般作为整个求解程序的一个子函数，表示ode求解问题Matlab提供了odefile的模板，采用typeodefile命令显示其详细内容，然后将其复制到脚本编辑窗口，在合适的位置填入所需内容一般而言，对于程序通用性要求不高的场合，只需将原有模型写成标准形式，然后“翻译”成Matlab语言即可odefile的编写规定1.ode文件的最简单格式必须有一个自变量t和函数y作为输入变量，一个y的导函数作为输出变量。其中自变量t不论在ode文件中是否使用都必须作为第一输入变量，y则必须作为第二输入变量，位置不能颠倒。2.可以向ode文件中传递参数，数目不受限制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];解算指令的使用方法调用格式：1.[T,Y]=ode45(@fun,TSPAN,Y0)2.[T,Y]=ode45(@fun,TSPAN,Y0,options)3.[T,Y]=ode45(@fun,TSPAN,Y0,options,P1,P2,…)4.[T,Y,TE,YE,IE]=ode45(@fun,TSPAN,Y0,options,P1,P2,…)说明：输出变量T为返回时间列向量；解矩阵Y的每一行对应于T的一个元素，列数与求解变量数相等。@fun为函数句柄，为根据待求解的ODE方程所编写的ode文件（odefile）；TSPAN＝[T0TFINAL]是微分系统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使用梯形规则适中ode23tbTR-BDF2（隐式Runge-Kutta法）低ode解算指令的选择(1)求解初值问题：1)0(2'yyxyy)10(x比较ode45和ode23的求解精度和速度1.根据常微分方程要求的求解精度与速度要求ode45和ode23的比较functionCha5demo1%Comparisonofresultsobtainedfromode45andode23solverformatlongy0=1;tic,[x1,y1]=ode45(@fun,[0,1],y0);t_ode45=toctic,[x2,y2]=ode23(@fun,[0,1],y0);t_ode23=tocplot(x1,y1,'b-',x2,y2,'m-.'),xlabel('x'),ylabel('y'),legend('ODE45','ODE23','location','Northwest')disp('ComparativeResultsatx=1:');fprintf('\nODE45\t\t\ty=%.8f\nODE23\t\t\ty=%.8f\nPrecisiveResult...=%.8f\n',y1(end),y2(end),1.7320508)%------------------------------------------------------------------------functionf=fun(x,y)f=y-2*x/y;ode解算指令的选择(2)2.根据常微分方程组是否为刚性方程如果系数矩阵A的特征值连乘积小于零，且绝对值最大和最小的特征值之比（刚性比）很大，则称此类方程为刚性方程1)0(,2)0(10099.9901.0212'221'1yyyyyyy00)()('yxyxbAyy•刚性方程在化学工程和自动控制领域的模型中比较常见。刚性比：100/0.01＝10000ode解算指令的选择(2)•刚性方程的物理意义：常微分方程组所描述的物理化学变化过程中包含了多个子过程，其变化速度相差非常大的数量级，于是常微分方程组含有快变和慢变分量。•常微分方程组数值积分的稳定步长受模值最大的特征值控制，即受快变量分量约束，特征值大则允许步长小；而过程趋于稳定的时间又由模值最小的特征值控制，特征值小则积分到稳定的时间则长。•Matlab提供了不同种类的刚性方程求解指令：ode15sode23sode23tode23tb，可根据实际情况选用functionCha5demo3figureode23s(@fun,[0,100],[0;1])holdon,ode45(@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];刚性常微分方程组求解解算指令的输出控制1.[T,Y]=ode45(@fun,TSPAN,Y0,options,P1,P2,…)2.sol=ode45(@fun,TSPAN,Y0)将解输出指结构体sol中；SXINT=deval(SOL,XINT)用于计算解在XINT处的值，XINT必须位于区间[SOL.x(1)SOL.x(end)]内。3.在无输出变量时，将调用默认的odeplot输出解的图形。4.除了以odeplot形式输出外，还可以以odephas2，和odephas3的形式输出解向量的二维和三维相平面图。5.采用以下语句options=odeset(‘outputfcn’,‘odephas2’)可以将输出方法改变为平面输出6.odeprint输出求解过程每一步的解解算指令的options选项1.RelTol－相对误差，它应用于解向量的所有分量。在每一步积分过程中，第i个分量误差e(i)满足：e(i)=max(RelTol*abs(y(i),AbsTol(i))。2.AbsTol－绝对误差，若是实数，则应用于解向量的所有分量，若是向量，则它的每一个元素应用于对应位置解向量元素。3.OutputFcn－可调用的输出函数名。每一步计算完后，这个函数将被调用输出结果，可以选择的值为：odeplot，odephas2，odephas3，odeprint。4.OutputSel－输出序列选择。指定解向量的哪个分量被传递给OutputFcn。5.MaxSetp－步长上界，缺省值为求解区间的1/10。6.InitialStep－初始步长，缺省时自动设置。7.采用odeset改变原有选项的值在间歇反应器中进行液相反应制备产物B，反应网络如图5-1所示。反应可在180~260℃的温度范围内进行，反应物X大量过剩，而C,D和E为副产物。各反应均为一级动力学关系：r＝－kC，式中间歇反应器中串联－平行复杂反应系统已知参数：k01＝5.78052×1010，k02＝3.92317×1012，k03＝1.64254×104，k04＝6.264×108，Ea1=124670，Ea2=150386，Ea3=77954，Ea4=111528。初始浓度：CA=1kmol/m3，其余物质浓度为0。已知是产物B收率最大的最优反应温度为224.6℃试计算1)在最优反应温度下各组分浓度随时间的动态变化；2)最优反应时间；3)输出产物D对反应物浓度A的关系图。RTEaekk0间歇反应器中串联－平行复杂反应系统数学模型AACkkdtdC)(21BABCkCkdtdC31CACCkCkdtdC42DBDCkCkdtdC53DCECkCkdtdC54间歇反应器中串联－平行复杂反应系统functionCha5demo4T=224.6+273.15;R=8.31434;k0=[5.78052E+103.92317E+121.64254E+46.264E+8];Ea=[12467015038677954111528];%Initialconcentration:C0(