常微分方程求解优秀课件.ppt

《常微分方程求解优秀课件.ppt》由会员分享，可在线阅读，更多相关《常微分方程求解优秀课件.ppt（38页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、常微分方程求解第1页，本讲稿共38页本章知识要点数值计算常微分方程初值问题常微分方程边值问题MATLAB微分方程求解常微分方程的相关函数ode45 ode23bvp4c第2页，本讲稿共38页微分方程在化工模型中的应用间歇反应器的计算活塞流反应器的计算全混流反应器的动态模拟定态一维热传导问题逆流壁冷式固定床反应器一维模型固定床反应器的分散模型第3页，本讲稿共38页Matlab常微分方程求解问题分类初值问题：定解附加条件在自变量的一端 一般形式为：初值问题的数值解法一般采用步进法，如Runge-Kutta法边值问题：在自变量两端均给定附加条件一般形式：边值问题可能有解、也可能无解，可能有唯一解、也

2、可能有无数解边值问题有3种基本解法迭加法打靶法松弛法第4页，本讲稿共38页Matlab求解常微分方程初值问题方法将待求解转化为标准形式，并“翻译”成Matlab可以理解的语言，即编写odefile文件选择合适的解算指令求解问题根据求解问题的要求，设置解算指令的调用格式第5页，本讲稿共38页Matlab求解初值问题函数指令含义指令含义解算ode23普通2-3阶法解ODEodefileODE文件格式ode45普通4-5阶法解ODE选项odeset创建、更改ODE选项的设置ode113普通变阶法解ODEodeget读取ODE选项的设置ode23t解适度刚性ODE输出odeplotODE的输出时间序列

3、图ode15s变阶法解刚性ODEodephas2ODE的二维相平面图ode23s低阶法解刚性ODEodephas3ODE的三维相空间图ode23tb低阶法解刚性ODEodeprint在Matlab指令窗显示结果第6页，本讲稿共38页odefileZ所谓的odefile实际上是一个Matlab函数文件，一般作为整个求解程序的一个子函数，表示ode求解问题ZMatlab提供了odefile的模板，采用type odefile命令显示其详细内容，然后将其复制到脚本编辑窗口，在合适的位置填入所需内容Z一般而言，对于程序通用性要求不高的场合，只需将原有模型写成标准形式，然后“翻译”成Matlab语言即可

4、第7页，本讲稿共38页odefile的编写规定1.ode文件的最简单格式必须有一个自变量t和函数y作为输入变量，一个y的导函数作为输出变量。其中自变量t不论在ode文件中是否使用都必须作为第一输入变量，y则必须作为第二输入变量，位置不能颠倒。2.可以向ode文件中传递参数，数目不受限制第8页，本讲稿共38页odefile的编写function f=fun(x,y)f=y-2*x/y;求解初值问题：自变量在前，因变量在后ode输入函数输出变量为因变量导数的表达式初值问题：function f=fun(x,y)f=yy2;第9页，本讲稿共38页常微分方程组odefile的编写常微分方程组与单个常微

5、分方程求解方法相同，只需在编写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;第10页，本讲稿共38页高阶微分方程高阶微分方程odefile的编写的编写本例的难度：求解：y(0)=0，y(0)=1，方程系数非线性可在odefile中定义方程高阶，非标准形式方程变形：令y1y；y2y则原方程等价于：function f=fun(t,y)a=-exp(-t)+cos(2*p

6、i*t)*exp(-2*t);b=cos(2*pi*t);f=y(2)-a*y(2)2-b*y(1)+exp(t)*b;第11页，本讲稿共38页解算指令的使用方法解算指令的使用方法调用格式：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,)说明：v输出变量T为返回时间列向量；解矩阵Y的每一行对应于T的一个元素，列数与求解变量数相等。vfun为函数句柄，为

7、根据待求解的ODE方程所编写的ode文件（odefile）；vTSPANT0 TFINAL是微分系统yF(t,y)的积分区间；Y0为初始条件voptions用于设置一些可选的参数值，缺省时，相对于第一种调用格式。options中可以设置的参数参见odesetvP1，P2，的作用是传递附加参数P1，P2，到ode文件。当options缺省时，应在相应位置保留，以便正确传递参数。第12页，本讲稿共38页常微分方程初值问题解算指令比较常微分方程初值问题解算指令比较解算指令算法精度ode45四五阶Runge-Kutta法较高ode23二三阶Runge-Kutta法低ode113可变阶Adams-Bas

8、hforth-Moulton法ode15s基于数值差分的可变阶方法（BDFs，Gear）低中ode23s二阶改进的Rosenbrock法低ode23t使用梯形规则适中ode23tbTR-BDF2（隐式Runge-Kutta法）低第13页，本讲稿共38页ode解算指令的选择(1)求解初值问题：比较ode45和ode23的求解精度和速度 1.根据常微分方程要求的求解精度与速度要求 第14页，本讲稿共38页ode45和ode23的比较function Cha5demo1%Comparison of results obtained from ode45 and ode23 solverformat

9、longy0=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(Comparative Results at x=1:);fprintf(nODE45ttt y=%.8fnODE23ttt y=%.8fnPrecisive Result.=%.8fn,y1(end),y2(end),1.7320508)%-fun

10、ction f=fun(x,y)f=y-2*x/y;第15页，本讲稿共38页ode解算指令的选择(2)2.根据常微分方程组是否为刚性方程 如果系数矩阵A的特征值连乘积小于零，且绝对值最大和最小的特征值之比（刚性比）很大，则称此类方程为刚性方程 刚性方程在化学工程和自动控制领域的模型中比较常见。刚性比：100/0.0110000 第16页，本讲稿共38页ode解算指令的选择(2)刚性方程的物理意义：常微分方程组所描述的物理化学变化过程中包含了多个子过程，其变化速度相差非常大的数量级，于是常微分方程组含有快变和慢变分量。常微分方程组数值积分的稳定步长受模值最大的特征值控制，即受快变量分量约束，特征

11、值大则允许步长小；而过程趋于稳定的时间又由模值最小的特征值控制，特征值小则积分到稳定的时间则长。Matlab提供了不同种类的刚性方程求解指令：ode15s ode23s ode23t ode23tb，可根据实际情况选用第17页，本讲稿共38页function Cha5demo3figureode23s(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=d

12、y1dx;dy2dx;刚性常微分方程组求解第18页，本讲稿共38页解算指令的输出控制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(out

13、putfcn,odephas2)可以将输出方法改变为平面输出6.odeprint输出求解过程每一步的解第19页，本讲稿共38页解算指令的options选项1.RelTol相对误差，它应用于解向量的所有分量。在每一步积分过程中，第i个分量误差e(i)满足：e(i)=max(RelTol*abs(y(i),AbsTol(i)。2.AbsTol绝对误差，若是实数，则应用于解向量的所有分量，若是向量，则它的每一个元素应用于对应位置解向量元素。3.OutputFcn可调用的输出函数名。每一步计算完后，这个函数将被调用输出结果，可以选择的值为：odeplot，odephas2，odephas3，odepr

14、int。4.OutputSel输出序列选择。指定解向量的哪个分量被传递给OutputFcn。5.MaxSetp步长上界，缺省值为求解区间的1/10。6.InitialStep初始步长，缺省时自动设置。7.采用odeset改变原有选项的值第20页，本讲稿共38页在间歇反应器中进行液相反应制备产物B，反应网络如图5-1所示。反应可在180260的温度范围内进行，反应物X大量过剩，而C,D和E为副产物。各反应均为一级动力学关系：rkC，式中 间歇反应器中串联平行复杂反应系统已知参数：k015.780521010，k023.923171012，k031.64254104，k046.264108，Ea1

15、=124670，Ea2=150386，Ea3=77954，Ea4=111528。初始浓度：CA=1kmol/m3，其余物质浓度为0。已知是产物B收率最大的最优反应温度为224.6试计算1)在最优反应温度下各组分浓度随时间的动态变化；2)最优反应时间；3)输出产物D对反应物浓度A的关系图。第21页，本讲稿共38页间歇反应器中串联平行复杂反应系统数学模型 第22页，本讲稿共38页间歇反应器中串联平行复杂反应系统function Cha5demo4T=224.6+273.15;R=8.31434;k0=5.78052E+10 3.92317E+12 1.64254E+4 6.264E+8;Ea=12

16、4670 150386 77954 111528;%Initial concentration:C0(i),kmol/m3C0=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)%plotplot(t,C(:,1),r-,t,C(:,2),k:,t,C(:,3),b-.,t,C(:,4),k-);xlabel(Time(s);ylabel(Concentration(kmol/m3);legend

17、(A,B,C,D)CBmax=max(C(:,2);%CBmaxyBmax=CBmax/C0(1)%yBmax:index=find(C(:,2)=CBmax);t_opt=t(index)%t_opt:the optimum batch time,s%-function dCdt=MassEquations(t,C,k0,Ea,R,T)%Reaction rate constants,1/sk=k0.*exp(-Ea/(R*T);k(5)=2.16667E-04;%Reaction rates,kmoles/m3 srA=-(k(1)+k(2)*C(1);rB=k(1)*C(1)-k(3)*

18、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);%Mass balancesdCdt=rA;rB;rC;rD;rE;第23页，本讲稿共38页Matlab求解边值问题方法：bvp4c函数1.把待解的问题转化为标准边值问题 2.因为边值问题可以多解，所以需要为期望解指定一个初始猜测解。该猜测解网（Mesh）包括区间a,b内的一组网点（Mesh points）和网点上的解S(x)3.根据原微分方程构造残差函数4.利用原微分方程和边界条件，借助迭代不断产生新的S(x)，使残差不断减小，从而获得满足精度要

19、求的解 第24页，本讲稿共38页Matlab求解边值问题求解边值问题的基本指令的基本指令solinitbvpinit（x,v,parameters）生成bvp4c调用指令所必须的“解猜测网”solbvp4c（odefun，bcfun，solinit，options，p1，p2，）给出微分方程边值问题的近似解 sxintdeval（sol，xint）计算微分方程积分区间内任何一点的解值第25页，本讲稿共38页初始解生成函数：初始解生成函数：bvpinit()solinitbvpinit（x,v,parameters）x指定边界区间a,b上的初始网络，通常是等距排列的（1M）一维数组。注意：使x(

20、1)=a，x(end)=b；格点要单调排列。v是对解的初始猜测solinit（可以取别的任意名）是“解猜测网（Mesh）”。它是一个结构体，带如下两个域：solinit.x是表示初始网格有序节点的（1M）一维数组，并且solinit.x(1)一定是a，solinit.x(end)一定是b。M不宜取得太大，10数量级左右即可。solinit.y是表示网点上微分方程解的猜测值的（NM）二维数组。solinit.y(:,i)表示节点solinit.x(i)处的解的猜测值。第26页，本讲稿共38页solbvp4c（odefun，bcfun，solinit，options，p1，p2，）输入参数：ode

21、fun是计算导数的M函数文件。该函数的基本形式为：dydxodefun(x,y,parameters,p1,p2，)，在此，自变量x是标量，y，dydx是列向量。其基本形式与初值问题的odefile形式相同bcfun是计算边界条件下残数的M函数文件。其基本形式为：resbcfun(ya,yb,parameters,p1,p2,)，文件输入宗量ya,yb是边界条件列向量,分别代表y在a和b处的值。res是边界条件满足时的残数列向量。注意：例如odefun函数的输入宗量中包含若干“未知”和“已知”参数，那么不管在边界条件计算中是否用到，它们都应作为bcfun的输入宗量。输入宗量options是用来

22、改变bvp4c算法的控制参数的。在最基本用法中，它可以缺省，此时一般可以获得比较满意的边值问题解。如需更改可采用bvpset函数，使用方法同odeset函数。输入宗量p1，p2等表示希望向被解微分方程传递的已知参数。如果无须向微分方程传递参数，它们可以缺省。边值问题求解指令：bvp4c()第27页，本讲稿共38页输出参数：输出变量sol是一个结构体usol.x是指令bvp4c所采用的网格节点；usol.y是y(x)在sol.x网点上的近似解值；usol.yp是y(x)在sol.x网点上的近似解值；usol.parameters是微分方程所包含的未知参数的近似解值。当被解微分方程包含未知参数时，

23、该域存在。边值问题求解指令：bvp4c()第28页，本讲稿共38页边值问题的求解边值问题的求解原方程组等价于以下标准形式的方程组：function Cha5demo5solinit=bvpinit(linspace(0,1,10),1 0);sol=bvp4c(ODEfun,BCfun,solinit);x=0:0.05:0.5;y=deval(sol,x);yP=0-0.41286057-0.729740656.-0.95385538-1.08743325-1.13181116;xP=0:0.1:0.5;plot(xP,yP,o,x,y(1,:),r-),legend(Analytical

24、Solution,Numerical Solution)%-function dydx=ODEfun(x,y)dydx=y(2);y(1)+10;%-function bc=BCfun(ya,yb)bc=ya(1);yb(1);求解两点边值问题：令：边界条件为：第29页，本讲稿共38页边值问题的求解边值问题的求解原方程组等价于以下标准形式的方程组：function Cha5demo6solinit=bvpinit(linspace(0,1,10),0 1);sol=bvp4c(ODEfun,BCfun,solinit);x=0:0.1:1;y=deval(sol,x);yP=1 1.0743

25、1.1695 1.2869 1.4284.1.5965 1.7947 2.0274 2.3004 2.6214 3;xP=0:0.1:1.0;plot(xP,yP,o,x,y(1,:),r-)legend(Analytical Solution,Numerical Solution,location,Northwest)legend boxoff%-function dydx=ODEfun(x,y)dydx=y(2);(1+x2)*y(1)+1;%-function bc=BCfun(ya,yb)bc=ya(1)-1;yb(1)-3;求解：令：边界条件为：第30页，本讲稿共38页边值问题的求解

26、边值问题的求解function Cha5demo7c=1;solinit=bvpinit(linspace(0,4,10),1;1);sol=bvp4c(ODEfun,BCfun,solinit,c);x=0:0.1:4;y=deval(sol,x);plot(x,y(1,:),b-,sol.x,sol.y(1,:),ro)legend(解曲线,初始网格点解)%-function dydx=ODEfun(x,y,c)dydx=y(2);-c*abs(y(1);%-function bc=BCfun(ya,yb)bc=ya(1);yb(1)+2;求解：令：c1 此程序有什么错误？第31页，本讲稿

27、共38页边值问题的求解边值问题的求解function Cha5demo8lmb=15;solinit=bvpinit(linspace(0,pi,10),1;1,lmb);opts=bvpset(Stats,on);sol=bvp4c(ODEfun,BCfun,solinit,opts);lambda=sol.parametersx=0:pi/60:pi;y=deval(sol,x);plot(x,y(1,:),b-,sol.x,sol.y(1,:),ro)legend(解曲线,初始网格点解)%-function dydx=ODEfun(x,y,lmb)q=5;dydx=y(2);-(lmb-

28、2*q*cos(2*x)*y(1);%-function bc=BCfun(ya,yb,lmb)bc=ya(1)-1;ya(2);yb(2);求解：边界条件：本例中，微分方程与参数的数值有关。一般而言，对于任意的值，该问题无解，但对于特殊的值（特征值），它存在一个解，这也称为微分方程的特征值问题。对于此问题，可在bvpinit中提供参数的猜测值，然后重复求解BVP得到所需的参数，返回参数为sol.parameters 第32页，本讲稿共38页边值问题的求解边值问题的求解function Cha5demo9solinit=bvpinit(linspace(0,1,5),ex1init);opti

29、ons=bvpset(Stats,on,RelTol,1e-5);sol=bvp4c(ODEfun,BCfun,solinit,options);x=sol.x;y=sol.y;.%-function dydx=ODEfun(x,y)dydx=0.5*y(1)*(y(3)-y(1)/y(2)-0.5*(y(3)-y(1)(0.9-1000*(y(3)-y(5)-0.5*y(3)*(y(3)-y(1)/y(4)0.5*(y(3)-y(1)100*(y(3)-y(5);%-function bc=BCfun(ya,yb)bc=ya(1)-1 ya(2)-1 ya(3)-1 ya(4)+10 yb(

30、3)-yb(5);%-function v=ex1init(x)v=1；1；-4.5*x2+8.91*x+1；-10；-4.5*x2+9*x+0.91;求解：边界条件：初始猜测解为：第33页，本讲稿共38页边值问题的求解边值问题的求解function Cha5demo11infinity=6;solinit=bvpinit(linspace(0,infinity,5),0 0 1);options=bvpset(stats,on);sol=bvp4c(ex5ode,ex5bc,solinit,options);eta=sol.x;f=sol.y;fprintf(n);fprintf(Cebec

31、i&Keller report f(0)=0.92768.n)fprintf(Value computed here is f(0)=%7.5f.n,f(3,1)clf resetplot(eta,f(2,:);axis(0 infinity 0 1.4);title(Falkner-Skan equation,positive wall shear,beta=0.5.)xlabel(eta),ylabel(df/deta),shg%-function dfdeta=ex5ode(eta,f)beta=0.5;dfdeta=f(2)；f(3)；-f(1)*f(3)-beta*(1-f(2)2)

32、;%-function res=ex5bc(f0,finf)res=f0(1)；f0(2)；finf(2)-1;求解：边界条件：0.5 如果取1，计算结果如何第34页，本讲稿共38页常微分方程边值问题的应用已知在球形氧化铝催化剂进行环己烷脱氢反应，催化剂直径为5mm，操作温度为700K，在此温度下，k4s-1，D=0.05cm2/s。计算催化剂颗粒内环己烷的浓度分布及催化剂的有效因子。数学模型：质量传递方程 边界条件：0r1 等温有效因子为：Thiele模数 第35页，本讲稿共38页function Cha5demo12 global fai CsCs=1;rp=5e-3/2;k=4;D=0.

33、05e-4;fai=rp*sqrt(k/D)%fai:Thiele模数a=0;b=1;solinit=bvpinit(linspace(a,b,100),0 0);sol=bvp4c(ODEs,BCfun,solinit);plot(sol.x,sol.y(1,:),b-)xlabel(Dimensionless Radius)，ylabel(Dimensionless Concentration)I1=quadl(IntFunc1,a,b,sol.x,sol.y(1,:)I2=quadl(IntFunc2,a,b)eta=I1/I2fprintf(t催化剂的有效因子为:=%.4fn,eta)

34、%-function dydr=ODEs(r,y)global fai Csdy1dr=y(2);if r=0 dy2dr=0;else dy2dr=fai2*y(1)/Cs-2/r*y(2);enddydr=dy1dr;dy2dr;%-function bc=BCfun(ya,yb)bc=yb(1)-1;ya(2);%-function f=IntFunc1(r,ri,yi)f=spline(ri,yi,r).*r.2;%-function f=IntFunc2(r)global Csf=Cs*r.2;第36页，本讲稿共38页化工数学模型中的微分方程化工数学模型中的微分方程1.集中参数模型和分布参数模型集中参数模型是指因变量不随空间坐标变化，模型中自变量仅为时间；分布参数模型则指因变量不仅与时间有关，而且与空间有关。集中参数模型为常微分方程，分布参数模型为偏微分方程。2.化工模型中最常见的常微分方程边值问题模型为涉及传热和扩散的问题3.在建立微分方程模型时，请注意初始和边界条件，这样才是一个完整的定解问题4.化工模型中的微分方程种类很多，但其求解过程往往并不复杂。第37页，本讲稿共38页明天休息1天！自行3人一组进行组队，选择组长，组长上报组员名单！第38页，本讲稿共38页