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; endC=Cv;endt=[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))endtv = linspace(min(t), max(t));Cfit = kinetics(k, tv);figure(1)plot(t, c, 'p')hold onhlp = plot(tv, Cfit);hold offgridxlabel('Time')ylabel('Concentration')legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')end
Best Answer