Solving ODEs for chemical rate equations

I am trying to model a chemical equation of the form rate = k [a]^1 *[b]^-1. I have managed witht the help of some you tube resources to model a an equation for one reactant but am unsure of how to include the second. Here is the code i have got so far and any help would be appreciated.

clear all 
close all 
clc

K = 0.06;
timespan = [0 30]';
A0 = 0.5;

first = @(t,A) -K*A
[t,A_calc] = ode45(first,timespan,A0)
plot(t,A_calc, 'o','Linewidth',1,'Markersize',10)

[Torsten]

Isn't your equation

dC_A/dt = -k*A

the equation for

rate = k*[A]

?

[Oliver Wilson]

Yes sorry if that was unclear. I am aiming for rate = k [a]^1 *[b]^-1 and am not quite sure how to add the second concentration

Answer 1

[Star Strider]

Here is prototype code for integrating the solution for a different kinetic parameter estimation problem with four rate equations —

t=[0.1
    0.2
    0.4
    0.6
    0.8
    1
    1.5
    2
    3
    4
    5
    6];

c=[0.902	0.06997	0.02463	0.00218
    0.8072	0.1353	0.0482	0.008192
    0.6757	0.2123	0.0864	0.0289
    0.5569	0.2789	0.1063	0.06233
    0.4297	0.3292	0.1476	0.09756
    0.3774	0.3457	0.1485	0.1255
    0.2149	0.3486	0.1821	0.2526
    0.141	0.3254	0.194	0.3401
    0.04921	0.2445	0.1742	0.5277
    0.0178	0.1728	0.1732	0.6323
    0.006431	0.1091	0.1137	0.7702
    0.002595	0.08301	0.08224	0.835]; 

theta0 = rand(10,1);

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

for k1 = 1:numel(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end

		Theta( 1) =  0.76483
		Theta( 2) =  0.23491
		Theta( 3) =  0.20880
		Theta( 4) =  0.49182
		Theta( 5) =  0.62211
		Theta( 6) =  0.00000
		Theta( 7) =  0.90287
		Theta( 8) =  0.07145
		Theta( 9) =  0.02840
		Theta(10) =  0.00000

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')

function C=kinetics(theta,t)
    c0=theta(end-3:end);
    [T,Cv]=ode45(@DifEq,t,c0);
    %
    function dC=DifEq(t,c)
        dcdt=zeros(4,1);
        dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
        dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
        dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
        dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
        dC=dcdt;
    end
    C=Cv;
end

Each rate equation needs to be coded separately (here in the ‘DifEq’ function), then the integration can proceed.

.

[Oliver Wilson]

Thank you very much.

Just so i can check my understanding of your code is dcdt (1) equivalent to the rate equation of the following form rate = K[a]^-1[b]

[Star Strider]

As always, my pleasure!

I would code that as:

dcdt(1) = K(1)*C(2)/C(1);

where ‘K’ is a vector of parameters to be estimated (‘theta’ in the code I posted), ‘C(1)’ is the concentration of ‘a’ and ‘C(2)’ is the concentration of ‘b’. Using the -1 exponent to indicate division is inefficient. Just do the division directly.

Code the other equations similarly.

I keep track of the vector elements that represent concentrations in a comment line, for example:

% C(1) = a, C(2) = b, C(3) = c

and if necessary, do the same thing for the parameters:

% K(1) = K_a, K(2) = K_b

in order to not lose track of them. That also makes it easier to display them correctly later.

.