Difference between revisions of "Collection and analysis of rate data - 2013"

From Introduction to Reactor Design: 3K4
Jump to navigation Jump to search
m
Line 59: Line 59:
# Does not quite look like a straight line, but let us fit a least-squares model anyway
# Does not quite look like a straight line, but let us fit a least-squares model anyway


mod.first <- lm( log(CA[1] / CA) ~ t)
mod.first <- lm( log(CA[1] / CA) ~ t + 0) # force intercept to be zero
summary(mod.first)
summary(mod.first)
#Call:
#Call:
#lm(formula = log(CA[1]/CA) ~ t)
#lm(formula = log(CA[1]/CA) ~ t + 0)
#
#
#Residuals:
#Residuals:
#       1       2       3       4       5       6  
#         1         2         3         4         5         6  
#-0.09425 0.02134 0.10646 0.05647 -0.04646 -0.04357
# 1.874e-16 8.988e-02 1.493e-01 7.361e-02 -5.502e-02 -7.784e-02
#
#
#Coefficients:
#Coefficients:
#           Estimate Std. Error t value Pr(>|t|)  
#   Estimate Std. Error t value Pr(>|t|)     
#(Intercept) 0.0942475  0.0604645  1.559 0.194066    
#t 0.0042174 0.0002555   16.51 1.49e-05 ***
#t           0.0037033 0.0003994   9.272 0.000753 ***
#---
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  
#
#
#Residual standard error: 0.08354 on 4 degrees of freedom
#Residual standard error: 0.09474 on 5 degrees of freedom
#Multiple R-squared: 0.9555, Adjusted R-squared: 0.9444
#Multiple R-squared: 0.982, Adjusted R-squared: 0.9784
 
#F-statistic: 272.5 on 1 and 5 DF,  p-value: 1.489e-05


# Now assume a second order system, rA = -k CA^2
# Now assume a second order system, rA = -k CA^2
Line 83: Line 83:
plot(t, 1/CA)
plot(t, 1/CA)
grid()
grid()
abline(mod.first)


# Definitely more linear. Let us fit the least squares model to check:
# Definitely more linear. Let us fit the least squares model to check:
Line 107: Line 108:


So we got a higher \(R^2\) on the second order model. The reaction rate can be considered second order. What is the reaction rate constant, \(k_A\) for this system?
So we got a higher \(R^2\) on the second order model. The reaction rate can be considered second order. What is the reaction rate constant, \(k_A\) for this system?
=== Using MATLAB or Octave ===
<syntaxhighlight lang="MATLAB">
t = [0, 50, 100, 150, 200, 250];    % [minutes]
CA = [154, 114, 87, 76, 70, 58];    % [mol/m^3]
% Plot the data first. See the expected trend. It it decreasing nonlinearly, so
% it is either first or second order (but not zeroth order)
figure
plot(t, CA, 'o')
grid on
% Some template code for the linear regression
x = t;
y = log(CA(1)./CA);
n = numel(y);
% Assume it is a first order, rA = -k CA; so plot log(CA0/CA) against t
figure
plot(x, y, 'o')
grid on
% Does not quite look like a straighTSt line, but let us fit a least-squares model anyway
% Fit the model y = b x + 0 (i.e. force the intercept to zero)
X = [x(:)];
y = y(:);              % make sure it is a column vector     
b = inv(X'*X)*X'*y;    % recall this from the 3E4 course
%b = regress(y, X);    % only if you have the Statistics Toolbox in MATLAB
% Plot the data along with the fitted line:
figure
plot(x, y, 'o')
hold('on')
grid('on')
plot(x, X*b, 'r')
xlabel('x values')
ylabel('y values')
legend({'Original data', 'Fitted line'}, 'Location', 'Best')
% Additional calculations
predicted = X*b;
resids = y - predicted;            % resids = e = y - Xb
RSS = resids' * resids;            % residual sum of squares
TSS = sum((y - mean(y)).^2);      % total sum of squares
R2_first_order = 1 - RSS/TSS
% Now try a second order model: -rA = k CA^2; so plot 1/CA against t
plot(t, 1./CA, 'o')
grid('on')
% Definitely more linear. Let us fit the least squares model,
% 1/CA = 1/CA0 + kt, or y = b0 + b1 x, where y = 1/CA and x = t, and there is an intercept
% Some template code for the linear regression
x = t;
y = 1./CA;
n = numel(y);
X = [ones(n,1) x(:)];
y = y(:);              % make sure it is a column vector     
b = inv(X'*X)*X'*y;    % recall this from the 3E4 course
% b = regress(y, X);    % only if you have the Statistics Toolbox in MATLAB
% Plot the data along with the fitted line:
figure
plot(x, y, 'o')
hold('on')
grid('on')
plot(x, X*b, 'r')
xlabel('x values')
ylabel('y values')
legend({'Original data', 'Fitted line'}, 'Location', 'Best')
% Additional calculations
predicted = X*b;
resids = y - predicted;            % resids = e = y - Xb
RSS = resids' * resids;            % residual sum of squares
TSS = sum((y - mean(y)).^2);      % total sum of squares
R2_second_order = 1 - RSS/TSS
</syntaxhighlight>

Revision as of 16:08, 4 March 2013

Class date(s): 28 February
Download video: Link (plays in Google Chrome) [198 M]

  • F2011: Chapter 7
  • F2006: Chapter 5

28 February 2013 (07C)

R code for the example covered in class:

t <- c(0, 50, 100, 150, 200, 250)    # [minutes]
# or you can write t <- seq(0, 250, 50)

CA <- c(154, 114, 87, 76, 70, 58)  # [mol/m^3]

# Plot the data first. See the expected trend. It it decreasing nonlinearly, so it is either first or second order (but not zeroth order)
plot(t, CA)
grid()

# Assume it is a first order, rA = -k CA
# So plot log(CA0/CA) against t
plot(t, log(CA[1]/CA))
grid()

# Does not quite look like a straight line, but let us fit a least-squares model anyway

mod.first <- lm( log(CA[1] / CA) ~ t + 0)  # force intercept to be zero
summary(mod.first)

#Call:
#lm(formula = log(CA[1]/CA) ~ t + 0)
#
#Residuals:
#         1          2          3          4          5          6 
# 1.874e-16  8.988e-02  1.493e-01  7.361e-02 -5.502e-02 -7.784e-02 
#
#Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
#t 0.0042174  0.0002555   16.51 1.49e-05 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
#
#Residual standard error: 0.09474 on 5 degrees of freedom
#Multiple R-squared: 0.982,	Adjusted R-squared: 0.9784 
#F-statistic: 272.5 on 1 and 5 DF,  p-value: 1.489e-05 

# Now assume a second order system, rA = -k CA^2
# So plot 1/CA against t
plot(t, 1/CA)
grid()
abline(mod.first)

# Definitely more linear. Let us fit the least squares model to check:
mod.second <- lm( I(1/CA) ~ t)
summary(mod.second)
#Call:
#lm(formula = I(1/CA) ~ t)
#
#Residuals:
#         1          2          3          4          5          6 
#-2.751e-04 -5.219e-05  6.146e-04  2.227e-04 -7.051e-04  1.951e-04 
#
#Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 6.769e-03  3.692e-04   18.33 5.21e-05 ***
#t           4.111e-05  2.439e-06   16.86 7.26e-05 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
#
#Residual standard error: 0.0005101 on 4 degrees of freedom
#Multiple R-squared: 0.9861,	Adjusted R-squared: 0.9826 
#F-statistic: 284.2 on 1 and 4 DF,  p-value: 7.259e-05

So we got a higher \(R^2\) on the second order model. The reaction rate can be considered second order. What is the reaction rate constant, \(k_A\) for this system?

Using MATLAB or Octave

t = [0, 50, 100, 150, 200, 250];    % [minutes]
CA = [154, 114, 87, 76, 70, 58];    % [mol/m^3]

% Plot the data first. See the expected trend. It it decreasing nonlinearly, so 
% it is either first or second order (but not zeroth order)
figure
plot(t, CA, 'o')
grid on

% Some template code for the linear regression
x = t;
y = log(CA(1)./CA);
n = numel(y);

% Assume it is a first order, rA = -k CA; so plot log(CA0/CA) against t
figure
plot(x, y, 'o')
grid on

% Does not quite look like a straighTSt line, but let us fit a least-squares model anyway
% Fit the model y = b x + 0 (i.e. force the intercept to zero) 
X = [x(:)]; 
y = y(:);              % make sure it is a column vector       
b = inv(X'*X)*X'*y;    % recall this from the 3E4 course 
%b = regress(y, X);     % only if you have the Statistics Toolbox in MATLAB
 
% Plot the data along with the fitted line:
figure
plot(x, y, 'o')
hold('on')
grid('on')
plot(x, X*b, 'r')
xlabel('x values')
ylabel('y values')
legend({'Original data', 'Fitted line'}, 'Location', 'Best')
 
% Additional calculations
predicted = X*b;
resids = y - predicted;            % resids = e = y - Xb
RSS = resids' * resids;            % residual sum of squares
TSS = sum((y - mean(y)).^2);       % total sum of squares
R2_first_order = 1 - RSS/TSS


% Now try a second order model: -rA = k CA^2; so plot 1/CA against t
plot(t, 1./CA, 'o')
grid('on')

% Definitely more linear. Let us fit the least squares model,
% 1/CA = 1/CA0 + kt, or y = b0 + b1 x, where y = 1/CA and x = t, and there is an intercept
% Some template code for the linear regression
x = t;
y = 1./CA;
n = numel(y);
X = [ones(n,1) x(:)]; 
y = y(:);              % make sure it is a column vector       
b = inv(X'*X)*X'*y;    % recall this from the 3E4 course 
% b = regress(y, X);     % only if you have the Statistics Toolbox in MATLAB

% Plot the data along with the fitted line:
figure
plot(x, y, 'o')
hold('on')
grid('on')
plot(x, X*b, 'r')
xlabel('x values')
ylabel('y values')
legend({'Original data', 'Fitted line'}, 'Location', 'Best')
 
% Additional calculations
predicted = X*b;
resids = y - predicted;            % resids = e = y - Xb
RSS = resids' * resids;            % residual sum of squares
TSS = sum((y - mean(y)).^2);       % total sum of squares
R2_second_order = 1 - RSS/TSS