Differential model won't calculate with powers <1

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?

[Jan]

You forgot to mention how you try to calculate what. These are the equations only.

[Jort Puiman]

Thank you for your quick remark, I edited my question. Hope it is clear right now!

Answer 1

[Bjorn Gustavsson]

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

[Jort Puiman]

Great answer! Thank you! I will take a look at it if it works and come back to you tomorrow. This will help me greatly!

[Bjorn Gustavsson]

But to my inderstanding these are all "decay-type reactions" - like for the concentrations of elements in a radioactive decay-chain?

[Jort Puiman]

Yes you could compare it to a radioactive decay-chain. We start with a set concentration of Ca, and it oxidises to Cb, then Cc, etc... I however want to get as much of Cd, and as few of Ce. That's the general idea.