Assignment 4 - 2010

Due date(s):	08 November 2010
(PDF)	Assignment questions
Other instructions	Hand-in at class.

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

.. rubric:: Questions 1, 2 and 3 will be posted at the end of the semester in assignment 4(a).

.. rubric:: Solutions by Kevin and Ali

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?

Solution

---

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

1. . Recall that an :math:`n-1` order Lagrange interpolating polynomial, for a set of :math:`n` data points, is given by the following relation:

.. math::

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}}

For our system the quadratic Lagrange interpolating polynomial is the following:

.. math::

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)} + 5.37 \cdot \frac{(x-20)(x-80)}{(60-20)(60-80)} + 17.4 \cdot \frac{(x-20)(x-60)}{(80-20)(80-60)}\\

1. . 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::

\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}

So the Newton interpolating polynomial is:

.. math::

P_{2}(x) = 1.40 + 0.09925(x-x_{1}) + 0.008371(x-x_{1})(x-x_{2})

1. . 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.

.. 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 \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::

\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_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}`: .. math:: {\bf Xa} &= {\bf y} \\ \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\{ \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 :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 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. #. The question was a bit unclear; the intention was to compare all estimates of :math:`\nu(p=40)`: * 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. * Newton polynomial estimate = :math:`1.40 + 0.09925(40-20) + 0.008371(40-20)(40-60) = 0.0366` millipascals. * 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. * Linear interpolation = 3.385 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] ======'"`UNIQ--h-3--QINU`"'========= The following data are collected from a bioreactor experiment, during the growth phase. ====='"`UNIQ--h-4--QINU`"'============= ===== ===== ===== ===== ===== 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 ====='"`UNIQ--h-5--QINU`"'============= ===== ===== ===== ===== ===== 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. 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]

========

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

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

</rst>