《常微分方程(组)的数值解法》由会员分享，可在线阅读，更多相关《常微分方程(组)的数值解法（20页珍藏版）》请在金锄头文库上搜索。
1、2.3 常微分方程（组）的数值解法,知识要点,常微分方程初值问题-ode45,0de23,微分方程在化工模型中的应用,间歇反应器的计算 活塞流反应器的计算 全混流反应器的动态模拟 定态一维热传导问题 逆流壁冷式固定床反应器一维模型 固定床反应器的分散模型,Matlab常微分方程求解问题分类,初值问题： 定解附加条件在自变量的一端 一般形式为： 初值问题的数值解法一般采用步进法，如Runge-Kutta法,边值问题： 在自变量两端均给定附加条件 一般形式： 边值问题可能有解、也可能无解，可能有唯一解、也可能有无数解 边值问题有3种基本解法 迭加法 打靶法 松弛法,Matlab求解常微分方程初值问
2、题方法,将待求解微分方程（组）转化为标准形式， “翻译”成Matlab可以理解的语言 编写odefile文件,选择合适的解算指令求解问题,根据求解问题的要求，设置解算指令的调用格式,Matlab求解初值问题函数,odefile,odefile是一个Matlab函数文件，一般作为整个求解程序的一个子函数，表示ode求解问题 对于程序通用性要求不高的场合，只需将原有模型写成标准形式，然后“翻译”成Matlab语言即可,odefile的编写-常微分方程,function f=fun(x,y) f=y-2*x/y;,求解初值问题：,自变量在前，因变量在后,ode输入函数,输出变量为因变量导数的表达式,
3、初值问题：,function f=fun(x,y) f=yy2;,odefile的编写-常微分方程组,常微分方程组与单个常微分方程求解方法相同，只需在编写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;,高阶微分方程odefile的编写,本例的难度：,求解：,y(0)=0，y(0)=1，,方程系数非线性,可在odefile中定义,方程高阶，非标
4、准形式,方程变形：令y1y；y2y 则原方程等价于：,function f=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）； TSPANT0 TFINAL是微分系统yF(t,y)的积分区间；Y0为初始条件,常微分方程初值问题解算
5、指令比较,ode解算指令的选择(1),求解初值问题：,比较ode45和ode23的求解精度和速度,1.根据常微分方程要求的求解精度与速度要求,ode45和ode23的比较-1,function xODE clear all clc format long y0 = 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) % - function dydx = f(x,y) dydx = y - 2*x/y;,ode45和ode23的比较-2,function
6、 xODEs clear all clc y0 = 0 1; x1,y1 = ode45(f,0:0.2:1,y0); x2,y2 = ode23(f,0:0.2:1,y0); disp(Results by using ode45()：) disp( x y(1) y(2) disp(x1 y1); disp(Results by using ode23()：) disp( x y(1) y(2) disp(x2 y2); % - function dydx = f(x,y) % Define simultaneous ODE equations dydx = 3*y(1) + 2*y(2)
7、; 4*y(1)+y(2);,ode解算指令的选择(2),2.根据常微分方程组是否为刚性方程,如果系数矩阵A的特征值连乘积小于零，且绝对值最大和最小的特征值之比（刚性比）很大，则称此类方程为刚性方程,刚性方程在化学工程和自动控制领域的模型中比较常见。,刚性比：100/0.0110000,ode解算指令的选择(2),刚性方程的物理意义：常微分方程组所描述的物理化学变化过程中包含了多个子过程，其变化速度相差非常大的数量级，于是常微分方程组含有快变和慢变分量。 Matlab提供了不同种类的刚性方程求解指令：ode15s ode23s ode23t ode23tb，可根据实际情况选用,function
8、 demo1 figure ode23s(fun,0,100,0;1) hold on, pause 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 = dy1dx; dy2dx;,刚性常微分方程组求解,间歇反应器中串联平行复杂反应系统（例4-1）,function BatchReactor clear all clc T = 224.6 + 273.15; % React
9、or temperature, Kelvin R = 8.31434; % Gas constant, kJ/kmol K % Arrhenius constant, 1/s k0 = 5.78052E+10 3.92317E+12 1.64254E+4 6.264E+8; % Activation energy, kJ/kmol Ea = 124670 150386 77954 111528; % 初始浓度C0(i), kmol/m3 C0 = 1 0 0 0 0; tspan = 0 1e4; t,C = ode45(MassEquations, tspan, C0,k0,Ea,R,T);
10、 % 绘图 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: the maximum concentration of B, kmol/m3 yBmax = CBmax/C0(1) % yBmax: the maximum yield of B index = find(C(:,2)=CBmax); t_opt = t(index) % t_opt
11、: the optimum batch time, s % - function dCdt = MassEquations(t,C,k0,Ea,R,T) % Reaction rate constants, 1/s k = k0.*exp(-Ea/(R*T); k(5) = 2.16667E-04; % Reaction rates, kmoles/m3 s rA = -(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
12、)*C(4); % Mass balances dCdt = rA; rB; rC; rD; rE;,三个串联的CSTR等温反应器（例4-3）,function IsothermCSTRs clear all clc CA0 = 1.8; % kmol/m3 CA10 = 0.4; % kmol/m3 CA20 = 0.2; % kmol/m3 CA30 = 0.1; % kmol/m3 k = 0.5; % 1/min tau = 2; stoptime = 2.9; % min t,y = ode45(Equations,0 stoptime,CA10 CA20 CA30,k,CA0,ta
13、u); disp( Results:) disp( t CA1 CA2 CA3) 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) % - function dydt = 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;,
