Solving ODEs for chemical rate equations

7 views (last 30 days)
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)
  2 Comments
Oliver Wilson
Oliver Wilson on 15 Nov 2022
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

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 15 Nov 2022
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.
.
  2 Comments
Star Strider
Star Strider on 15 Nov 2022
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.
.

Sign in to comment.

More Answers (0)

Categories

Find more on Chemistry in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!