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?

Question

Lenny Guitar 2021-1-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/727483-kinetics-model-why-i-dont-get-3-graphs-of-these-modeling-but-i-get-just-one-for-c-catalyst-is-ther

评论： Star Strider 2021-1-27

采纳的回答： Star Strider

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2021-1-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/727483-kinetics-model-why-i-dont-get-3-graphs-of-these-modeling-but-i-get-just-one-for-c-catalyst-is-ther#answer_606703

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Lenny Guitar 2021-1-26

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 2021-1-26

在 MATLAB Online 中打开

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 2021-1-26

编辑：Lenny Guitar 2021-1-26

在 MATLAB Online 中打开

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 2021-1-27

在 MATLAB Online 中打开

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]');