Need Help with surface plot. Getting error that dimensions of S do not match

For the reaction A->B->C, with k1 and k2 as rate constants, corresponding to a prefactor “c” of 10^-2 per second, find the ratio of [B] to [C] at 1 second, for an A0 of 1 Mole per liter, if T varies from 100 to 1000K, and Ea1 and Ea2 vary from 10 to 100 kJ/mol.

Make a surface plot to show the results,

1)

With the selectivity ratio “S”, as the z/color axis, T as the x axis, and the ratio of Ea1 to Ea2 on the y axis

2)

With the selectivity ratio “S”, as the z/color axis, T as the x axis, and the difference between Ea1 to Ea2 on the y axis

This is the code we have so far and the plot is not coming up.

clear all
clc
close all
 
tspan=logspace(-10,6,(16*10)+1)
y0=[1;0;0]% IC of Concentration
 
Ea1=linspace(10,15);
Ea2=linspace(10,15);
T=linspace(100,200);
 
leng1=length(T);
leng2=length(Ea1);
leng3=length(Ea2);
 
for i=1:leng1
    for j=1:leng2
            for k=1:leng3
            Ea11=Ea1(i)
            Ea22=Ea2(j)
            TT=T(k)
            
            [t,CC]=ode23(@(t,CC) dosimplereaction(t,CC,Ea1,Ea2,T), tspan,y0); % the code for this function is shown on the bottom
            
            end
    end
end
      
 
A=CC(:,1); B=CC(:,2); C=CC(:,3);
S=B/C
 
Y=Ea1/Ea2;
 
figure(2)
hold on
surf(T,Y,S)
shading 'interp'
colormap('jet(256)')
ylabel 'Ea1/Ea2'
xlabel 'T in Kelvin'
colorbar
caxis([0 10])

% A-B-C reaction, simple 1st order , irreversible
function dCCdt= dosimplereaction(t,CC,Ea1,Ea2,T)
 
prefact=10^-2;
A=CC(1); B=CC(2); C=CC(3);
R=8.314;
 
term1=prod([R,T]);
term2a=prod([-Ea1,term1^-1]);
term2b=prod([-Ea2,term1^-1]);
term3a=prod([prefact,exp(term2a)]);
term3b=prod([prefact,exp(term2b)]);
 
k1=term3a;
k2=term3b;
 
dCCdt(1)=-k1*A;
dCCdt(2)=k1*A-k2*B;
dCCdt(3)=k2*B;
 
dCCdt=dCCdt';
end

[KSSV]

You cannot do like that.....you need to have idea on storing the the data into a matrix and use surf at the end.

[Emily]

How would you go about storing the data in a matrix?

Answer 1

[Walter Roberson]

clear all
clc
close all
 
tspan=logspace(-10,6,(16*10)+1)
y0=[1;0;0]% IC of Concentration
 
Ea1=linspace(10,15);
Ea2=linspace(10,15);
T=linspace(100,200);
 
leng1=length(T);
leng2=length(Ea1);
leng3=length(Ea2);
 
for i=1:leng1
    for j=1:leng2
            for k=1:leng3
            Ea11=Ea1(i)
            Ea22=Ea2(j)
            TT=T(k)
            
            [t,CC]=ode23(@(t,CC) dosimplereaction(t,CC,Ea1,Ea2,T), tspan,y0); % the code for this function is shown on the bottom
            
            end
    end
end
      
 
A=CC(:,1); B=CC(:,2); C=CC(:,3);
S = B ./ C;         %FIXED
 
Y = Ea1 ./ Ea2;     %FIXED
 
figure(2)
hold on
surf(T,Y,S)
shading 'interp'
colormap('jet(256)')
ylabel 'Ea1/Ea2'
xlabel 'T in Kelvin'
colorbar
caxis([0 10])

% A-B-C reaction, simple 1st order , irreversible
function dCCdt= dosimplereaction(t,CC,Ea1,Ea2,T)
 
prefact=10^-2;
A=CC(1); B=CC(2); C=CC(3);
R=8.314;
 
term1=prod([R,T]);
term2a=prod([-Ea1,term1^-1]);
term2b=prod([-Ea2,term1^-1]);
term3a=prod([prefact,exp(term2a)]);
term3b=prod([prefact,exp(term2b)]);
 
k1=term3a;
k2=term3b;
 
dCCdt(1)=-k1*A;
dCCdt(2)=k1*A-k2*B;
dCCdt(3)=k2*B;
 
dCCdt=dCCdt';
end

[Chelsea]

Thank you so much

[Walter Roberson]

tspan=logspace(-10,6,(16*10)+1)

That is a pretty big time span! 10^-10 to 10^6 in 161 steps!

for i=1:leng1
    for j=1:leng2
            for k=1:leng3
            Ea11=Ea1(i)
            Ea22=Ea2(j)
            TT=T(k)
            
            [t,CC]=ode23(@(t,CC) dosimplereaction(t,CC,Ea1,Ea2,T), tspan,y0); % the code for this function is shown on the bottom
            
            end
    end
end

You do all of those calculations, but ever time you throw away all of the results, except for the very last one when i=len1, j=len2, k=len3 . Why are you wasting time doing all of those ode23 calls when you are only interested in the results of the last one??? And why do you assign to TT but not use the value of TT ?

[Chelsea]

We tried using your correction but there is still an error with the surface plot. It says that S is a vector when it must be a matrix