常微分方程求解课件.ppt

常微分方程数值解 化工学院软件应用教科组 2006-10 本章知识要点 Matlab常微分方程求解问题分类 Matlab求解初值问题函数 odefile 所谓的odefile实际上是一个Matlab函数文件，一般作为整个求解程序的一个子函数，表示ode求解问题 Matlab提供了odefile的模板，采用type odefile命令显示其详细内容，然后将其复制到脚本编辑窗口，在合适的位置填入所需内容 一般而言，对于程序通用性要求不高的场合，只需将原有模型写成标准形式，然后“翻译”成Matlab语言即可 odefile的编写规定 常微分方程组odefile的编写 高阶微分方程odefile的编写 解算指令的输出控制 解算指令的options选项 边值问题求解指令：bvp4c () 边值问题求解指令：bvp4c () 常微分方程边值问题的应用 已知在球形氧化铝催化剂进行环己烷脱氢反应，催化剂直径为5mm，操作温度为700K，在此温度下，k＝4s-1，D=0.05cm2/s。计算催化剂颗粒内环己烷的浓度分布及催化剂的有效因子。 明天休息1天！ 自行3人一组进行组队，选择组长， 组长上报组员名单！ 边值问题的求解 function Cha5demo8 lmb=15; solinit = bvpinit(linspace(0,pi,10),[1;1],lmb); opts=bvpset(Stats,on); sol = bvp4c(@ODEfun,@BCfun,solinit,opts); lambda=sol.parameters x = [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-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 边值问题的求解 function Cha5demo9 solinit = bvpinit(linspace(0,1,5),@ex1init); options = 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(3) - yb(5)]; %------------------------------------------------------------------------- function v = ex1init(x) v = [1；1； -4.5*x^2+8.91*x+1