TMA681 Tillämpad matematik K, 2001

STUDIO 3. THE TANK REACTOR: ARRHENIUS’ LAW.

Important: write a readable report of your work in each studio session. You will need this at the end of the course when you prepare for the written exam; several exam questions will be based directly on the studio work which is only documented in your own notes and the instructions.

Some exercises are called homework in order to save time in the studio classroom.

1. The method of least squares

(TLA 4.9) Consider the linear system of equations

(1)  Ax = b,

where A ∈ R^{m×n}, x ∈ R^{n}, b ∈ R^{m}. If m > n (more equations than unknowns), then the system is “overdetermined” and such a system has no solution in general. Geometrically, the reason for this is that (in general) the vector b lies outside the range (“värderummet”) of A,

R(A) = {y ∈ R^{m} : y = Ax for some x ∈ R^{n}}.

Since the dimension of R(A) is ≤ n and m > n, we realize that the space R(A) does not “fill out” the whole space R^{m}. Therefore it is likely that a given vector b ∈ R^{m} will lie outside R(A), see Figure 1, and then Ax cannot be equal to b and (1) has no solution.

PSfrag replacements

b

b − Ax̂

Ax̂ Av R(A)

Figure 1. Orthogonal projection onto the range of A.

In this situation we seek an approximate solution which makes the residual

b − Ax

as small as possible. More precisely, we seek a vector x̂ ∈ R^{n}, which minimizes the square of the norm (length) of the residual:

(2)  f(x̂) = min f(x), f(x) = ∥b − Ax∥^{2}.

Recall the scalar product ⟨x, y⟩ = y^{t}x and the norm ∥x∥ = \sqrt{⟨x, x⟩} = \sqrt{x^{t}x} of column vectors. From (TLA 6.9) we know that there is a unique vector ŷ = Ax̂ ∈ R(A) such that the distance ∥b − Ax∥ is minimal, i.e., ∥b − Ax̂∥ ≤ ∥b − Ax∥ for all x. The vector ŷ = Ax̂ is the orthogonal projection of b onto R(A), see Figure 1. It is characterized by the condition that b − ŷ = b − Ax̂ is orthogonal to all vectors Av ∈ R(A). This means that

Date: March 28, 2001, Niklas Ericsson, Stig Larsson, Mathematics, Chalmers University of Technology.

0 = ⟨b − Ax̂, Av⟩ = (Av)^{t}(b − Ax̂) = v^{t}A^{t}(b − Ax̂) = v^{t}(A^{t}b − A^{t}Ax̂).

Since this holds for all v ∈ R^{n}, we may take v = A^{t}b − A^{t}Ax̂ to get

(3)  A^{t}Ax̂ − A^{t}b = 0.

Thus, we can compute x̂ by solving the linear system

(4)  A^{t}Ax = A^{t}b.

Note that the coefficient matrix A^{t}A is n × n and symmetric. The system (4) has at least one solution (namely x̂).

In order to see that the minimization problem (2) is equivalent to solving the linear system (4), we write x = x̂ + v, v = x − x̂, and compute

\begin{aligned}f(x) & = f(x̂ + v) = ∥(b − Ax̂) − Av∥^{2} = ⟨(b − Ax̂) − Av, (b − Ax̂) − Av⟩ \\ & = ⟨b − Ax̂, b − Ax̂⟩ − 2⟨b − Ax̂, Av⟩ + ⟨Av, Av⟩ \\ & = ∥b − Ax̂∥^{2} + 2(Av)^{t}(Ax̂ − b) + ∥Av∥^{2} \\ & = f(x̂) + 2v^{t}A^{t}(Ax̂ − b) + ∥Av∥^{2} \\ & = f(x̂) + 2v^{t}(A^{t}Ax̂ − A^{t}b) + ∥Av∥^{2}. \end{aligned} ; (5)

Taking (3) into account we get

f(x) = f(x̂) + ∥Av∥^{2} ≥ f(x̂).

This shows that x minimizes f(x), if and only if ∥Av∥ = 0, in which case Ax = Ax̂ + Av = Ax̂ and A^{t}Ax = A^{t}Ax̂ = A^{t}b. Therefore, x minimizes f(x) if and only if x is a solution of (4).

We can also interpret this in terms of the general minimization problem in (MA3 2.9–10). Recall Taylor’s formula:

(6)  f(x) = f(x̂ + v) = f(x̂) + v^{t}f′(x̂) + \frac{1}{2}v^{t}f′′(x̂)v + R(x).

Noting that ∥Av∥^{2} = ⟨Av, Av⟩ = (Av)^{t}(Av) = v^{t}(A^{t}A)v, we re-write (5) as

(7)  f(x) = f(x̂ + v) = f(x̂) + 2v^{t}(A^{t}Ax̂ − A^{t}b) + v^{t}(A^{t}A)v.

Comparing (6) with (7), we identify the Jacobi matrix (gradient vector) f′(x̂) = 2(A^{t}Ax̂ − A^{t}b), the Hesse matrix f′′(x̂) = 2A^{t}A, and the remainder R(x) = 0. Recall that stationary points are given by the system of equations f′(x) = 2(A^{t}Ax − A^{t}b) = 0, which is the same as (4). Note also that the Hesse matrix is constant (with respect to x) and positive semidefinite: v^{t}f′′(x)v = 2v^{t}(A^{t}A)v = 2∥Av∥^{2} ≥ 0. Exercise 1. Suppose that the variables y and x are related by y = kx + m. In order to determine the coefficients k and m we make measurements of y and x:

\begin{matrix}x & 5 & 6 & 7 & 8 & 9 & 10 \\ y & 19.5888 & 23.4043 & 25.5754 & 29.1231 & 31.9575 & 35.8116 \end{matrix}

This leads to an overdetermined system of the form

kx_{1} + m = y_{1}
⋮
kx_{6} + m = y_{6}

or, in matrix form Av = y,

[\begin{matrix}x_{1} & 1 \\ ⋮ & ⋮ \\ x_{6} & 1 \end{matrix}] [\begin{matrix}k \\ m \end{matrix}] = [y_{1} ; ⋮ ; y_{6}] .

Solve this system by the least squares method in Matlab. Hint: set up the column vectors x, y and the matrix A=[x ones(size(x))], then form the matrices B=A'*A and g=A'*y. Solve the system Bv = g by the command v=B\g.

Plot the data points (x_{i}, y_{i}) and the fitted function y = kx+m in the same figure. The following

commands are useful: plot(x,y,'or'), fplot('ykxm',[x(1) x(6)]). Here ykxm.m is a function

file that implements the function y = kx + m. Don’t forget to declare global k m both inside the function file and in the main program.

Actually, Matlab’s backslash command v=A\y automatically uses the least squares method when the system Av = y is overdetermined. Try this also!

2. The tank reactor

The rate coefficient depends on the temperature according to the Arrhenius law:

k = k_{0} exp(−E/(RT)) [s^{−1}] (8)

where R [8.31 J/(mol K)] is the gas constant, E [J/mol] is the activation energy and k_{0} [s^{−1}] is the rate constant of the reaction. The following rates have been measured:

\begin{matrix}T [K] & 343 & 353 & 363 & 373 & 383 & 393 & 403 \\ k [s^{−1}] & 2.8 · 10^{−5} & 5.6 · 10^{−5} & 11.2 · 10^{−5} & 22.4 · 10^{−5} & 44.8 · 10^{−5} & 89.6 · 10^{−5} & 179.2 · 10^{−5} \end{matrix}

The task is now to determine the coefficients k_{0} and E by fitting the rate law (8) to these data. Last week we wrote (8) in dimensionless form

(9)  kτ = δe^{γ(1−1/X)}, where γ = \frac{E}{RT_{f}} , δ = k_{0}τe^{−γ}, X = \frac{T}{T_{f}} , τ = \frac{V}{q_{ref}} .

Introducing new variables r = kτ and ξ = 1 − 1/X we get

(10)  r = δe^{γξ}.

The task is now to fit this function to the given data points (ξ_{i}, r_{i}). Exercise 2. (Linear least squares method.) Form the logarithm of (10) so that you get a linear relation of the form y = kx + m, namely, (see also (TLA 5.2))

(11)  log(r) = γξ + log(δ).

(Note that the natural logarithm is denoted log(x) in English and in Matlab, but ln(x) in Swedish.) Solve for γ and δ by using the least squares method as in Exercise 1. Begin by forming column vectors X, r, xi and so on. Plot the data points (X_{i}, r_{i}) and the fitted function r = δe^{γ(1−1/X)} in the same figure. Finally, determine k_{0} and E. Homework 1. (A nonlinear least squares method.) Alternatively, we can form the residual r − δe^{γξ} from the nonlinear relation (10) and minimize the square of its norm

(12)  g(δ, γ) = ∑_{i} (r_{i} − δe^{γξ_{i}})^{2}.

Write a Matlab function that implements this function and use Matlab’s program fminsearch to minimize it. Why does this method give a slightly different result? Hint: the Matlab function norm may be useful for computing the right side of (12). Exercise 3. Insert the new values for δ and γ in your Matlab programs from Studio 2. Repeat all the computations. Let U_{1} and U_{2} be equal to constant values Ū_{1} and Ū_{2}. Does the solution

X(s) = [\begin{matrix}X_{1}(s) \\ X_{2}(s) \end{matrix}] approach an equilibrium X̄ = [X̄_{1} ; X̄_{2}] as s → ∞? In fact, you should be able to find

two equilibrium points by choosing different initial values X_{0}, say, X_{0} = [\begin{matrix}0.5 \\ 1 \end{matrix}] and X_{0} = [\begin{matrix}0.5 \\ 1.1 \end{matrix}].

Next week we will look for an equilibrium at X̄_{1} = 0.5 and analyze the stability of this desired operating point.