Assignment 5 - 2010 - Solution

Due date(s):	17 November 2010
(PDF)	Assignment questions
(PDF)	Solutions, prepared by Ali and Kevin
Other instructions	Hand-in at class.

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

Question 1 [1.5]

=====

Richardson's extrapolation allows us to accelerate convergence. The technique was used to decrease the error for the forward difference approximation for numerical differentiation.

Using Richardson's extrapolation technique, but this time with the central difference approximation -- an :math:`O(h^2)` method -- derive an improved (accelerated accuracy) approximation. What is the order of this approximation?

Solution

---

The central difference approximation was derived from the Taylor Series expansion:

.. math::

f'(x_0+h) &= f(x_0) + f'(x_0)h + \frac{f(x_0)h^2}{2} + \frac{f(x_0)h^3}{6} + \frac{f'(x_0)h^4}{24} + O(h^5) \\ f'(x_0-h) &= f(x_0) - f'(x_0)h + \frac{f(x_0)h^2}{2} - \frac{f(x_0)h^3}{6} + \frac{f'(x_0)h^4}{24} + O(h^5)

Subtracting one equation from the other, and simplifying:

.. math::

f'(x_0) &= \frac{f(x_0+h) - f(x_0-h)}{2h} + \frac{f(x_0)}{12}h^2 + O(h^4) \\

So the central difference approximate is an O(h^2) approximation, but this series is only a function of the *even* powers of :math:`h`. So applying Richardson's extrapolation, we may write:

.. math::

\begin{array}{rll} \text{A} \qquad f'(x_0) &= \displaystyle \frac{f(x_0+h) - f(x_0-h)}{2h} &+ c_2h^2 + O(h^4) \\ &= Q_2(h) &+ c_2h^2 + O(h^4) \\ \text{B} \qquad f'(x_0) &= \displaystyle \frac{f(x_0+h/2) - f(x_0-h/2)}{h} &+ \displaystyle c_2\frac{h^2}{4} + O(h^4) \\ &= Q_2(h/2) &+ \displaystyle c_2\frac{h^2}{4} + O(h^4) \end{array}

Subtracting: :math:`4\text{B} - \text{A}`:

.. math::

3f'(x_0) &= 4Q_2(h/2) - Q_2(h) + O(h^4) \\ f'(x_0) &= \displaystyle \frac{4Q_2(h/2) - Q_2(h)}{3} + O(h^4) \\ f'(x_0) &= Q_4(h) + O(h^4)

So this improved (accelerated) approximation for :math:`f'(x)` will be of the order :math:`O(h^4)`.

.. raw:: latex

\newpage

Question 2 [2]

==

.. note:: Adapted from the course textbook, problem 24.11

A colleague has designed a new transdermal patch to deliver insulin through the skin to diabetic patients, eliminating the need for painful injections. She has collected the following data on the mass flux, :math:`f`, of insulin being delivered through the patch (and skin), as a function of time, :math:`t`. [Mass flux is the flow rate through a unit area.]

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

math:`t` [hour] 0 1 2 3 4 5 6 10 14 18 24

---

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

math:`f` [mg/(cm\ :sup:`2`.hr)] 15 14 12 11 9 8 7.2 5 2.5 2 1

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

1. . Provide an estimate, using the trapezoidal rule, of the mass of drug delivered to the patient over a 24 hour period with a patch of 10 cm\ :sup:`2`.

1. . Can you provide an alternate, more accurate estimate, using another technique we learned about? If so, provide the answer; if not, explain why.

Solution

---

Note that not all the points are equally spaced. From :math:`t=0` to :math:`t=6` the step size is :math:`h=1`. From :math:`t=6` to :math:`t=18` the step size is :math:`h=4`. And for the last two points, :math:`h=6`. So, we can apply the composite trapezoidal to each of these three partitions:

.. math::

\begin{array}{llrl} \text{Partition 1}& h=1: & I_1=\int_{0}^{6}{f(t)dt}&\approx \frac{1}{2}(f(0)+2f(1)+2f(2)+2f(3)+2f(4)+2f(5)+f(6))\\ &&&=\frac{1}{2}(15+2(14) +2(12)+2(11)+2(9)+2(8)+7.2) = 65.10\\ \text{Partition 2}& h=4: & I_2=\int_{6}^{18}{f(t)dt}&\approx \frac{4}{2}(f(6)+2f(10)+2f(14)+f(18))\\ &&&=\frac{4}{2}(7.2+2(5)+2(2.5)+2)=48.40\\ \text{Partition 3}& h=6: & I_3=\int_{18}^{24}{f(t)dt}&\approx \frac{6}{2}(f(18)+f(24))\\ &&&=\frac{6}{2}(2+1)=9\\ &&\int_{0}^{24}{f(t)dt} &\approx I_1+I_2+I_3=\textbf{122.50} \end{array}

Note the above result is the mass flux. To obtain the mass, we multiply this by the surface area: :math:`M=122.50\times 10=1225.00`

A more accurate result can be obtained by using Simpson's rule, where possible, instead of the trapezoidal rule. Simpson's 1/3 rule can be applied on every three points if they are equally spaced. According to the given table, we can apply composite Simpson's 1/3 rule from :math:`t=0` to :math:`t=6` and use Simpson's 3/8 rule from :math:`t=6` to :math:`t=18`. For the last two points, we again use the trapezoidal method. Summing the results of the above partitions we get:

.. math::

\begin{array}{llrl} \text{Partition 1}& h=1:& I_1= \int_{0}^{6}{f(t)dt}&\approx\frac{1}{3}(f(0)+4f(1)+2f(2)+4f(3)+2f(4)+4f(5)+f(6))\\ &&&=\frac{1}{3}(15+4(14) +2(12)+4(11)+2(9)+4(8)+7.2)=65.40 \\ \text{Partition 2}& h=4:& I_2=\int_{6}^{18}{f(t)dt}&\approx\frac{3\times 4}{8}(f(6)+3f(10)+3f(14)+f(18))\\ &&&=\frac{3\times 4}{8}(7.2+3(5)+3(2.5)+2)=47.55\\ \text{Partition 3}& h=6: & I_3=\int_{18}^{24}{f(t)dt}&\approx \frac{6}{2}(f(18)+f(24))\\ &&&=\frac{6}{2}(2+1)=9\\ &&\int_{0}^{24}{f(t)dt}&\approx I_1+I_2+I_3=\textbf{121.95} \end{array}

Again, we multiply the above value by 10 cm\ :sup:`2` to get the mass of the drug: :math:`M=121.95\times 10=1219.50`

.. raw:: latex

\newpage

Question 3 [3]

==

.. note:: This question is a little unconventional; but we strongly recommended you struggle with it, without asking for help from TA's or other group members.

The following reaction is taking place in a large, well-mixed reactor: :math:`{\sf N} \rightarrow {\sf B}` with rate expression:

.. math::

r = \displaystyle \frac{k_1 C_{\sf N}}{k_2 + C_{\sf N}} \qquad\qquad\,\, k_1 = 30.0\,\,\text{g/(m$^3$.day)}\qquad\qquad\,\, k_2 = 20.4\,\,\text{g/m$^3$}

1. . Derive a dynamic mass balance for species :math:`{\sf N}`, given that:

* the flow into the large tank is :math:`F^\text{in}(t) = 500` m\ :sup:`3`/day * the inlet stream only contains species :math:`{\sf N}`, approximately constant at 100 g/m\ :sup:`3` * the tank is operated at constant volume, :math:`V = 2000` m\ :sup:`3`

1. . Using techniques learned in this course up to 11 November, calculate the concentration of :math:`{\sf N}` as a function of time over the first 3 days (roughly), as the reactor is starting up. The tank was initially filled to 2000 m\ :sup:`3` with material at a concentration of :math:`C_{\sf N} = 100` g/m\ :sup:`3`.

Show your calculations, source code and plots. Please do not simply attach your source code in an appendix; your code should be worked into your solution, as part of your explanations. Also note that you should **not** use methods such as Euler's methods, or MATLAB/Python functions such as ``scipy.integrate.ode`` or ``ode45``.

Solution

---

The main purpose of this question was for you to struggle with the concept of time-based functions, derivatives and integrals in the expectation that when you see ODEs in part F of the course that they will make sense. Any serious attempt at this question will get a :math:`\beta`, while any solutions that calculate a reasonable plot of concentration over time, with concentration decreasing in time will get an :math:`\alpha`.

1. . The species mass balance for :math:`{\sf N}` is:

* In = :math:`C_{\sf N}^\text{in}(t)F^\text{in}(t)` with units of g/day * Out = :math:`C_{\sf N}(t)F^\text{in}(t)` * Depletion of :math:`{\sf N}` = :math:`- \displaystyle \frac{k_1 C_{\sf N}(t)}{k_2 + C_{\sf N}(t)} V(t)` also with units of g/day * Accumulation = :math:`\displaystyle \frac{V(t)C_{\sf N}(t)}{dt}` with units of g/day

Assumptions made in the above are that the tank is well mixed and that the inflow and outflow are the same, which should be true if operating at constant volume. In the bullet points we had all the terms as a function of time. When forming the dynamic balance below we have removed the time-dependence from any terms which are constant.

.. math::

\frac{d\left[V(t)C_{\sf N}(t)\right]}{dt} = V\frac{dC_{\sf N}(t)}{dt} &= F^\text{in}\left[C_{\sf N}^\text{in} - C_{\sf N}(t) \right] - \frac{k_1 C_{\sf N}(t)}{k_2 + C_{\sf N}(t)} V \\ \frac{dC_{\sf N}(t)}{dt} &= \frac{F^\text{in}}{V}\left[C_{\sf N}^\text{in} - C_{\sf N}(t) \right] - \frac{k_1 C_{\sf N}(t)}{k_2 + C_{\sf N}(t)}

1. . Integrating the dynamic balance between time :math:`t=0` and :math:`t=3` days, and initial concentration :math:`C_{\sf N}(t=0) = 100` g/m\ :sup:`3` can be expressed as:

.. math::

\int_{100}^{C_{\sf N}(3)}{f(C_{\sf N}) dC_{\sf N}} &= \int_{0}^{3} dt \\ \int_{100}^{C_{\sf N}(3)}{ \left[\frac{F^\text{in}}{V}\left[C_{\sf N}^\text{in} - C_{\sf N}(t) \right] - \frac{k_1 C_{\sf N}(t)}{k_2 + C_{\sf N}(t)} \right]^{-1} dC_{\sf N}} &= 3.0

This is not easy to integrate analytically, so the tools we have learned in the course can be used now to approximate the integral on the left hand side starting from :math:`C_{\sf N}(t=0) = 100` g/m\ :sup:`3` up to some unknown concentration, :math:`C_{\sf N}(3)` g/m\ :sup:`3`. We will add to the approximation until the left integral is equal to right integral of 3.0 days. This process seems suited to using the trapezoidal rule, or Simpson's method. We can use panels to approximate the integral, and keep adding panels until our approximation to the integral is roughly 3.0 days.

Either way, we need to decide where to place our divisions to form the integrating panels. Do we add them to the left or to the right of the starting point :math:`C_{\sf N}(t=0) = 100` g/m\ :sup:`3`? To help us decide this, we can evaluate the derivative at time :math:`t=0` days:

.. math::

\frac{C_{\sf N}(t)}{dt} &= \frac{500.0}{2000.0}\left[100 - 100\right] - \frac{30 \times 100}{20.4 + 100.0} \\ &= 0 - 24.9\,\, \text{g/(m$^3$.day)}

This gives two pieces of information. Firstly, that the concentration will decrease by about 25 units per day at the startup of the reactor. Secondly it helps us decide on the number of panels. We can use panels about 5 units wide, giving about 5 divisions per day. i.e. :math:`C_{\sf N} = [100, 95, 90, 85, \ldots ]`. So we will keep adding panels until the approximation is roughly equal to 3.0 days.

Using the trapezoidal approximation:

* From :math:`C_{\sf N}` = 100 to 95 g/m\ :sup:`3`: area = :math:`\frac{(95 - 100)}{2} \cdot \left[f(95)+f(100)\right]` = :math:`\frac{-5}{2} \cdot \left[-0.04013 - -0.04265 \right] = 0.20695`, total area so far = 0.21 days. * From :math:`C_{\sf N}` = 95 to 90 g/m\ :sup:`3`: area = :math:`\frac{(90 - 95)}{2} \cdot \left[f(90)+f(95)\right]` = :math:`\frac{-5}{2} \cdot \left[-0.04265 - -0.04554 \right] = 0.2205`, total area so far = 0.43 days. * From :math:`C_{\sf N}` = 90 to 85 g/m\ :sup:`3`: area = 0.2361, area so far = 0.6636 days. * From :math:`C_{\sf N}` = 85 to 80 g/m\ :sup:`3`: area = 0.2545, area so far = 0.9181 days. * From :math:`C_{\sf N}` = 80 to 75 g/m\ :sup:`3`: area = 0.2765, area so far = 1.195 days (which we expected from our initial derivative) * From :math:`C_{\sf N}` = 75 to 70 g/m\ :sup:`3`: area = 0.3031, area so far = 1.498 days. * From :math:`C_{\sf N}` = 70 to 65 g/m\ :sup:`3`: area = 0.3364, area so far = 1.834 days. * From :math:`C_{\sf N}` = 65 to 60 g/m\ :sup:`3`: area = 0.3793, area so far = 2.213 days. * From :math:`C_{\sf N}` = 60 to 55 g/m\ :sup:`3`: area = 0.4369, area so far = 2.650 days. * From :math:`C_{\sf N}` = 55 to 50 g/m\ :sup:`3`: area = 0.5190, area so far = 3.169 days.

The following code was used:

.. twocolumncode:: :code1: ../che3e4/Assignments/Assignment-5/code/integrate_q3.m :language1: matlab :header1: MATLAB code :code2: ../che3e4/Assignments/Assignment-5/code/integrate_q3.py :language2: python :header2: Python code

So we now have the concentrations in the tank as a function of time, over a period of 3 days. We can plot these as shown below (the calculations above are to the left of the red line). The points beyond the red line were calculated when continuing beyond 3 days.

.. figure:: ../che3e4/Assignments/Assignment-5/images/assign5_q3_plot.png :scale: 55 :align: center

The approximation error will obviously get worse and worse with this crude panel size, so we should start using a smaller panel size as we keep going.

The same approach as above could have been performed using Simpson's 1/3 rule for more accurate estimates.

Question 4 [1.5]

=====

.. note:: From the 2009 final exam, 3 points out of 50; 3 hours.

You are given a set of algorithms for numerical integration. You are also given a set of problems. You are to match every problem (1, 2, 3) with an algorithm (A, B, or C). You may use each algorithm **only once**. Motivate each of your answers **carefully**.

• • Problem**

• • 1**: Determine the integral of a very nonlinear function with high accuracy.

• • 2**: Determine the integral of a computationally expensive, black-box function with reasonable accuracy.

• • 3**: Determine the integral of a cubic model for a pure component's specific heat capacity: e.g. :math:`C_p(T) = a + bT + cT^2 + dT^3`

• • Algorithms**

• • A**: Simpson's 1/3 rule

• • B**: Romberg method

• • C**: 4-point Gauss quadrature

Solution

---

Most of the grade in this question is for your justification. Below is one potential set of links:

• • 1**: Determine the integral of a very nonlinear function with high accuracy: use Romberg-method; because Romberg builds on the trapezoidal method, potentially many function evaluations are made. The Romberg method improves its accuracy by order :math:`O(h^2p)` in every row in the Romberg table. So just 4 rows in the table would have an error of :math:`O(h^8)`. We can guarantee the error, even if the function is nonlinear.

• • 2**: Determine the integral of a computationally expensive, black-box function with reasonable accuracy: use 4-point Gauss quadrature; because Gauss quadrature has low error for very few function evaluations. A 4-point quadrature would be exact for 7th-order polynomials, so we would get very good accuracy with only 4 function evaluations.

• • 3**: Determine the integral of a cubic model for a pure component's specific heat capacity: e.g. :math:`C_p(T) = a + bT + cT^2 + dT^3`: I would use Simpson's 1/3 method because it is exact for cubic polynomials and lower - so we would get an exact answer for this type of polynomial.

Question 5 [2]

===

.. note:: From the 2009 final exam, 2+1.5 points, out of 50; 3 hours. There were other parts to the question, but they are not relevant here.

Consider a CSTR, fed with one inlet stream at volumetric flow rate :math:`F^{\rm in}` and containing species :math:`{\sf A}` at concentration :math:`C_{\sf A}^{\rm in}=1` mol/L. There is no outlet from the reactor. A fluid of constant density is considered throughout this question.

The reactor is initially empty, and the fluid is fed into the tank via the inlet stream at the volumetric flow rate given by :math:`F^{\rm in}(t) = 1-\exp(-t)` [in m\ :sup:`3`/min].

1. . Determine the volume of liquid :math:`V(t)` accumulated at :math:`t_1=5` min using Romberg integration (perform **3 steps** of Romberg's method; that is, 3 columns in the graphical depiction).

1. . Repeat the calculation using two-point Gauss quadrature, and compare the results. Which method is likely to provide the best estimate?

Solution

---

This question asks to evaluate the cumulative volume in the tank after a certain time:

.. math::

\frac{dV}{dt} &= F^{\rm in}(t) \\ \int_{0}^{V}{dV} &= V = \int_{0}^{t}{F(t) dt} = \int_{0}^{t=5}{1-e^{-t} dt}

1. . Using Romberg integration:

The Romberg table can be constructed as follows:

.. math::

\begin{array}{llll} h & I_1, O(h^2) & I_2, O(h^4) & I_3, O(h^6) \\ \hline h=5 & I_1(5) = 2.4832 & \\ h=2.5 & I_1(2.5) = 3.5364 & I_2(2.5) = 3.8874 &\\ h=1.25 & I_1(1.25) = 3.8807 & I_2(1.25) = 3.9955 & I_3(1.25) = 4.0027\\ \hline \end{array}

The following expanded entries were used in the above table:

* :math:`I_1(5) = \frac{5}{2}\left[f(0)+f(5)\right] = 2.4832` * :math:`I_1(2.5) = \frac{5}{2 \times 2}\left[f(0)+2f(2.5) + f(5)\right] = 3.5364` * :math:`I_1(1.25) = \frac{5}{4 \times 2}\left[f(0)+ 2f(1.25) + 2f(2.5) + 2f(3.75) + f(5)\right] = 3.8807` * :math:`I_2(2.5) = \frac{4I_1(2.5) - I_1(5)}{4-1} = \frac{4 \times 3.5364 - 2.4832}{3} = 3.8875` using the messy formula with :math:`j=1` and :math:`k=1` * :math:`I_2(1.25) = \frac{4I_1(1.25) - I_1(2.5)}{4-1} = \frac{4 \times 3.8807 - 3.5364}{3} = 3.9955` with :math:`j=1` and :math:`k=2` * :math:`I_3(1.25) = \frac{16I_2(1.25) - I_2(2.5)}{16-1} = \frac{16 \times 3.9955 - 3.8875}{15} = 4.0027` with :math:`j=2` and :math:`k=2`

The final estimate for the integral is the last column in the last row, so :math:`I_3(1.25) = 4.003`

1. . Using Gauss quadrature: our function to be integrated is :math:`f(t) = 1-e^{-t}` between limits :math:`a=0` and :math:`b=5`. Since Gaussian quadrature is integrated between -1 and +1, our formula for the change in variables is:

	.. math::

t &= \frac{b+a}{2} + \frac{b-a}{2}\xi \\ t &= \frac{5}{2} + \frac{5}{2}\xi~~\text{where}\,\,-1 \leq \xi \leq +1

.. math::

I_2 &= \omega_1 f(\xi_1) + \omega_2 f(\xi_2)\\ &= \bigg[\omega_1 f(t_1) + \omega_2 f(t_2)\bigg] \frac{dt}{d\xi}\\ &= \bigg[1.0 f\left(\frac{5}{2} - \frac{5}{2}\cdot\frac{1}{\sqrt{3}}\right) + 1.0 f\left(\frac{5}{2} + \frac{5}{2}\cdot\frac{1}{\sqrt{3}}\right)\bigg] \frac{5}{2}\\ &= \bigg[1.0 f(1.0566) + 1.0 f(3.9434)\bigg] \frac{5}{2}\\ &= \bigg[0.6524 + 0.9806 \bigg] \frac{5}{2}\\ &= 4.083

Romberg quadrature taken to three rows will have error in the order of :math:`O(h^6)`, while 2-point Gaussian quadrature will be exact for polynomials of order 3 and below: i.e. the error of Gaussian quadrature will be some function of the fourth derivative of :math:`f(t)`. Errors from Gaussian quadrature are not a function of :math:`h`, since there is no :math:`h` in Gaussian quadrature. Since the function we are integrating is not polynomial, I would expect error in both methods.

The true value of the integral is easy to find: :math:`\int_{0}^{t=5}{1-e^{-t} dt} = t + e^{-t}\big|_{0}^{5} = 4 + e^{-5} = 4.0067379`; using this we can see the Romberg method had lower error.

.. raw:: latex

\vspace{0.5cm} \hrule \begin{center}END\end{center}

</rst>