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

# 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 $$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:

$\begin{split}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)\end{split}$

Subtracting one equation from the other, and simplifying:

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

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

$\begin{split}\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}\end{split}$

Subtracting: $$4\text{B} - \text{A}$$:

$\begin{split}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)\end{split}$

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

# 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, $$f$$, of insulin being delivered through the patch (and skin), as a function of time, $$t$$. [Mass flux is the flow rate through a unit area.]

 $$t$$ [hour] 0 1 2 3 4 5 6 10 14 18 24 $$f$$ [mg/(cm2.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 cm2.
2. 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 $$t=0$$ to $$t=6$$ the step size is $$h=1$$. From $$t=6$$ to $$t=18$$ the step size is $$h=4$$. And for the last two points, $$h=6$$. So, we can apply the composite trapezoidal to each of these three partitions:

$\begin{split}\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}\end{split}$

Note the above result is the mass flux. To obtain the mass, we multiply this by the surface area: $$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 $$t=0$$ to $$t=6$$ and use Simpson's 3/8 rule from $$t=6$$ to $$t=18$$. For the last two points, we again use the trapezoidal method. Summing the results of the above partitions we get:

$\begin{split}\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}\end{split}$

Again, we multiply the above value by 10 cm2 to get the mass of the drug: $$M=121.95\times 10=1219.50$$

# 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: $${\sf N} \rightarrow {\sf B}$$ with rate expression:

$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 $${\sf N}$$, given that:

• the flow into the large tank is $$F^\text{in}(t) = 500$$ m3/day
• the inlet stream only contains species $${\sf N}$$, approximately constant at 100 g/m3
• the tank is operated at constant volume, $$V = 2000$$ m3
2. Using techniques learned in this course up to 11 November, calculate the concentration of $${\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 m3 with material at a concentration of $$C_{\sf N} = 100$$ g/m3.

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 $$\beta$$, while any solutions that calculate a reasonable plot of concentration over time, with concentration decreasing in time will get an $$\alpha$$.

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

• In = $$C_{\sf N}^\text{in}(t)F^\text{in}(t)$$ with units of g/day
• Out = $$C_{\sf N}(t)F^\text{in}(t)$$
• Depletion of $${\sf N}$$ = $$- \displaystyle \frac{k_1 C_{\sf N}(t)}{k_2 + C_{\sf N}(t)} V(t)$$ also with units of g/day
• Accumulation = $$\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.

$\begin{split}\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)}\end{split}$
2. Integrating the dynamic balance between time $$t=0$$ and $$t=3$$ days, and initial concentration $$C_{\sf N}(t=0) = 100$$ g/m3 can be expressed as:

$\begin{split}\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\end{split}$

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 $$C_{\sf N}(t=0) = 100$$ g/m3 up to some unknown concentration, $$C_{\sf N}(3)$$ g/m3. 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 $$C_{\sf N}(t=0) = 100$$ g/m3? To help us decide this, we can evaluate the derivative at time $$t=0$$ days:

$\begin{split}\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)}\end{split}$

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

The following code was used:

MATLAB code Python code
% The function on the right hand side of the integral: i.e. the
% mass balance.  Please note: this is not a good way to write ODE
% functions in MATLAB.  I'm doing it here for convenience, so that
% the entire solution is in one m-file.
dCN_dt = @(t, CN) 500/2000*(100-CN) - 30*CN/(20.4 + CN);

%  The function we will integrate on the left hand side is just the
%  inverse of the derivative
func = @(CN) 1./dCN_dt(0, CN);

% Step down in steps of -5 concentration units
concs = 100.0:-5:25.0;
integral = 0.0;

time = zeros(size(concs));
x1 = concs(1);
k=2;
for x2 = concs(2:end)
panel_area = (x2-x1)/2.0 * (func(x2) + func(x1));
integral = integral + panel_area;
time(k) = integral;
%print k, panel_area, integral, x1, x2
x1 = x2;
k = k+1;
end

plot(time, concs, 'ko-')
grid('on')
xlabel('Time [days]')
ylabel('Concentration of N [g/m^3]')
axis([-0.1, 15.0, 25.0, 104.0])

% Use a proper integrator (should not be used in this assignment:
% we will use these techniques in assignment 6 only).
CN_t_zero = 100.0;
t_start = 0.0;
t_final = 3.0;

[t, CN] = ode45(dCN_dt, [t_start, t_final], [CN_t_zero]);

import numpy as np
from scipy import integrate
from matplotlib.pylab import *

def dCN_dt(t, CN):
""" The derivative function: dCN/dt = ...
INPUTS:

t:  time (this particular function is not dependent on time)
CN: the value of C_N(t)
"""
V = 2000.0    # m^3
F = 500.0     # m^3/day
k1 = 30.0     # g/(m^3.day)
k2 = 20.4     # g/m^3
Cin = 100.0   # g/m*3
return F/V*(Cin-CN) - k1*CN/(k2 + CN)

def func(CN):
""" The function we will integrate on the left hand side."""
return 1/dCN_dt(0, CN)

concs = np.arange(100.0, 25.0, -5.0)
integral = 0.0

time = concs.copy()
time[0] = 0.0
x1 = concs[0]

for k, x2 in enumerate(concs[1:]):
panel_area = (x2-x1)/2.0 * (func(x2) + func(x1))
integral += panel_area
time[k+1] = integral
print k, panel_area, integral, x1, x2
x1 = x2

fig = figure()
plot(time, concs, 'ko-', ms=10)
grid('on')
xlabel('Time [days]')
ylabel('Concentration of N [g/m^3]')
axis([-0.1, 15.0, 25.0, 104.0])
axvline(x=3, color='r')
fig.savefig('../images/assign5_q3_plot.png')

# Use a proper integrator (should not be used in this assignment:
# we will use these techniques in assignment 6 only).
r = integrate.ode(dCN_dt).set_integrator('vode', method='bdf')
CN_0 = 100.0
t_0 = 0.0
r.set_initial_value(CN_0, t_0)
t1 = 3.0
dt = 0.10
while r.successful() and r.t < t1:
r.integrate(r.t+dt)
print (r.t, r.y)


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.

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. $$C_p(T) = a + bT + cT^2 + dT^3$$

Algorithms

A: Simpson's 1/3 rule

B: Romberg method

## 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 $$O(h^2p)$$ in every row in the Romberg table. So just 4 rows in the table would have an error of $$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. $$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 $$F^{\rm in}$$ and containing species $${\sf A}$$ at concentration $$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 $$F^{\rm in}(t) = 1-\exp(-t)$$ [in m3/min].

1. Determine the volume of liquid $$V(t)$$ accumulated at $$t_1=5$$ min using Romberg integration (perform 3 steps of Romberg's method; that is, 3 columns in the graphical depiction).
2. 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:

$\begin{split}\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}\end{split}$
1. Using Romberg integration:

The Romberg table can be constructed as follows:

$\begin{split}\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}\end{split}$

The following expanded entries were used in the above table:

• $$I_1(5) = \frac{5}{2}\left[f(0)+f(5)\right] = 2.4832$$
• $$I_1(2.5) = \frac{5}{2 \times 2}\left[f(0)+2f(2.5) + f(5)\right] = 3.5364$$
• $$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$$
• $$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 $$j=1$$ and $$k=1$$
• $$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 $$j=1$$ and $$k=2$$
• $$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 $$j=2$$ and $$k=2$$

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

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

$\begin{split}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\end{split}$
$\begin{split}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\end{split}$

Romberg quadrature taken to three rows will have error in the order of $$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 $$f(t)$$. Errors from Gaussian quadrature are not a function of $$h$$, since there is no $$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: $$\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.