Differential model won't calculate with powers <1
2 views (last 30 days)
Show older comments
Jort Puiman
on 8 Mar 2021
Commented: Jort Puiman
on 10 Mar 2021
Hi everyone,
I am working on a model which describes multiple reactions over time, using simple power law models (R = k*[C]^n).
ode1 = diff(Ca) == -k1*Ca;
cond1 = Ca(0) == 5;
ode2 = diff(Cb) == k1*Ca - k2*Cb;
cond2 = Cb(0) == 0;
ode3 = diff(Cc) == k2*Cb - k3*Cc;
cond3 = Cc(0) == 0;
ode4 = diff(Cd) == k3*Cc - k4*Cd;
cond4 = Cd(0) == 0;
ode5 = diff(Ce) == k4*Cd;
cond5 = Ce(0) == 0;
t works just like I want to with n = 1, however, our data suggests that n < 1. I tried adding powers to my concentrations, but then, Matlab has a hard time calculating it, and it never finishes.
I want to calculate the concentrations of all components over time. All constants (k1, k2, k3, k4) and reaction orders will be added later. However, the model should be able to calculate the concentrations when the order is lower than 1.
Do you have any suggestions on how to tackle this problem? Do I have to use other functions?
2 Comments
Jan
on 8 Mar 2021
You forgot to mention how you try to calculate what. These are the equations only.
Accepted Answer
Bjorn Gustavsson
on 8 Mar 2021
It is not given that a general nonlinear system of ODEs have analytical solutions. If you cannot get the symbolic toolbox to find one for you, the best (possibly?) advice is to turn to numerical solutions. For that have a look at ode45 and its siblings, and the large number of ODE-examples. This should be very straightforward to implement for numerical integration. The ODE-function might look something like this:
function dCdt = your_power_reactions_ODE(t,C,k,g)
Ca = C(1);
Cb = C(2);
Cc = C(3);
Cd = C(4);
Ce = C(5);
dCadt = -k(1)*Ca^g(1);
dCbdt = k(1)*Ca^g(1) - k(2)*Cb^g(2);
dCcdt = k(2)*Cb^g(2) - k(3)*Cc^g(3);
dCddt = k(3)*Cc^g(3) - k(4)*Cd^g(4);
dCedt = k(4)*Cd^g(4);
% etc...
end
That should be pretty standard to integrate. It might be a numerical stiff set of ODEs - I've not spent any time to judge that. You then simply set the initial conditions and the time-period and then run:
C0 = [5 0 0 0 0];
t_span = [0 17];
k_all = [1 log(2),2^pi,12];
g_all = 1./[1:4];
[t_out,C_out] = ode45(@(t,C) your_power_reactions_ODE(t,C,k_all,g_all),...
t_span,C0);
If it is too stiff, try ode15s, ode113, or ode23s...
HTH
3 Comments
Bjorn Gustavsson
on 8 Mar 2021
But to my inderstanding these are all "decay-type reactions" - like for the concentrations of elements in a radioactive decay-chain?
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!