Kinetics model-Why i dont get 3 graphs of these modeling, but i get just one for C catalyst. Is there problem with definition of k?

function Pd

t_span=[0 1000];

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

[t,c]=ode45(@kinetics,t_span,c0);

plot(t,c)

xlabel('Time,[min]');

ylabel('Concentration,[mol/cm^3]');

legend('cct-CDT','ctt-CDT','ttt-CDT','cc-CDD', 'ct-CDD', 'tt-CDD', 'c-CDE', 't-CDE', 'CDA')

end

function dz=kinetics(t,c);

k=[0.105 0.097 0.065 0.071 0.060 0.056 0.065 0.250 0.061 0.07 0.0110 0.120; %%C

0.015 0.007 0.010 0.008 0.025 0.016 0.005 0.003 0.007 0.015 0.001 0.035; %%CaCO3

0.090 0.078 0.070 0.058 0.125 0.045 0.069 0.007 0.073 0.080 0.006 0.039]; %%Al2O3

dz=zeros(9,3);

r1=k(1)*c(1);

r2=k(2)*c(1);

r3=k(3)*c(2);

r4=k(4)*c(2);

r5=k(5)*c(3);

r6=k(1)*c(1)-k(6)*c(4);

r7=k(2)*c(1)+k(3)*c(2)-k(7)*c(5)-k(8)*c(5)-k(9)*c(5);

r8=k(4)*c(2)+k(5)*c(3)-k(10)*c(6);

r9=k(6)*c(4)+k(7)*c(5)-k(11)*c(7);

r10=k(9)*c(5)+k(10)*c(6)-k(12)*c(8);

r11=k(11)*c(7)+k(8)*c(5)+k(12)*c(8);

dz(1)=-r1-r2;

dz(2)=-r3-r4;

dz(3)=-r5;

dz(4)=r6;

dz(5)=r7;

dz(6)=r8;

dz(7)=r9;

dz(8)=r10;

dz(9)=r11;

end

Answer 1

[Star Strider]

If you want to fit your differential equations to the ‘k’ matrix to estimate the parameters of your system, your current approach is not going to do that. See for example Parameter Estimation for a System of Differential Equations for a relaible way of estimating the parameters you want.

There are several other examples as well (most of which I wrote). I can likely help you estimate the parameters if you describe what you want to do more clearly.

[Lenny Guitar]

Thank you so much. I need to crete model, where will be this three function. As per below, I wrote the functions for two catalysts in self-contained scripts, but I would need to combine it into one function. The equations are the same, the only thing that is different is the velocity constant k. The function for the third catalyst is exactly the same, just different is in k, but this is stated in the first post.

t_span=[0 60];

%........A B C D E F G H I

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

[t,c]=ode45(@kineticsC,t_span,c0);

plot(t,c)

xlabel('Time,[min]');

ylabel('Concentration,[mol/cm^3]');

legend('cct-CDT','ctt-CDT','ttt-CDT','cc-CDD', 'ct-CDD', 'tt-CDD', 'c-CDE', 't-CDE', 'CDA')

function dz=kinetics(t,c)

k1=0.105;

k2=0.097;

k3=0.065;

k4=0.071;

k5=0.060;

k6=0.056;

k7=0.065;

k8=0.250;

k9=0.061;

k10=0.07;

k11=0.011;

k12=0.120;

dz=zeros(9,1);

r1=k1*c(1);

r2=k2*c(1);

r3=k3*c(2);

r4=k4*c(2);

r5=k5*c(3);

r6=k1*c(1)-k6*c(4);

r7=k2*c(1)+k3*c(2)-k7*c(5)-k8*c(5)-k9*c(5);

r8=k4*c(2)+k5*c(3)-k10*c(6);

r9=k6*c(4)+k7*c(5)-k11*c(7);

r10=k9*c(5)+k10*c(6)-k12*c(8);

r11=k11*c(7)+k8*c(5)+k12*c(8);

dz(1)=-r1-r2;

dz(2)=-r3-r4;

dz(3)=-r5;

dz(4)=r6;

dz(5)=r7;

dz(6)=r8;

dz(7)=r9;

dz(8)=r10;

dz(9)=r11;

end

t_span=[0 300];

%........A B C D E F G H I

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

[t,c]=ode45(@kineticsCaCO3,t_span,c0);

plot(t,c)

xlabel('Time,[min]');

ylabel('Concentration,[mol/cm^3]');

legend('cct-CDT','ctt-CDT','ttt-CDT','cc-CDD', 'ct-CDD', 'tt-CDD', 'c-CDE', 't-CDE', 'CDA')

function dz=kineticsCaCO3(t,c)

k1=0.015;

k2=0.007;

k3=0.01;

k4=0.008;

k5=0.025;

k6=0.016;

k7=0.005;

k8=0.003;

k9=0.007;

k10=0.015;

k11=0.001;

k12=0.035;

dz=zeros(9,1);

r1=k1*c(1);

r2=k2*c(1);

r3=k3*c(2);

r4=k4*c(2);

r5=k5*c(3);

r6=k1*c(1)-k6*c(4);

r7=k2*c(1)+k3*c(2)-k7*c(5)-k8*c(5)-k9*c(5);

r8=k4*c(2)+k5*c(3)-k10*c(6);

r9=k6*c(4)+k7*c(5)-k11*c(7);

r10=k9*c(5)+k10*c(6)-k12*c(8);

r11=k11*c(7)+k8*c(5)+k12*c(8);

dz(1)=-r1-r2;

dz(2)=-r3-r4;

dz(3)=-r5;

dz(4)=r6;

dz(5)=r7;

dz(6)=r8;

dz(7)=r9;

dz(8)=r10;

dz(9)=r11;

end

[Star Strider]

This appears to be slightly more efficient, and uses the same axes to present both integrations:

km = [  0.1050    0.0150
        0.0970    0.0070
        0.0650    0.0100
        0.0710    0.0080
        0.0600    0.0250
        0.0560    0.0160
        0.0650    0.0050
        0.2500    0.0030
        0.0610    0.0070
        0.0700    0.0150
        0.0110    0.0010
        0.1200    0.0350];

t_span=linspace(0, 60, 120);
%........A B C D E F G H I
c0=[1 1 1 0 0 0 0 0 0];
for k = 1:size(km,2)
    [t{k},c{k}]=ode45(@(t,c)kinetics(t,c,km(:,k)),t_span,c0);
end

function dz=kinetics(t,c,kv)
k1=kv(1);
k2=kv(2);
k3=kv(3);
k4=kv(4);
k5=kv(5);
k6=kv(6);
k7=kv(7);
k8=kv(8);
k9=kv(9);
k10=kv(10);
k11=kv(11);
k12=kv(12);
dz=zeros(9,1);
r1=k1*c(1);
r2=k2*c(1);
r3=k3*c(2);
r4=k4*c(2);
r5=k5*c(3);
r6=k1*c(1)-k6*c(4);
r7=k2*c(1)+k3*c(2)-k7*c(5)-k8*c(5)-k9*c(5);
r8=k4*c(2)+k5*c(3)-k10*c(6);
r9=k6*c(4)+k7*c(5)-k11*c(7);
r10=k9*c(5)+k10*c(6)-k12*c(8);
r11=k11*c(7)+k8*c(5)+k12*c(8);
dz(1)=-r1-r2;
dz(2)=-r3-r4;
dz(3)=-r5;
dz(4)=r6;
dz(5)=r7;
dz(6)=r8;
dz(7)=r9;
dz(8)=r10;
dz(9)=r11;
end

figure
hold on
for k = 1:size(km,2)
    plot3(ones(size(t{k}))*(k-1)*0.01,t{k},c{k})
    xlabel('Time,[min]');
    ylabel('Concentration,[mol/cm^3]');
    legend('cct-CDT','ctt-CDT','ttt-CDT','cc-CDD', 'ct-CDD', 'tt-CDD', 'c-CDE', 't-CDE', 'CDA')
    text([0 0.01], [25 25], [1 1], {'C','CaCO_3'})
end
hold off
grid on
view(60,30)

producing:

That should get you started. Make appropriate changes to get the result you want.

[Lenny Guitar]

Firstly, I would like to tjank you. But there is still problem with curly bracket in ode45 and plot3. It said error:Brace indexing is not supported for variables of this type. Is there any possibility, how to rewrite t{k},c{k} to not to report error, please. Thank you so much.

km = [  0.1050    0.0150    0.090 
        0.0970    0.0070    0.078
        0.0650    0.0100    0.070
        0.0710    0.0080    0.058
        0.0600    0.0250    0.125
        0.0560    0.0160    0.045
        0.0650    0.0050    0.069
        0.2500    0.0030    0.007
        0.0610    0.0070    0.073
        0.0700    0.0150    0.080
        0.0110    0.0010    0.006
        0.1200    0.0350    0.039];
t_span=linspace(0, 60, 120);
%........A B C D E F G H I
c0=[1 1 1 0 0 0 0 0 0];
for k = 1:size(km,3)
    [t{k},c{k}]=ode45(@(t,c)kinetics(t,c,km(:,k)),t_span,c0);
end
function dz=kinetics(t,c,kv)
k1=kv(1);
k2=kv(2);
k3=kv(3);
k4=kv(4);
k5=kv(5);
k6=kv(6);
k7=kv(7);
k8=kv(8);
k9=kv(9);
k10=kv(10);
k11=kv(11);
k12=kv(12);
dz=zeros(9,1);
r1=k1*c(1);
r2=k2*c(1);
r3=k3*c(2);
r4=k4*c(2);
r5=k5*c(3);
r6=k1*c(1)-k6*c(4);
r7=k2*c(1)+k3*c(2)-k7*c(5)-k8*c(5)-k9*c(5);
r8=k4*c(2)+k5*c(3)-k10*c(6);
r9=k6*c(4)+k7*c(5)-k11*c(7);
r10=k9*c(5)+k10*c(6)-k12*c(8);
r11=k11*c(7)+k8*c(5)+k12*c(8);
dz(1)=-r1-r2;
dz(2)=-r3-r4;
dz(3)=-r5;
dz(4)=r6;
dz(5)=r7;
dz(6)=r8;
dz(7)=r9;
dz(8)=r10;
dz(9)=r11;
end
figure
hold on
for k = 1:size(km,2)
    plot3(ones(size(t{k}))*(k-1)*0.01,t{k},c{k})
    xlabel('Time,[min]');
    ylabel('Concentration,[mol/cm^3]');
    legend('cct-CDT','ctt-CDT','ttt-CDT','cc-CDD', 'ct-CDD', 'tt-CDD', 'c-CDE', 't-CDE', 'CDA')
    text([0 0.01], [25 25], [1 1], {'C','CaCO_3','Al2O3'})
end
hold off
grid on
view(60,30)

[Star Strider]

Note that ‘km’ has only 2 dimensions (row, column), not 3 (row, column, page), and that is the source of the error you are seeing.

This:

for k = 1:size(km,2)
    [t{k},c{k}]=ode45(@(t,c)kinetics(t,c,km(:,k)),t_span,c0);
end

works, and presents all variables:

No other changes to the code itself are necessary.

Note — You can add as many columns as you want to ‘km’, providing they have the same number of rows. My code automatically adapts to the column-size of ‘km’.

The only other change is to the text call, that now enlarges to include the third simulation and becomes:

text([0 0.01 0.02], [25 25 25], [1 1 1], {'C','CaCO_3','Al_2O_3'})

to produce the figure I plotted here.

—————

(As an undergraduate Chemistry major — albeit almost in the phlogiston era — I find this interesting! I would welcome a description of what you are simulating here, simply for my edification.)

.

[Lenny Guitar]

Thank you so much for your help. I send you a description of these problem. I had to create a program that will monitor the concentrations of the individual components of the reaction mixture. The reaction scheme is demonstrated as follows.

Based on the graphs, I have to compare the selectivity with respect to the yield of t-CDE.

[Star Strider]

If you only want to see the effects on ‘t-CDE’, add this plot at the end of the current code:

cm = cell2mat(c);
figure
plot(t_span, cm(:,8:9:end))
legend({'C','CaCO_3','Al_2O_3'})
grid
xlabel('Time,[min]');
ylabel('t-CDE Concentration,[mol/cm^3]');

.

[Star Strider]

If my Answer helped you solve your problem, please Accept it!

.