Tutorial 7 - 2010 - Solution
| Due date(s): | 11 November 2010 |
 (PDF)
|
Tutorial questions |
| (PDF)
|
Solutions by Kevin |
| Other instructions | Hand-in at class. |
<rst>
<rst-options: 'toc' = False/>
<rst-options: 'reset-figures' = False/>
.. |m3| replace:: m\ :sup:`3`
.. highlight:: python
.. rubric:: Tutorial objectives: Numerical differentiation and plotting.
Question 1 [2]
====
This question uses the experimental data from the class notes. A batch system was charged with reactant :math:`{\sf A}` at time :math:`t=0` at a concentration of 1.0 mol/L. As the assumed first-order reaction evolved, the concentrations were recorded:
.. csv-table:: :header: "Time (minutes)", "Concentration (mol/L)"
0, 1.00 5, 0.84 10, 0.72 20, 0.57 30, 0.47 45, 0.37 60, 0.30
The dynamic balance for this system can be shown to be :math:`\displaystyle \frac{dC_{\sf A}(t)}{dt} = -k C_{\sf A}(t)`.
- . Reproduce the above table and add a column that provides an approximation of :math:`\displaystyle \frac{dC_{\sf A}(t)}{dt}` at each time step.
- . Add a fourth column that lists the order of your approximation listed in part 1.
- . Use a 2nd order Lagrange polynomial to approximate the :math:`C_{\sf A}(t)` data, and use this polynomial to estimate :math:`C_{\sf A}\,'(t)` at time :math:`t=15` minutes.
- . Can you provide an :math:`O(h^2)` estimate of :math:`\displaystyle \frac{dC_{\sf A}(t)}{dt}` at :math:`t=0`? If so, please calculate it.
Solution
.. NOTE: For code that calculates these answers, please see: Notes/Figures/E-Diff-Int/plot_concentrations.py
A table showing the derivative approximation, and the order:
.. csv-table:: :header: "Time (minutes)", "Concentration (mol/L)", Type of approximation, Approximation, Order
0, 1.00, Forward, :math:`\frac{0.84 - 1.00}{5 - 0} = -0.032`, :math:`O(h)` 5, 0.84, Central, :math:`\frac{0.72 - 1.00}{2(10 - 5)} = -0.028`, :math:`O(h^2)` 10, 0.72, Central, :math:`\frac{0.57 - 1.00}{2(20 - 10)} = -0.0215`, :math:`O(h^2)` 20, 0.57, Central, :math:`\frac{0.47 - 0.72}{2(30 - 20)} = -0.0125`, :math:`O(h^2)` 30, 0.47, Backward, :math:`\frac{0.47 - 0.57}{30 - 20} = -0.01` , :math:`O(h)` 45, 0.37, Central, :math:`\frac{0.30 - 0.47}{2(60 - 45)} = -0.005667`, :math:`O(h^2)` 60, 0.30, Backward, :math:`\frac{0.30 - 0.37}{60 - 45)} = -0.004667`, :math:`O(h)`
Note that the table uses an :math:`O(h^2)` estimate where possible; you should always prefer an :math:`O(h^2)` estimate over an :math:`O(h)` estimate, if it is possible to calculate one.
In general, one should fit a polynomial using the points closest to where you want to interpolate or use that polynomial. [There are exceptions of course: if one of the points was clearly noisy, relative to the others, you can skip that point and use another instead.] In this case one should fit the polynomial using points :math:`t=(5, 10, 20)`, and their corresponding concentrations, :math:`C_{\sf A} = (0.84, 0.72, 0.57)`. Although in this case it didn't change the answer by much using points at :math:`t=(0, 10, 20)`.
.. math::
P_2(x) &= y_{1}\ell_{1}(x) + y_{2}\ell_{2}(x) + y_{3}\ell_{3}(x) \\ P'_2(x) &= y_{1}\ell'_{1}(x) + y_{2}\ell'_{2}(x) + y_{3}\ell'_{3}(x) \\ &= 0.84 \cdot \frac{2x - 10 - 20}{( 5-10)( 5-20)} + 0.72 \cdot \frac{2x - 5 - 20}{(10- 5)(10-20)} + 0.57 \cdot \frac{2x - 5 - 10}{(20- 5)(20-10)}\\ P'_2(15) &= 0.84 \cdot 0 + 0.72 \cdot \frac{5}{-50} + 0.57 \cdot \frac{15}{150} = -0.15 The polynomial, superimposed on the data points is shown below:
.. figure:: ../che3e4/Tutorials/Tutorial-7/images/tut7_q1_plot.png :scale: 50 :align: center
An :math:`O(h^2)` estimate can be found at :math:`t=0` when using Richardson's extrapolation with :math:`h=10` (full step), and :math:`h=5` (half step).
* :math:`h = 10`: :math:`Q_1(h) = \frac{0.72 - 1.00}{10 - 0} = -0.028` * :math:`h = 5`: :math:`Q_1(h/2) = \frac{0.84 - 1.00}{5 - 0} = -0.032` * Combine these to calculate a :math:`Q_2(h) = 2Q_1(h/2) - Q_1(h) = 2(-0.032) - (-0.028) = -0.036`
Question 2 [1.5]
====
A second order reaction in a constant-volume batch reactor has the following time-dependent behaviour:
.. math::
C_{\sf A}(t) = \frac{C_{\sf A}(0)}{C_{\sf A}(0)kt + 1} \qquad\,\text{with}\qquad\, k = 0.002 \frac{\text{m}^3}{\text{mol.s}} \qquad\,\text{and}\qquad\,C_{\sf A}(0) = 2.5 \frac{\text{mol}}{\text{m}^3}
Starting with :math:`h = 1.0` and decreasing in powers of 10, which step size, :math:`h=\Delta t`, provides the lowest error estimate of :math:`\displaystyle \left.\frac{dC_{\sf A}}{dt}\right|_{t=100}` at :math:`t=100` seconds?
Present your results in tabular form with these 5 columns:
* Step size, :math:`h = \Delta t` * Forward difference approximation * Error in the forward difference approximation * Central difference approximation * Error in the central difference approximation
Solution
.. note:: the last column was intended to be the error for the *central difference* approximation; an earlier version of the tutorial asked for the *backward difference* error. Either values will be accepted for the solution.
Given the above equation for :math:`C_{\sf A}(t)`, the analytical derivative can be found as:
.. math::
\frac{dC_{\sf A}(t)}{dt} = -k\frac{C^2_{\sf A}(0)}{\left(C_{\sf A}(0)kt + 1\right)^2} = -k C^2_{\sf A}(t)
which is expected for a second-order reaction system. So the true derivative at :math:`t=100` seconds is:
.. math::
\frac{dC_{\sf A}(100)}{dt} = -0.002\frac{2.5^2}{(2.5^2 \times 0.002 \times 100 + 1)^2} = -0.005556
The required table is shown below. Two extra columns for the relative errors have been added, since that is arguably more useful than the absolute error (which was asked for in the question).
.. math::
\begin{array}{l|lll|lll} \hline \hline h &\text{Forward diff} & \text{Fwd diff error} & \text{Fwd relative error[%]} & \text{Central diff} & \text{Central diff error} & \text{Central relative error[%]} \\ \hline \hline 1 & -0.0055371 & 1.8457\times 10^{-5} & -0.33223 & -0.0055556 & 6.1729\times 10^{-8} & -0.0011111\\ 0.1 & -0.0055537 & 1.8512\times 10^{-6} & -0.033322 & -0.0055556 & 6.1728\times 10^{-10} & -1.1111\times 10^{-5}\\ 0.01 & -0.0055554 & 1.8518\times 10^{-7} & -0.0033332 & -0.0055556 & 6.1648\times 10^{-12} & -1.1097\times 10^{-7}\\ 1\times 10^{-3} & -0.0055555 & 1.8518\times 10^{-8} & -0.00033333 & -0.0055556 & 1.4742\times 10^{-13} & {\bf -2.6535\times 10^{-9}}\\ 1\times 10^{-4} & -0.0055556 & 1.8496\times 10^{-9} & -3.3293\times 10^{-5} & -0.0055556 & 1.1466\times 10^{-12} & -2.0639\times 10^{-8}\\ 1\times 10^{-5} & -0.0055556 & 1.6206\times 10^{-10}& -2.917\times 10^{-6} & -0.0055556 & 1.558\times 10^{-11} & -2.8043\times 10^{-7}\\ 1\times 10^{-6} & -0.0055556 & 1.558\times 10^{-11} & {\bf-2.8043\times 10^{-7}} & -0.0055556 & 1.558\times 10^{-11} & -2.8043\times 10^{-7}\\ 1\times 10^{-7} & -0.0055556 & 4.5967\times 10^{-10}& -8.274\times 10^{-6} & -0.0055556 & 4.5967\times 10^{-10} & -8.274\times 10^{-6}\\ 1\times 10^{-8} & -0.0055556 & 4.5967\times 10^{-10}& -8.274\times 10^{-6} & -0.0055556 & 4.5967\times 10^{-10} & -8.274\times 10^{-6}\\ 1\times 10^{-9} & -0.0055556 & 4.5967\times 10^{-10}& -8.274\times 10^{-6} & -0.0055556 & 4.5967\times 10^{-10} & -8.274\times 10^{-6}\\ 1\times 10^{-10} & -0.0055578 & 2.2209\times 10^{-6} & -0.039976 & -0.0055567 & 1.1107\times 10^{-6} & -0.019992\\ 1\times 10^{-11} & -0.0055733 & 1.7764\times 10^{-5} & -0.31975 & -0.0055622 & 6.6618\times 10^{-6} & -0.11991\\ 1\times 10^{-12} & -0.0055511 & 4.4404\times 10^{-6} & -0.079928 & -0.0054401 & 0.00011546 & -2.0783\\ \hline \hline \end{array}
What's remarkable in the above table is how accurate the central different approximation is, even using very large values of :math:`h`. The point here is that the central difference should be used instead of the forward or backward difference, wherever possible.
Furthermore, when using the central difference approximation with very small step sizes we start to accumulate round-off error, but no faster than the error we would have accumulated with the forward difference approximation.
Bonus question [0.5]
=========
Show the two error columns from question 2 in graphical form, using log-axes for both the x- and y-axis. Search for help over the internet if you don't know how to plot log-plots in MATLAB, or with Python's ``matplotlib`` library. Add a legend, grid lines and axis labels to the plot.
Solution
.. note:: A clear plot, showing the points, lines connecting the points (different colours or dashed/solid lines), labels, legend on a log-log axis are required to obtain full grade for this question.
A log-log plot with the required labels, for the absolute errors is shown here. The code to generate the plot is below.
.. raw:: latex
\newpage
.. figure:: ../che3e4/Tutorials/Tutorial-7/images/tut7_qbonus_errors.png :scale: 60 :align: center
.. twocolumncode:: :code1: ../che3e4/Tutorials/Tutorial-7/code/table_of_errors.m :language1: matlab :header1: MATLAB code :code2: ../che3e4/Tutorials/Tutorial-7/code/table_of_errors.py :language2: python :header2: Python code
</rst>