Difference between revisions of "Assignment 4 - 2010"

From Process Model Formulation and Solution: 3E4
Jump to navigation Jump to search
m (Undo revision 817 by Kevindunn (talk))
m
 
(2 intermediate revisions by the same user not shown)
Line 1: Line 1:
{{OtherSidebar
{{OtherSidebar
| due_dates = 08 November 2010
| due_dates = 09 December 2010
| dates_alt_text = Due date(s)
| dates_alt_text = Due date(s)
| questions_PDF = Assignment-4b-2010.pdf
| questions_PDF = Assignment-4a-2010.pdf
| questions_text_alt = Assignment questions
| questions_text_alt = Assignment questions
| solutions_PDF =  
| solutions_PDF =  
| solutions_text_alt = Solutions, prepared by Ali and Kevin
| solutions_text_alt = Solutions, prepared by Ali and Elliot
| other_instructions = Hand-in at class.
| other_instructions = Hand-in at class.
}}
}}
Line 14: Line 14:
.. |m3| replace:: m\ :sup:`3`
.. |m3| replace:: m\ :sup:`3`


.. rubric:: Questions 1, 2 and 3 will be posted at the end of the semester in assignment 4(a).
Question 1 [2]
 
.. rubric:: Solutions by Kevin and Ali
 
Question 4 [2]
==============
==============


.. Similar to Tutorial 8, 2009
A new type of `thermocouple <http://en.wikipedia.org/wiki/Thermocouple>`_ is being investigated by your group.  These devices produce an *almost* linear voltage  (millivolt) response at different temperatures. In practice though it is used the other way around: use the millivolt reading to predict the temperature.  The process of fitting this linear model is called *calibration*. 


The viscosity of sulphuric acid, :math:`\nu`, varies with purity, :math:`p` in the following manner:
#. Use the following data to calibrate a linear model:


============================= ===== ===== =====  
================= ==== ==== ==== ==== ==== ==== ==== ==== ==== ====
:math:`p` [%]                 20    60    80
Temperature [K]   273  293  313  333  353  373  393  413  433  453
----------------------------- ----- ----- -----  
----------------- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
:math:`\nu` [millipascal]     1.40  5.37 17.4
Reading [mV]   0.01 0.12 0.24 0.38 0.51 0.67 0.84 1.01 1.15 1.31  
============================= ===== ===== =====  
================= ==== ==== ==== ==== ==== ==== ==== ==== ==== ====
Show the linear model and provide the predicted temperature when reading 1.00 mV.
#. Are you satisfied with this model, based on the coefficient of determination (:math:`R^2`) value? 


#. Express :math:`\nu(p)` as a quadratic function using Lagrange interpolating polynomials. Do not simplify the polynomial.
#. What is the model's standard error?  Now, are you satisfied with the model's prediction ability, given that temperatures can usually be recorded to an accuracy of :math:`\pm 0.1` K with most thermocouples.
#. Express :math:`\nu(p)` as a quadratic function using Newton interpolating polynomials. Do not simplify the polynomial.
#. Fit a cubic spline through the data points: clearly show your :math:`{\bf Xa = y}` linear system of equations, then solve them using computer software; finally report your spline coefficients.
#. Use computer software to plot:


* the Newton interpolating polynomial
Question 2 [2]
* the cubic spline,
==============
* and the 3 data points on the same graph.
#. What is the estimated viscosity at :math:`p = 40\%` purity using linear interpolation?
#. Which of the estimation procedures that you used above has the closest estimate to the true value of 2.51 millipascal?


Solution
Batch reactors are often used to calculate reaction rate data, since the dynamic balances are a direct function of the reaction rate:
---------


Redo the solution with these new values. The solution below is for old values.
.. math::


#. Recall that an :math:`n-1` order Lagrange interpolating polynomial, for a set of :math:`n` data points, is given by the following relation:
\frac{dC_{\sf A}(t)}{dt} = r_{\sf A}
assuming constant volume and assuming all energy-related terms are negligible. We are unsure about the order of our reaction rate, :math:`\alpha`, and we also don't know the reaction rate constant, :math:`k` in the expression :math:`r_A = - k C_{\sf A}^{\alpha}`.


.. math::
Integrating the above differential equation from time :math:`t=0` to time :math:`t=t` gives:


P_{n-1}(x) = y_{1}\ell_{1}(x) + y_{2}\ell_{2}(x) + \ldots + y_{n}\ell_{n}(x) = \sum_{i=1}^{n}y_{i}\ell_{i}(x) \qquad \text{where} \qquad \ell_{i}(x) = \prod_{j=1, j\neq i}^{n} \frac{x - x_{i}}{x_{i} - x_{j}}
.. math::


For our system the quadratic Lagrange interpolating polynomial is the following:
\int_{C_{\sf A,0}}^{C_{\sf A}}{ \frac{1}{C_{\sf A}^\alpha} dC_{\sf A}} &= -k \int_{0}^{t}{dt} \\
                          \left.          \frac{C_{\sf A}^{1-\alpha}}{1-\alpha} \right|_{C_{\sf A,0}}^{C_{\sf A}} &= -kt  \qquad \text{where}\,\, \alpha \neq 1


.. math::
The following concentrations are collected as a function of time


P_{2}(x) &= y_{1}\ell_{1}(x) + y_{2}\ell_{2}(x) + \ldots + y_{n}\ell_{n}(x) \\
====================== ==== ==== ==== ==== ==== ==== ==== ==== ====
        &= 1.40 \cdot \frac{(x-60)(x-80)}{(20-60)(20-80)}
Time [hours]          0    0.5  1   2   3    5    7.5  10  17
      + 5.37 \cdot \frac{(x-20)(x-80)}{(60-20)(60-80)}
---------------------- ---- ---- ---- ---- ---- ---- ---- ---- ----
      + 17.4 \cdot \frac{(x-20)(x-60)}{(80-20)(80-60)}\\
Concentration [mol/L]  3.1.5  0.92 0.47 0.29 0.16 0.09 0.06 0.03
====================== ==== ==== ==== ==== ==== ==== ==== ==== ====
#. Recall that an :math:`n-1` order Newton interpolating polynomial, for a set of :math:`n` data points, is given by the following relation:


.. math::
P_{n-1}(x) = b_{0}+b_{1}(x-x_{1})+b_{2}(x-x_{1})(x-x_{2})+ \ldots +b_{n-1}(x-x_{1})\cdots(x-x_{n-1})
The coefficients, :math:`b_i` are best calculated in table form, known as "Newton's divided differences":


.. math::
#. Use these concentration-time data to determine the values of :math:`\alpha` and :math:`k`, using a least squares model.
\begin{array}{cccc}
\hline \hline
p_{i} & \nu_{i}      & \nu[p_{i},p_{i-1}]                      & \nu[p_{i}, p_{i-1}, p_{i-2}] \\ \hline \hline
20    & \textbf{1.40} & \frac{5.37-1.40}{60-20}=\textbf{0.09925} & \frac{0.6015-0.09925}{80-20}=\textbf{0.008371} \\
60    & 5.37          & \frac{17.4-5.37}{80-60}=0.6015          &  \\
80    & 17.4          &  &  \\ \hline \hline
\end{array}


#. Once you have the rate constants, simulate the batch reactor (use the ``ode45`` or Python integrator) and superimpose the data from the above table on your concentration-time trajectory.  Do they agree?  Which regions of the plot have the worst agreement (this is where you would collect more data in the future).


So the Newton interpolating polynomial is:


.. math::
Question 3 [3]
==============


P_{2}(x) = 1.40 + 0.09925(x-x_{1}) + 0.008371(x-x_{1})(x-x_{2})
The yield from your lab-scale bioreactor, :math:`y`, is a function of reactor temperature, impeller speed and reactor type (one with with baffles and one without). You have collected these data from various experiments.
#. The cubic spline through 3 data points consists of two cubic polynomials -- one connecting points 1 and 2, and another connecting points 2 and 3.
.. csv-table::
  :header: Temp = :math:`T` [°C], Speed = :math:`S` [RPM], Baffles = :math:`B` [Yes/No], Yield = :math:`y` [g]
  :widths: 30, 30, 30, 30


.. math::
82,      4300,      No,      51
90,      3700,      Yes,    30
88,      4200,      Yes,    40
86,      3300,      Yes,    28
80,      4300,      No,      49
78,      4300,      Yes,    49
82,      3900,      Yes,    44
83,      4300,      No,      59
65,      4200,      No,      61
73,      4400,      No,      59
60,      4400,      No,      57
60,      4400,      No,      62
101,    4400,      No,      42
92,      4900,      Yes,    38
#. Use a built-in software function in MATLAB or Python (please don't use Excel) to fit a linear model that predicts the bioreactor yield from the above variables. Your model should be of the form:


s_3(x) = \left\{
\begin{array}{l}
a_1 + b_1x + c_1x^2 + d_1x^3 \qquad\,\,x_1 \leq x \leq x_2 \\
a_2 + b_2x + c_2x^2 + d_2x^3 \qquad\,\,x_2 \leq x \leq x_3
\end{array}
\right.
So we have 8 unknown coefficients in these 2 polynomials.  The 8 linear equations that can be used to solve for the coefficients are:
.. math::
.. math::
\begin{array}{ll}
\hat{y} = a_0 + a_T T + a_S S + a_B B
\text{Starting knot} & s_3(x_1) = y_1 = a_1 + b_1x_1 + c_1x_1^2 + d_1x_1^3\\
\text{End knot} & s_3(x_3) = y_3 = a_2 + b_2x_3 + c_2x_3^2 + d_2x_3^3\\
\text{Continuity} &  s_3(x_2) = \left\{
\begin{array}{l} y_2 = a_1 + b_1x_2 + c_1x_2^2 + d_1x_2^3 \\
y_2 = a_2 + b_2x_2 + c_2x_2^2 + d_2x_2^3
\end{array}
      \right.\\
\text{1st derivative continuity} &
s'_3(x_2) = b_1 + 2c_1x_2 + 3d_1x_2^2 = b_2 + 2c_2x_2 + 3d_2x_2^2 \\
\text{2nd derivative continuity} &
s''_3(x_2) = 2c_1 + 6d_1x_2 = 2c_2 + 6d_2x_2 \\
\text{$s''_3(x)=0$ at the ends}  &
s''_3(x_1) = 0 = 2c_1 + 6d_1x_1 \\
        &
s''_3(x_3) = 0 = 2c_2 + 6d_2x_3
\end{array}
Expressing this as a linear system of equations, :math:`{\bf Xa = y}`:
*Hint*: convert the No/Yes variables to integers with ``0=No`` and ``1=Yes``.
.. math::
#. Interpret the meaning of each coefficient in the model and comment on the standard error.
 
{\bf Xa} &= {\bf y} \\
#. Calculate and show the :math:`\mathbf{X}^T\mathbf{X}` and :math:`\mathbf{X}^T\mathbf{y}` matrices for this linear model in order to calculate the 4 coefficients in :math:`{\bf a} = \left[a_0, a_T, a_S, a_B\right]` yourself.   Your answer must include the MATLAB or Python code that computes these matrices and the :math:`{\bf a}` solution vector.
\left( \begin{array}{cccccccc}  
 
0 & 0  & 2    & 6x_1  \\
1 & x_1 & x_1^2 & x_1^3  \\
1 & x_2 & x_2^2 & x_2^3  \\
0 & 1  & 2x_2  & 3x_2^2 & 0 & -1  & -2x_2 & -3x_2^2 \\
0 & 0  & 2    & 6x_2  & 0 & 0  & -2    & -6x_2  \\
  &    &      &        & 1 & x_2 & x_2^2 & x_2^3  \\
  &    &      &        & 1 & x_3 & x_3^2 & x_3^3  \\
  &    &      &        &  &    & 2    & 6x_3   
\end{array}\right)
\left(
\begin{array}{c}
a_1 \\
b_1 \\
c_1 \\
d_1 \\
a_2 \\
b_2 \\
c_2 \\
d_2
\end{array}
\right)
&=  
\left(
\begin{array}{c}
0  \\
y_1 \\
y_2 \\
0  \\
0  \\
y_2 \\
y_3 \\
0
\end{array}
\right) \\
\left( \begin{array}{cccccccc}
0 & 0   & 2    & 120    \\
1 & 20  & 400  & 8000  \\
1 & 60  & 3600  & 216000 \\
0 & 1  & 120  & 10800  & 0 & -1  & -120  & -10800  \\
0 & 0  & 2    & 360    & 0 & 0  & -2    & -360  \\
  &    &      &        & 1 & 60  & 3600  & 216000 \\
  &    &      &        & 1 & 80  & 6400  & 512000  \\
  &    &      &        & 0 & 0  & 2    & 480   
\end{array}\right)
\left(
\begin{array}{c}
a_1 \\
b_1 \\
c_1 \\
d_1 \\
a_2 \\
b_2 \\
c_2 \\
d_2
\end{array}
\right)
&=
\left(
\begin{array}{c}
0    \\
1.40 \\
5.37 \\
0    \\
0    \\
5.37 \\
17.4 \\
0
\end{array}
\right)
Solving this system of equations with computer software gives:
.. math::


s_3(x) = \left\{
Question 4 [2]
\begin{array}{l}
==============
1.93 + 0.0579x - 6.28\times10^{-3} x^2 + 1.05\times10^{-4}x^3 \qquad\,\,x_1 \leq x \leq x_2 \\
69.73 - 3.33x + 5.02\times10^{-2} x^2 - 2.09\times10^{-4}x^3 \qquad\,\,x_2 \leq x \leq x_3
\end{array}
\right.
#. A plot of the Lagrange (or Newton) polynomial, the spline, the linear interpolation (not required) superimposed on the 3 original data points is shown below.


.. figure:: ../che3e4/Assignments/Assignment-4/images/q1_plot.png
.. Similar to Tutorial 8, 2009
:scale: 45
:align: center


The code used to generate this plot is similar to that shown on `course software tutorial <http://modelling3e4.connectmv.com/wiki/Software_tutorial/Creating_and_saving_plots>`_.
The viscosity of sulphuric acid, :math:`\nu`, varies with purity, :math:`p` in the following manner:


#. The linear interpolation for :math:`\nu(p=40)` has slope = :math:`\displaystyle \frac{5.37 - 1.40}{60 - 20} = 0.09925` and intercept = :math:`1.40 - 0.09925 \times 20 = -0.585`.  So :math:`\nu(p=40) = 0.09925 \times 40 - 0.585 = 3.385` millipascal.
============================= ===== ===== =====
:math:`p` [%]                20    60    80
----------------------------- ----- ----- -----
:math:`\nu` [millipascal]    1.40  5.37 17.4
============================= ===== ===== =====  


#. The question was a bit unclear; the intention was to compare all estimates of :math:`\nu(p=40)`:
#. Express :math:`\nu(p)` as a quadratic function using Lagrange interpolating polynomials. Do not simplify the polynomial.
#. Express :math:`\nu(p)` as a quadratic function using Newton interpolating polynomials. Do not simplify the polynomial.
#. Fit a cubic spline through the data points: clearly show your :math:`{\bf Xa = y}` linear system of equations, then solve them using computer software; finally report your spline coefficients.
#. Use computer software to plot:


* Lagrange polynomial estimate = :math:`1.40 \cdot \frac{(40-60)(40-80)}{(20-60)(20-80)}+ 5.37 \cdot \frac{(40-20)(40-80)}{(60-20)(60-80)} + 17.4 \cdot \frac{(40-20)(40-60)}{(80-20)(80-60)} = 1.40(0.3333) + 5.37(1.0) + 17.4(-0.333) = 0.0366` millipascals.
* the Newton interpolating polynomial  
* Newton polynomial estimate = :math:`1.40 + 0.09925(40-20) + 0.008371(40-20)(40-60) = 0.0366` millipascals.
* the cubic spline,
* Cubic spline estimate = :math:`1.93 + 0.0579\times 40 - 6.28\times10^{-3} \times 40^2 + 1.05\times10^{-4} \times 40^3 = 0.918` millipascals.
* and the 3 data points on the same graph.
* Linear interpolation = 3.385 millipascal.
#. What is the estimated viscosity at :math:`p` = 40% purity using linear interpolation?
#. Which of the estimation procedures that you used above has the closest estimate to the true value of 2.51 millipascal?
All estimates had large error when compared to the true value of 2.51, because the viscosity is changing fairly rapidly in this neighbourhood.  The estimates would be improved if there was another data point close by.  The Lagrange/Newton estimates are especially poor, since they are negative in the range from approximately 30% to 40% purity.


Question 5 [2]
Question 5 [2]
Line 242: Line 141:


Show your matrix derivation for the linear system of equations, and solve it using computer software.  Plot the cubic spline over the time range 0 to 8 hours.
Show your matrix derivation for the linear system of equations, and solve it using computer software.  Plot the cubic spline over the time range 0 to 8 hours.
Solution
---------
The cubic spline through 5 data points consists of four cubic polynomials -- connecting points 1 and 2, 2 and 3, 3 and 4, and 4 and 5.
.. math::
s_3(x) = \left\{
\begin{array}{l}
a_1 + b_1x + c_1x^2 + d_1x^3 \qquad\,\,x_1 \leq x \leq x_2 \\
a_2 + b_2x + c_2x^2 + d_2x^3 \qquad\,\,x_2 \leq x \leq x_3\\
a_3 + b_3x + c_3x^2 + d_3x^3 \qquad\,\,x_3 \leq x \leq x_4\\
a_4 + b_4x + c_4x^2 + d_4x^3 \qquad\,\,x_4 \leq x \leq x_5
\end{array}
\right.
So we have 16 unknown coefficients in these 4 polynomials.  The 16 linear equations that can be used to solve for the coefficients are:
.. math::
\begin{array}{ll}
\text{Starting knot} & s_3(x_1) = y_1 = a_1 + b_1x_1 + c_1x_1^2 + d_1x_1^3\\
\text{End knot} & s_3(x_5) = y_5 = a_2 + b_2x_5 + c_2x_5^2 + d_2x_5^3\\
\text{Continuity} &  s_3(x_2) = \left\{
\begin{array}{l} y_2 = a_1 + b_1x_2 + c_1x_2^2 + d_1x_2^3 \\
y_2 = a_2 + b_2x_2 + c_2x_2^2 + d_2x_2^3
\end{array}
      \right.\\
\text{Continuity} &  s_3(x_3) = \left\{
\begin{array}{l} y_3 = a_2 + b_2x_3 + c_2x_3^2 + d_2x_3^3 \\
y_3 = a_3 + b_3x_3 + c_3x_3^2 + d_3x_3^3
\end{array}
      \right.\\
\text{Continuity} &  s_3(x_4) = \left\{
\begin{array}{l} y_4 = a_3 + b_3x_4 + c_3x_4^2 + d_3x_4^3 \\
y_4 = a_4 + b_4x_4 + c_4x_4^2 + d_4x_4^3
\end{array}
      \right.\\
\text{1st derivative continuity} &
s'_3(x_2) = b_1 + 2c_1x_2 + 3d_1x_2^2 = b_2 + 2c_2x_2 + 3d_2x_2^2 \\
\text{1st derivative continuity} &
s'_3(x_3) = b_2 + 2c_2x_3 + 3d_2x_3^2 = b_3 + 2c_3x_3 + 3d_3x_3^2 \\
\text{1st derivative continuity} &
s'_3(x_4) = b_3 + 2c_3x_4 + 3d_3x_4^2 = b_4 + 2c_4x_4 + 3d_4x_4^2 \\
\text{2nd derivative continuity} &
s''_3(x_2) = 2c_1 + 6d_1x_2 = 2c_2 + 6d_2x_2 \\
\text{2nd derivative continuity} &
s''_3(x_3) = 2c_2 + 6d_2x_3 = 2c_3 + 6d_3x_3 \\
\text{2nd derivative continuity} &
s''_3(x_4) = 2c_3 + 6d_3x_4 = 2c_4 + 6d_4x_4 \\
\text{$s''_3(x)=0$ at the ends}  &
s''_3(x_1) = 0 = 2c_1 + 6d_1x_1 \\
        &
s''_3(x_5) = 0 = 2c_2 + 6d_2x_5
\end{array}
Expressing this as a linear system of equations, :math:`{\bf Xa = y}` (see the lecture notes for the matrix pattern):
.. math::
{\bf Xa} &= {\bf y} \\
\left( \begin{array}{cccccccccccccccc}
0 & 0  & 2    & 6x_1  & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & x_1 & x_1^2 & x_1^3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & x_2 & x_2^2 & x_2^3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 1  & 2x_2  & 3x_2^2 & 0 & -1  & -2x_2 & -3x_2^2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0  & 2    & 6x_2  & 0 & 0  & -2    & -6x_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0  \\
0  &  0  &  0    &  0    & 1 & x_2 & x_2^2 & x_2^3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0  \\
0  &  0  &  0    &    0    & 1 & x_3 & x_3^2 & x_3^3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0  \\
0 & 0 & 0 & 0 & 0 & 1  & 2x_3  & 3x_3^2 & 0 & -1  & -2x_3 & -3x_3^2 &  0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0  & 2    & 6x_3  & 0 & 0  & -2    & -6x_3 & 0 &  0 & 0 & 0  \\
0 & 0 & 0 & 0 & 0 & 0 &  0  &  0    &    1 & x_3 & x_3^2 & x_3^3 & 0 & 0 & 0 & 0  \\
0 & 0 & 0 & 0 & 0 &  0  &  0    &    0    & 1 & x_4 & x_4^2 & x_4^3 & 0 & 0 & 0 &  0  \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1  & 2x_4  & 3x_4^2 & 0 & -1  & -2x_4 & -3x_4^2  \\
0 &  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0  & 2    & 6x_4  & 0 & 0  & -2    & -6x_4    \\
0 & 0 & 0 & 0 &0 & 0 & 0 & 0 & 0 & 0 &  0  &  0    &    1 & x_4 & x_4^2 & x_4^3  \\
0 & 0 & 0 &  0 &0 & 0 & 0 & 0 & 0 &  0  &  0    &    0    & 1 & x_5 & x_5^2 & x_5^3  \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &  0  & 0      &    0    & 0  &  0  & 2    & 6x_5   
\end{array}\right)
\left(
\begin{array}{c}
a_1 \\
b_1 \\
c_1 \\
d_1 \\
a_2 \\
b_2 \\
c_2 \\
d_2 \\
a_3 \\
b_3 \\
c_3 \\
d_3 \\
a_4 \\
b_4 \\
c_4 \\
d_4
\end{array}
\right)
&=
\left(
\begin{array}{c}
0  \\
y_1 \\
y_2 \\
0  \\
0  \\
y_2 \\
y_3 \\
0  \\
0  \\
y_3 \\
y_4 \\
0  \\
0  \\
y_4 \\
y_5 \\
0  \\
      \end{array}
\right) \\
\left( \begin{array}{cccccccccccccccc}
0 & 0  & 2    & 0  & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 1  & 2  & 3 & 0 & -1  & -2 & -3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0  & 2    & 6  & 0 & 0  & -2    & -6 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0  \\
0  &  0  &  0    &  0    & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0  \\
0  &  0  &  0    &    0    & 1 & 2 & 4 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0  \\
0 & 0 & 0 & 0 & 0 & 1  & 4  & 12 & 0 & -1  & -4 & -12 &  0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0  & 2    & 12  & 0 & 0  & -2    & -12 & 0 &  0 & 0 & 0  \\
0 & 0 & 0 & 0 & 0 & 0 &  0  &  0    &    1 & 2 & 4 & 8 & 0 & 0 & 0 & 0  \\
0 & 0 & 0 & 0 & 0 &  0  &  0    &    0    & 1 & 4 & 16 & 64 & 0 & 0 & 0 &  0  \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1  & 8  & 48 & 0 & -1  & -8 & -48  \\
0 &  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0  & 2    & 24  & 0 & 0  & -2    & -24    \\
0 & 0 & 0 & 0 &0 & 0 & 0 & 0 & 0 & 0 &  0  &  0    &    1 & 4 & 16 & 64  \\
0 & 0 & 0 &  0 &0 & 0 & 0 & 0 & 0 &  0  &  0    &    0    & 1 & 6 & 36 & 216  \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &  0  & 0      &    0    & 0  &  0  & 2    & 36   
\end{array}\right)
\left(
\begin{array}{c}
a_1 \\
b_1 \\
c_1 \\
d_1 \\
a_2 \\
b_2 \\
c_2 \\
d_2 \\
a_3 \\
b_3 \\
c_3 \\
d_3 \\
a_4 \\
b_4 \\
c_4 \\
d_4
\end{array}
\right)
&=
\left(
\begin{array}{c}
0  \\
0.1 \\
0.341 \\
0  \\
0  \\
0.341 \\
1.102 \\
0  \\
0  \\
1.102 \\
4.95 \\
0  \\
0  \\
4.95 \\
11.24 \\
0  \\
      \end{array}
\right) \\
Solving this system of equations with computer software gives:
.. math::
s_3(x) = \left\{
\begin{array}{l}
0.1000 + 0.1457x + 0\times x^2 + 0.0953x^3 \qquad\,\,0 \leq x \leq 1 \\
0.1521  -0.0106x + 0.1562 x^2 + 0.0433x^3 \qquad\,\,1 \leq x \leq 2\\
0.5809 -0.6537x +  0.4778x^2 -0.0103x^3 \qquad\,\,2 \leq x \leq 4\\
3.6951 + -2.9895x + 1.0617x^2 + -0.0590x^3 \qquad\,\,4 \leq x \leq 6
\end{array}
\right.
Now we can estimate the number of cells at time 3, 5, and 7 hours. All we need is to see which of the above cubic polynomials are appropriate for each point:
:math:`2\leq x=3\leq 4\rightarrow` :math:`C(3)\approx 0.5809 -0.6537(3)  +  0.4778(3)^2 -0.0103(3)^3=2.6411`
:math:`4\leq x=5\leq 6\rightarrow` :math:`C(5)\approx 3.6951 + -2.9895(5) + 1.0617(5)^2 + -0.0590(5)^3=7.9180`
The last point 7 is not within the domain of interpolation for any of the intervals. So, we need to extrapolate the last polynomial (which is the closet to :math:`x=7`):
:math:`C(7)\approx 3.6951 + -2.9895(7) + 1.0617(7)^2 + -0.0590(7)^3=14.5620`
The cubic spline is plotted below, followed by the code used to generate the estimate values.
.. figure:: ../che3e4/Assignments/Assignment-4/images/spline_q5.png
:scale: 45
:align: center
**Note**: Polynomial coefficients in MATLAB and Python are sorted from the highest order to the lowest one, which is in the reverse direction we wrote the above polynomials. So, make sure you flip the coefficients (in the decreasing order) when using these polynomials in software.
.. twocolumncode::
    :code1: ../che3e4/Assignments/Assignment-4/code/cubic_spline_5pts.m
    :language1: matlab
    :header1: MATLAB code
    :code2: ../che3e4/Assignments/Assignment-4/code/cubic_spline_5pts_python.py
    :language2: python
    :header2: Python code


Bonus question [0.5]
Bonus question [0.5]
Line 466: Line 146:


Use the cubic spline from the previous question and find the time where the cell count was approximately 10.0 g/L.  Do not solve the equation by hand, but investigate `MATLAB's <http://www.mathworks.com/help/techdoc/ref/roots.html>`_ or `Python's <http://docs.scipy.org/doc/numpy/reference/generated/numpy.roots.html>`_ polynomial root finding function: ``roots(...)``
Use the cubic spline from the previous question and find the time where the cell count was approximately 10.0 g/L.  Do not solve the equation by hand, but investigate `MATLAB's <http://www.mathworks.com/help/techdoc/ref/roots.html>`_ or `Python's <http://docs.scipy.org/doc/numpy/reference/generated/numpy.roots.html>`_ polynomial root finding function: ``roots(...)``
Solution
---------
The value of :math:`C=10` falls between 4.95 (at :math:`t=4`) and 11.24 (at :math:`t=6`) in the table of experimental data . So, we use the last spline polynomial that represents the points in this range:
:math:`C(t)\approx 3.6951 + -2.9895(t) + 1.0617(t)^2 + -0.0590(t)^3=10.0 \rightarrow`
:math:`C(t)\approx -6.3049 + -2.9895(t) + 1.0617(t)^2 + -0.0590(t)^3=0`
This is a simple polynomial root finding problem which can be solved using either
* MATLAB: ``roots([-0.0590, 1.0617, -2.9895, -6.3049])``
* Python: ``np.roots([-0.0590, 1.0617, -2.9895, -6.3049])``
Since it is a cubic function, three roots are returned: :math:`t_1 = 13.7417, t_2 = 5.6336, t_3 = -1.3804`.
In order to decide which one is the desired solution, we again look at the tabulated data. Since the :math:`C` value falls between 4.95 and 11.24, the corresponding :math:`t` value must also fall between 4.0 and 6.0. So, the desired root is :math:`t_2 = 5.6336`.


.. raw:: latex
.. raw:: latex
Line 490: Line 152:
\hrule
\hrule
\begin{center}END\end{center}
\begin{center}END\end{center}
</rst>
</rst>

Latest revision as of 20:18, 2 December 2010

Due date(s): 09 December 2010
Nuvola mimetypes pdf.png (PDF) Assignment questions
Other instructions Hand-in at class.

<rst> <rst-options: 'toc' = False/> <rst-options: 'reset-figures' = False/> .. |m3| replace:: m\ :sup:`3`

Question 1 [2]

==

A new type of `thermocouple <http://en.wikipedia.org/wiki/Thermocouple>`_ is being investigated by your group. These devices produce an *almost* linear voltage (millivolt) response at different temperatures. In practice though it is used the other way around: use the millivolt reading to predict the temperature. The process of fitting this linear model is called *calibration*.

  1. . Use the following data to calibrate a linear model:

================= ==== ==== ==== ==== ==== ==== ==== ==== ==== ==== Temperature [K] 273 293 313 333 353 373 393 413 433 453 ----------------- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- Reading [mV] 0.01 0.12 0.24 0.38 0.51 0.67 0.84 1.01 1.15 1.31 ================= ==== ==== ==== ==== ==== ==== ==== ==== ==== ====


Show the linear model and provide the predicted temperature when reading 1.00 mV.

  1. . Are you satisfied with this model, based on the coefficient of determination (:math:`R^2`) value?
  1. . What is the model's standard error? Now, are you satisfied with the model's prediction ability, given that temperatures can usually be recorded to an accuracy of :math:`\pm 0.1` K with most thermocouples.

Question 2 [2]

==

Batch reactors are often used to calculate reaction rate data, since the dynamic balances are a direct function of the reaction rate:

.. math::

\frac{dC_{\sf A}(t)}{dt} = r_{\sf A}

assuming constant volume and assuming all energy-related terms are negligible. We are unsure about the order of our reaction rate, :math:`\alpha`, and we also don't know the reaction rate constant, :math:`k` in the expression :math:`r_A = - k C_{\sf A}^{\alpha}`.

Integrating the above differential equation from time :math:`t=0` to time :math:`t=t` gives:

.. math::

\int_{C_{\sf A,0}}^{C_{\sf A}}{ \frac{1}{C_{\sf A}^\alpha} dC_{\sf A}} &= -k \int_{0}^{t}{dt} \\

                          \left.          \frac{C_{\sf A}^{1-\alpha}}{1-\alpha} \right|_{C_{\sf A,0}}^{C_{\sf A}} &= -kt  \qquad \text{where}\,\, \alpha \neq 1

The following concentrations are collected as a function of time

================== ==== ==== ==== ==== ==== ==== ==== ====

Time [hours] 0 0.5 1 2 3 5 7.5 10 17


---- ---- ---- ---- ---- ---- ---- ---- ----

Concentration [mol/L] 3.5 1.5 0.92 0.47 0.29 0.16 0.09 0.06 0.03

================== ==== ==== ==== ==== ==== ==== ==== ====

  1. . Use these concentration-time data to determine the values of :math:`\alpha` and :math:`k`, using a least squares model.
  1. . Once you have the rate constants, simulate the batch reactor (use the ``ode45`` or Python integrator) and superimpose the data from the above table on your concentration-time trajectory. Do they agree? Which regions of the plot have the worst agreement (this is where you would collect more data in the future).


Question 3 [3]

==

The yield from your lab-scale bioreactor, :math:`y`, is a function of reactor temperature, impeller speed and reactor type (one with with baffles and one without). You have collected these data from various experiments.

.. csv-table:: :header: Temp = :math:`T` [°C], Speed = :math:`S` [RPM], Baffles = :math:`B` [Yes/No], Yield = :math:`y` [g] :widths: 30, 30, 30, 30

82, 4300, No, 51 90, 3700, Yes, 30 88, 4200, Yes, 40 86, 3300, Yes, 28 80, 4300, No, 49 78, 4300, Yes, 49 82, 3900, Yes, 44 83, 4300, No, 59 65, 4200, No, 61 73, 4400, No, 59 60, 4400, No, 57 60, 4400, No, 62 101, 4400, No, 42 92, 4900, Yes, 38

  1. . Use a built-in software function in MATLAB or Python (please don't use Excel) to fit a linear model that predicts the bioreactor yield from the above variables. Your model should be of the form:

.. math::

\hat{y} = a_0 + a_T T + a_S S + a_B B

*Hint*: convert the No/Yes variables to integers with ``0=No`` and ``1=Yes``.

  1. . Interpret the meaning of each coefficient in the model and comment on the standard error.
  1. . Calculate and show the :math:`\mathbf{X}^T\mathbf{X}` and :math:`\mathbf{X}^T\mathbf{y}` matrices for this linear model in order to calculate the 4 coefficients in :math:`{\bf a} = \left[a_0, a_T, a_S, a_B\right]` yourself. Your answer must include the MATLAB or Python code that computes these matrices and the :math:`{\bf a}` solution vector.


Question 4 [2]

==

.. Similar to Tutorial 8, 2009

The viscosity of sulphuric acid, :math:`\nu`, varies with purity, :math:`p` in the following manner:

======================== ===== =====
math:`p` [%] 20 60 80

----- ----- -----

math:`\nu` [millipascal] 1.40 5.37 17.4
======================== ===== =====
  1. . Express :math:`\nu(p)` as a quadratic function using Lagrange interpolating polynomials. Do not simplify the polynomial.
  2. . Express :math:`\nu(p)` as a quadratic function using Newton interpolating polynomials. Do not simplify the polynomial.
  3. . Fit a cubic spline through the data points: clearly show your :math:`{\bf Xa = y}` linear system of equations, then solve them using computer software; finally report your spline coefficients.
  4. . Use computer software to plot:
	*	the Newton interpolating polynomial 

* the cubic spline, * and the 3 data points on the same graph.

  1. . What is the estimated viscosity at :math:`p` = 40% purity using linear interpolation?
  2. . Which of the estimation procedures that you used above has the closest estimate to the true value of 2.51 millipascal?

Question 5 [2]

===

The following data are collected from a bioreactor experiment, during the growth phase.

============= ===== ===== ===== =====

Time [hours] 0 1.0 2.0 4.0 6.0


----- ----- ----- ----- -----

math:`C` [g/L] 0.1 0.341 1.102 4.95 11.24
============= ===== ===== ===== =====

Fit a natural cubic spline for these data and use it to estimate the number of cells at time 3, 5, and 7 hours.

Show your matrix derivation for the linear system of equations, and solve it using computer software. Plot the cubic spline over the time range 0 to 8 hours.

Bonus question [0.5]

========

Use the cubic spline from the previous question and find the time where the cell count was approximately 10.0 g/L. Do not solve the equation by hand, but investigate `MATLAB's <http://www.mathworks.com/help/techdoc/ref/roots.html>`_ or `Python's <http://docs.scipy.org/doc/numpy/reference/generated/numpy.roots.html>`_ polynomial root finding function: ``roots(...)``

.. raw:: latex

\vspace{0.5cm} \hrule \begin{center}END\end{center} </rst>