ODE45 code won't run

11 次查看（过去 30 天）
Najeem Afolabi
Najeem Afolabi 2020-12-25
评论： Najeem Afolabi 2020-12-27
So I have this complex reaction model that I made but it just won't run. It doesn't say that there is an error necessarily, it just tries to run for a really long time. I eventually stop the the attempt and it takes me to the ode45 function script and points at a line. This is the function and code that I'm using.
clear all; close all;
clc
%Initial Concentrations
CA0 = 0.5; %mol/L
CB0 = 0.1; %mol/L
Qoa = 0.000025; %L/s
Qob = 0.000025; %L/s
Q1 = Qob + Qoa;
CAii1 = Qoa*CA0/Q1; CB1 = Qob*CB0/Q1; CAiiB = 0; CAiB = 0; CAi = 0; CC1 = 0; CD1 = 0;
initCon = [CAii1, CB1, CAiiB, CAi, CAiB, CC1, CD1];
V = 10*10^-3;
k2 = 0.01;
%Residence Time
tau = V/Q1;
%Volume
tauRange = [0:0.2:tau]; %s
%ODE SOLVER
[tauRange, C1] = ode45(@diffpfr, tauRange, initCon, []);
%PLOT PFR
figure (1)
plot(tauRange, C1,'Linewidth',1)
function dcdv = diffpfr(tau, C)
a = 2; b = 0.25; c = 1; d = 0.25;
k1 = 1; k1_ = 1; k2_ = 0.25;k2 = 0.25;
%Rate Equations
rAii = - a * C(1)*C(2) + b * C(3);
rB = - a* C(1)* C(2) + b * C(3) - c * C(4) * C(2) + d * C(7);
rAiiB = -k1 * C(4) * C(5) - k2 * C(4) * C(6) + a * C(1)*C(2) - b * C(3);
rAi = k1 * C(4) * C(5) + k2 * C(4) * C(6) - c * C(4) * C(2) + d * C(7);
rAiB = c * C(4) * C(2) - d * C(7) - k1_ * C(7) - k2_ * C(7);
rC = k1 * C(3) + k1_ * C(7);
rD = k2 * C(3) + k2_ * C(7);
%differential equations
dcdv = [ rAii;rB; rAiiB; rAi; rAiB; rC; rD];
end
  2 个评论
Walter Roberson
Walter Roberson 2020-12-25
Experiment with using ode23s() as perhaps your ode is stiff.
Najeem Afolabi
Najeem Afolabi 2020-12-27
Yeah, I tried that too but still no luck unfortunately.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Cris LaPierre
Cris LaPierre 2020-12-27
编辑：Cris LaPierre 2020-12-27
Ran in:
There's something with your tau. Right now it is 200. When I make it <=30 I get it to run very quickly.
If instead I use ode15s, it works as is.
%Initial Concentrations
CA0 = 0.5; %mol/L
CB0 = 0.1; %mol/L
Qoa = 0.000025; %L/s
Qob = 0.000025; %L/s
Q1 = Qob + Qoa;
CAii1 = Qoa*CA0/Q1;
CB1 = Qob*CB0/Q1;
CAiiB = 0;
CAiB = 0;
CAi = 0;
CC1 = 0;
CD1 = 0;
initCon = [CAii1, CB1, CAiiB, CAi, CAiB, CC1, CD1];
V = 10*10^-3;
k2 = 0.01;
%Residence Time
tau = 200;%V/Q1;
%Volume
tauRange = [0 tau]; %s
%ODE SOLVER ### Switch to ode15s
[tauRange, C1] = ode15s(@diffpfr, tauRange, initCon);
%PLOT PFR
plot(tauRange, C1,'Linewidth',1)
function dcdv = diffpfr(tau, C)
a = 2; b = 0.25; c = 1; d = 0.25;
k1 = 1; k1_ = 1; k2_ = 0.25;k2 = 0.25;
%Rate Equations
rAii = - a*C(1)*C(2) + b*C(3);
rB = - a*C(1)*C(2) + b*C(3) - c*C(4)*C(2) + d*C(7);
rAiiB = -k1*C(4)*C(5) - k2*C(4)*C(6) + a*C(1)*C(2) - b*C(3);
rAi = k1*C(4)*C(5) + k2*C(4)*C(6) - c*C(4)*C(2) + d*C(7);
rAiB = c*C(4)*C(2) - d*C(7) - k1_*C(7) - k2_*C(7);
rC = k1*C(3) + k1_*C(7);
rD = k2*C(3) + k2_*C(7);
%differential equations
dcdv = [ rAii; rB; rAiiB; rAi; rAiB; rC; rD];
end
  2 个评论
Cris LaPierre
Cris LaPierre 2020-12-27
For more guidance on choosing a solver, see the Basic Solver Selection setion of the Choose and ODE Solver page.
Najeem Afolabi
Najeem Afolabi 2020-12-27
Ahh I see. Thank you very much. I'll apply this now.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by