87 views (last 30 days)

Hello. I am new to matlab and I have a set of kinetic equations that need to be solve. The reaction has 1 reactant in which it will decompose into 3 other product over time. I have the experimental results that shows the degradation of reactant and the concentration of products over a period of time, as below. The differential equations I need to solve are also as below.

Reaction Kinetics:

d[A]/dt = -k1[A]-k2[A]

d[B]/dt = k1[A]-k3[B]-k4[B]

d[C]/dt = k2[A]+k3[B]

d[D]/dt = k4[B]

Experimental result

Time (min) / Conc (g/L) [A] [B] [C] [D]

0 1.000 0 0 0

20 0.8998 0.0539 0.0039 0.0338

30 0.8566 0.0563 0.00499 0.0496

60 0.7797 0.0812 0.00715 0.0968

90 0.7068 0.0918 0.00937 0.1412

120 0.6397 0.0989 0.01028 0.1867

I have browse through various code and the one solved StarStrider for Igor Moura shows the result that I wish to achieve where plotted is the graph of the concentrations of reactant and products over time as well as solving the reaction rate constants "ki". However I have modified the code according to my equations and the graph came out weird and the results are not fitted well into with the experimental results. Can someone help me with this? Attach below is the code I have tried out.

function Igor_Moura

% 2016 12 03

% NOTES:

%

% 1. The ‘theta’ (parameter) argument has to be first in your

% ‘kinetics’ funciton,

% 2. You need to return ALL the values from ‘DifEq’ since you are fitting

% all the values

function C=kinetics(k,t)

c0=[1;0;0;0];

[T,Cv]=ode45(@DifEq,t,c0);

%

function dC=DifEq(t,c)

dcdt=zeros(4,1);

dcdt(1)=-k(1).*c(1)-k(2).*c(1);

dcdt(2)= k(1).*c(1)-k(3).*c(2)-k(4).*c(2);

dcdt(3)= k(2).*c(1)+k(3).*c(2);

dcdt(4)= k(4).*c(2);

dC=dcdt;

end

C=Cv;

end

t=[20

30

60

90

120];

c=[0.91257 0.02086 0.00939 0.00725

0.88232 0.02531 0.01664 0.00897

0.83324 0.02882 0.03927 0.01195

0.76289 0.03137 0.06834 0.01514

0.70234 0.03380 0.10542 0.01679 ];

k0=[1;1;1;1];

[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,t,c);

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

for k1 = 1:length(k)

fprintf(1, '\t\k(%d) = %8.5f\n', k1, k(k1))

end

tv = linspace(min(t), max(t));

Cfit = kinetics(k, tv);

figure(1)

plot(t, c, 'p')

hold on

hlp = plot(tv, Cfit);

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')

end

Star Strider
on 22 Sep 2020 at 14:23

Edited: Star Strider
on 22 Sep 2020 at 14:37

Using this version of ‘kinetics’ (that also estimates the initial conditions for ‘DifEq’),

function C=kinetics(k,t)

c0=k(5:8);

[T,Cv]=ode45(@DifEq,t,c0);

%

function dC=DifEq(t,c)

dcdt=zeros(4,1);

dcdt(1)=-k(1).*c(1)-k(2).*c(1);

dcdt(2)= k(1).*c(1)-k(3).*c(2)-k(4).*c(2);

dcdt(3)= k(2).*c(1)+k(3).*c(2);

dcdt(4)= k(4).*c(2);

dC=dcdt;

end

C=Cv;

end

I got a good fit with these parameters (note that the first 4 are ‘k(1)’ through ‘k4’ and the last 4 are the initial conditions):

Rate Constants:

Theta(1) = 0.00180

Theta(2) = 0.00049

Theta(3) = 0.02622

Theta(4) = 0.01094

Theta(5) = 0.89995

Theta(6) = 0.00011

Theta(7) = 0.00213

Theta(8) = 0.00436

producing this plot:

I used the ga (genetic algorithm) function to find them. I am cleaning up the code, and will post the entire code file in a few minutes.

EDIT — (22 Sep 2020 at 14:36)

The complete (cleaned) code:

function Daphne_Deidre_Wong_GA

% 2020 09 22

% NOTES:

%

% 1. The ‘k’ (parameter) argument has to be first in your

% ‘kinetics’ function,

% 2. You need to return ALL the values from ‘DifEq’ since you are fitting

% all the values

function C=kinetics(k,t)

c0=k(5:8);

[T,Cv]=ode45(@DifEq,t,c0);

%

function dC=DifEq(t,c)

dcdt=zeros(4,1);

dcdt(1)=-k(1).*c(1)-k(2).*c(1);

dcdt(2)= k(1).*c(1)-k(3).*c(2)-k(4).*c(2);

dcdt(3)= k(2).*c(1)+k(3).*c(2);

dcdt(4)= k(4).*c(2);

dC=dcdt;

end

C=Cv;

end

format long

t=[ 20

30

60

90

120];

c=[0.91257 0.02086 0.00939 0.00725

0.88232 0.02531 0.01664 0.00897

0.83324 0.02882 0.03927 0.01195

0.76289 0.03137 0.06834 0.01514

0.70234 0.03380 0.10542 0.01679];

ftns = @(theta) norm(c-kinetics(theta,t));

PopSz = 500;

Parms = 8;

opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-5, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);

t0 = clock;

fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

[theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)

t1 = clock;

fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

GA_Time = etime(t1,t0)

DT_GA_Time = datetime(num2str(GA_Time,'%.6f'),'InputFormat','ss.SSSSSS', 'Format','HH:mm:ss.SSS')

fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

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

for k1 = 1:length(theta)

fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))

end

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, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','E')

format short eng

end

.

Star Strider
on 5 Oct 2020 at 11:00

If you want to post your code, I will of course look at it and do my best to answer whatever questions you have about it.

I do not understand about ‘time constraint’ however. Please explain that in a bit more detail.

Star Strider
on 6 Oct 2020 at 11:06

Change ‘DT_GA_Time’ to this:

DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');

and the problem ino longer exists.

Opportunities for recent engineering grads.

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

Start Hunting!
## 7 Comments

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1019827

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1019827

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1019875

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1019875

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1020511

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1020511

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1020850

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1020850

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1022185

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1022185

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1022839

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1022839

## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1022866

⋮## Direct link to this comment

https://de.mathworks.com/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#comment_1022866

Sign in to comment.