MATLAB: Hello lovely people, i have been able to write this code which corresponds to the graphs in literature. I however wanna make a second graph of R2,R3,R6 vs R1. How do i make the plot command? Pls check the images attached.

graph generateplot of parameters

</matlabcentral/answers/uploaded_files/100288/Graphs%20from%20literature.PNG
FUNCTION FILE
 function [rate_out] = COMPREAC(t,x)
% Example  COMPREAC 
%COMPLEX REACTION SCHEME BETWEEN FORMALDEHYDE AND SODIUM PARA PHENOL SULPHONATE%
 %Second order rate constants, m3/kmol s%
K1=0.16;
K2=0.05;
K3=0.15;
K4=0.14; 
K5=0.03;
K6=0.058;
K7=0.05;
K8=0.05;
 %The concentration of the species are first assigned in the vector format
CA=x(1);
CB=x(2);
CC=x(3);
CD=x(4);
CF=x(5);
CG=x(6);
 %Kinetic rates, kmol/m3 s%
R1= K1*CA*CB;  
R2= K2*CA*CC;  
R3= K3*CC*CD ; 
R4= K4*CB*CD ; 
R5= K5*CC*CC ; 
R6= K6*CC*CB ; 
R7= K7*CA*CG ; 
R8= K8*CA*CF;
 %Batch mass balances% 
dCAdt= -R1-R2-R7-R8;
dCBdt= -R1-R4-R6 ;  
dCCdt= R1-R2-R3-2*R5-R6 ; 
dCDdt= R2-R3-R4;
dCEdt= R3+R8;   
dCFdt= R4+R5+R7-R8; 
dCGdt= R6-R7;
 RA=dCAdt;
RB=dCBdt;
RC=dCCdt;
RD=dCDdt;
RE=dCEdt;
RF=dCFdt;
RG=dCGdt;
 % The output of the function reaction rates are saved in vector format 
rate_out=[RA; RB; RC; RD; RE; RF; RG];
  end
 %Initial concentrations, kmol/m3}
 C0= [0.15 0.1 0 0 0 0 0];
RUN FILE
 % Time span
tspan=[0 1000];
 % Run ODE solver
 [t,x]=ode15s(@COMPREAC,tspan, C0);
 %Graph Plot
 plot(t,x);
 xlabel('Time[s]');
 ylabel('CA,CB,CC,CD,CE,CF,CG [kmol/m3 ]');
 legend('RA','RB','RC','RD','RE','RF','RG');

Best Answer

In the run file, after the line
[t,x]=ode15s(@COMPREAC,tspan, C0);
recalculate R1,R2,R3 and R6 from x in a for-loop:
K1 = 0.16;
K2 = 0.05;
K3 = 0.15;
K6 = 0.058;
n = numel(t);
R1 = zeros(n,1);
R2 = zeros(n,1);
R3 = zeros(n,1);
R4 = zeros(n,1);
for i=1:numel(t)  
 %The concentration of the species are first assigned in the vector format
  CA = x(i,1);
  CB = x(i,2);
  CC = x(i,3);
  CD = x(i,4);
%Kinetic rates, kmol/m3 s%
  R1(i) = K1*CA*CB;  
  R2(i) = K2*CA*CC;  
  R3(i) = K3*CC*CD ; 
  R6(i) = K6*CC*CB ;
end
After that, use plot:
plot(R1,R2,'-o',R1,R3,'-.',...)
Best wishes
Torsten.
Related Question