使用 ode45 在 MATLAB 中求解 4 ODE 系统

Question

我不太习惯 MATLAB，我正在尝试使用 MATLAB ode45 解决以下问题，但是它不起作用。

我正在使用半间歇式反应器解决反应工程中的问题。 反应由下式给出 A + B ---> C + D 将 A 放入反应器中，将 B 以 v0 = 0.05 L/s 的流量连续添加到反应器中。初始体积为 V0 = 5 L。该反应是基本反应。反应常数为 k = 2.2 L/mol.s。 初始浓度：A：0.05 M，B：0.025 M。

对反应器中各物种进行摩尔平衡，得到以下4个ODE，以及V的表达式（反应器的体积不断增加）

解决这个系统并根据时间绘制解决方案，我应该得到这个

请注意，C(C) 和 C(D) 的图是相同的。 让我们设置 tau = v0/V。

现在是 MATLAB 代码部分。

我在网上进行了广泛的搜索，根据我所学到的，我想出了以下代码。

首先，我编写了ODE系统的代码

function f = ODEsystem(t, y, tau, ra, y0)
f = zeros(4, 1);
f(1) = ra - tau*y(1);
f(2) = ra + tau*(y0(2) - y(2));
f(3) = -ra - tau*y(3);
f(4) = -ra - tau*y(4);
end

然后，在命令窗口中，

t = [0:0.01:5];
v0 = 0.05;
V0 = 5;
k = 2.2;
V = V0 + v0*t;
tau = v0./V;
syms y(t);
ra = -k*y(1)*y(2);
y0 = [0.05 0.025 0 0];
[t, y] = ode45(@ODEsystem(t, y, tau, ra, y0), t, y0); 
plot(t, y); 

但是，我明白了......

如果有人可以帮助我修复我的代码，请。这真的很烦人:)

Answer 1

ra

不应作为参数传递，而应在 ODE 系统内部计算。

V

同样不是常数。符号表达式应用于公式转换，而不是数值方法。人们还必须明确地计算所需数值的符号表达式。

function f = ODEsystem(t, y, k, v0, V0, cB0)
f = zeros(4, 1);
ra = -k*y(1)*y(2);
tau = v0/(V0+t*v0);
f(1) = ra - tau*y(1);
f(2) = ra + tau*(cB0 - y(2));
f(3) = -ra - tau*y(3);
f(4) = -ra - tau*y(4);
end

然后使用图形的时间跨度，从除

A

之外的所有浓度为零开始，仅对流入使用浓度

B

。

t = [0:1:500];
v0 = 0.05;
V0 = 5;
k = 2.2;
cB0 = 0.025;
y0 = [0.05 0 0 0];
[t, y] = ode45(@(t,y) ODEsystem(t, y, k, v0, V0, cB0), t, y0); 
plot(t, y); 

并获得参考图像的良好再现

Answer 2

我们需要求解以下非齐次微分方程组，其中伴生矩阵A如下图所示：

可以使用以下matlab代码来解决：

%Define parameters 
v0 = 0.05;
V0 = 5;
k = 2.2;
y0 = [0.05 0.025 0 0];
t = [0:0.01:300];

%Numerically solve DE
[t,y] = ode45(@(t,y) diag([-v0/(V0 + v0*t),-v0/(V0 + v0*t),-v0/(V0 + v0*t),-v0/(V0 + v0*t)])*y + [-k*y(1)*y(2);-k*y(1)*y(2) + v0/(V0 + v0*t)*(y0(2) - y(2));--k*y(1)*y(2);--k*y(1)*y(2)],t,y0);

%Plot steady state solution
plot(t, y); 

使用 matlab 在线获得以下输出：