29  Final

Time Allowed: 2 hours

Answer SECTION A (40%) and TWO (30% + 30%) of the three optional sections (B, C, D). If you have answered more than than the required two optional sections, you will be given credit for your best two optional sections. As a guide, I suggest you spend around 40 minutes on each section.

Please indicate your answers to multiple choice questions clearly with an ╳ to the left of your answers. Answer the optional sections on separate paper.

Please label any additional paper you use with your name and page number.

Electronic devices are not needed and not permitted in this examination. You may use your own notes, limited to 2 sides of letter sized paper. Answer the questions to Section A in the boxes provided. Additional paper may be requested.

29.1 A: Compulsory Section

  1. [2 points] Suppose that the speed of light in a vacuum is c = 2.9979 \times 10^8 meters per second (this is exact to 5 significant digits). In Physics, c is often approximated by 3 \times 10^8 meters per second. Write down a formula for the relative error in this approximation. Don’t simplify your expression.
  1. [2 points] Suppose that \xi = g(\xi) is a fixed point of a twice continuously differentiable function g with g'(\xi) = 0. Further suppose that x_{n+1} = g(x_n) converges to \xi. What is the order of convergence?

Section A continued.

  1. [2 points] Note that x := 1.001^{2025} and y := e^{2025 \log 1.001} are mathematically identical (x=y).
    In the following code, we use 16 bit floating point arithmetic to compute x (by repeated multiplication) and y (with the above formula):
b, n, x = Float16( 1.001 ), 2025, Float16( 1 );
for i = 1:n
    x = b*x
end
y = exp( n * log( b ) ) 
println( "x = ", x, "; y = ", y )
x = 5.91; y = 7.22

Which result do you trust more and why?

  1. [5 points]
    (i) Show that f(x) = \cos(x) - 2 x has a root \xi \in [0,\frac{\pi}{2}],
    (ii) Write the root as a fixed point of some contraction g:[0,\frac{\pi}{2}] \to [0, \frac\pi2],
    (iii) Does the iteration x_{n+1} := g(x_n) converge for all x_1 \in [0,\frac\pi2]?

Section A continued.

  1. [2 points] True or false: there exists a continuous function f with a root in [a,b] and f(a) f(b) > 0? Why?
  1. [2 points] True or false: all continuous functions f: [a,b] \to \mathbb R have fixed points \xi = f(\xi)? Why?
  1. [5 points] For a fixed constant \alpha > 0, consider the iteration x_{n+1} = 2x_n - \alpha x_n^2.
    (i) Suppose (x_n) \to \xi. Show that \xi is either 0 or \alpha^{-1}
    (ii) Show that (x_n) \to \alpha^{-1} for all x_1 sufficiently close to \alpha^{-1}.

Section A continued.

  1. [5 points] The conversion from degrees Fahrenheit to degrees Celsius is linear: f(x) = a x + b where x is the temperature in Fahrenheit and f(x) is the corresponding temperature in Celsius. Given f(-40) = -40 and f(32) = 0, what is 68^\circ F in Celsius?
  1. [5 points] Show that approximating the following integral with the composite Trapezoid rule on the n-1 intervals [k,k+1] for k=1,\dots,n-1 gives

\begin{align} I = \int_1^n \log(x) \mathrm{d}x \approx \log n! - \frac12 \log n \tag{A.1} \end{align}

Section A continued.

  1. [5 points] The error function is defined as

\begin{align} \mathrm{erf}(x) = \frac2{\sqrt\pi} \int_0^x e^{-t^2} \mathrm{d}t. \tag{A.2} \end{align}

What is the approximation to \mathrm{erf} if you use the (i) rectangular rule and the (ii) trapezoid rule on the single interval [0,x]?

Recall that Gauss quadrature rules are of the form \int_{-1}^1 f(x) \mathrm{d}x \approx \sum_{j=0}^n w_j f(x_j) where X = \{ x_0, \dots, x_n \} are the zeros of the Legendre polynomial P_{n+1}. This rule is exact for all polynomials of degree at most 2n+1.

  1. [5 points] Notice that P_{n+1}(x) = \prod_{j=0}^n(x-x_j) and consider q(x) := \prod_{j=1}^n(x-x_j). Using the fact that P_{n+1} is orthogonal to q, we have

\begin{align} 0 &= \int_{-1}^1 P_{n+1}(x) q(x) \mathrm{d}x % = \int_{-1}^1 (x-x_0) \prod_{j=1}^n (x-x_j)^2 \mathrm{d}x. \tag{A.3} \end{align}

Use this to show that x_0 \in [-1,1]. (The same argument applies to all the roots so X\subset [-1,1]).

Section A continued.

  1. [5 points] Use the polynomial \ell_X(x)^2 to show that the n+1 point Gauss quadrature rule is not exact for all degree 2n+2 polynomials.

Choose TWO sections from B, C, D (each worth 30%)

29.2 Section B: Lebesgue Constants

We define \|f\|_{L^\infty} := \max_{x\in[a,b]} |f(x)| and let \mathcal P_n be the set of all polynomials of degree at most n. We denote by I_Xf the degree n polynomial interpolation of f \colon [a,b] \to \mathbb R on X = \{x_0 < \dots < x_n\}: i.e. we have I_Xf \in \mathcal P_n and I_Xf(x_j) = f(x_j) for all j=0,\dots,n (in class, we simply called this p_n but here we want to explicitly specify f in the notation).

  1. [2 points] Write down the formula for I_Xf in terms of Lagrange polynomials,

The Lebesgue constant for X is the smallest number \Lambda_n = \Lambda_n(X) such that

\begin{align} \| I_Xf \|_{L^\infty} \leq \Lambda_n \| f \|_{L^\infty}. \tag{B.1} \end{align}

  1. [5 points] Show that

\begin{align} \| I_Xf \|_{L^\infty} \leq \Big( \max_{x\in[a,b]} \sum_{j=0}^n |\ell_j(x)| \Big) \| f \|_{L^\infty} \tag{B.2} \end{align}

and explain why this means that \Lambda_n \leq \max\limits_{x \in [a,b]} \sum_{j=0}^n |\ell_j(x)|.

  1. [3 points] Use question 1 to show that I_X(f + g) = I_X f + I_X g and I_X( \alpha f ) = \alpha I_X f for all functions f,g and constants \alpha.
  2. [5 points] Explain why I_X p_n = p_n for all polynomials p_n of degree at most n.

Section B continued.

Suppose that p_n^\star is the polynomial of degree n such that

\begin{align} \| f - p_n^\star \|_{L^\infty} = \min_{p_n \in \mathcal P_n} \| f - p_n \|_{L^\infty} \tag{B.3} \end{align}

i.e. p_n^\star is the best polynomial approximation to f in \mathcal P_n in the sense that \| f - p_n^\star\|_{L^\infty} \leq \| f - p_n\|_{L^\infty} for all polynomials p_n \in \mathcal P_n.

Notice, by the triangle inequality, we have

\begin{align} \| f - I_Xf \|_{L^\infty} % &= \| f - p_n^\star + p_n^\star - I_Xf \|_{L^\infty} \nonumber\\ % &\leq \| f - p_n^\star\|_{L^\infty} + \| p_n^\star - I_Xf \|_{L^\infty}. \tag{B.4} \end{align}

  1. [10 points] Use \textrm{(B.1)}, \textrm{(B.4)} and questions 3 and 4 to show that

\begin{align} \| f - I_Xf \|_{L^\infty} % \leq (1 + \Lambda_n) \| f - p_n^\star\|_{L^\infty}. \tag{B.5} \end{align}

  1. [5 points] In the following, we plot the Lebesgue constants for interpolation in (i) Equispaced points on [-1,1], and (ii) in Chebyshev nodes. Which set of constants grows faster as a function of n? Explain your answer with reference to \textrm{(B.5)}.

Choose TWO sections from B, C, D (each worth 30%)

29.3 Section C: Hermite cardinal functions and cubic splines

In lectures, we saw the Lagrange cardinal polynomials \ell_j(x) (polynomials with \ell_j(x_k) = \delta_{jk}) and, in A8, we saw the cardinal functions for piecewise linear interpolation (we called these “hat functions”). In this part, we consider the cardinal functions for Hermite interpolation on \{ t_0, t_1, t_2, t_3 \} where t_0 = t_1 := 0 and t_2 = t_3 := 1.

We define the Hermite cardinal functions H_0, H_1, K_0, K_1 : [0,1] \to \mathbb R as the cubic polynomials such that

\begin{align} H_0(0) = 1, \quad H_0(1) = 0, \quad H'_0(0) = 0, \quad H'_0(1) = 0, \tag{C.1}\\ H_1(0) = 0, \quad H_1(1) = 1, \quad H'_1(0) = 0, \quad H'_1(1) = 0, \tag{C.2}\\ K_0(0) = 0, \quad K_0(1) = 0, \quad K'_0(0) = 1, \quad K'_0(1) = 0, \tag{C.3}\\ K_1(0) = 0, \quad K_1(1) = 0, \quad K'_1(0) = 0, \quad K'_1(1) = 1. \tag{C.4} \end{align}

  1. [4 points] Explain in words why these conditions uniquely define H_0, H_1, K_0, K_1.

Hint: To help visualise, here is a plot of these functions


Left: plots of H_{i}. Here, you can see that H'_{i}(0) = H'_{i}(1) = 0 for i=0,1,
Right: plots of K_i. Here, we also plot x and x-1 in dashed line (so we can see that K_0'(0)=1 and K_1'(1)=1).

Section C continued.

  1. [4 points] Write down the finite difference table for H_0(t) (you may fill in the table below)

\begin{align*} \begin{matrix} & H_0[t_i] & & H_0[t_i, t_{i+1}] & & H_0[t_i, t_{i+1},t_{i+1}] & & H_0[t_0,t_1,t_2,t_3] \\ \\ t_0 = 0: & & & & \\ \\ t_1 = 0: & & & & \\ \\ t_2 = 1: & & & & \\ \\ t_3 = 1: & & & & \end{matrix} \tag{C.5} \end{align*}

  1. [4 points] Hence or otherwise show that H_0(t) = 2t^3 - 3 t^2 + 1.

A similar argument (you do not need to prove) gives

\begin{align} H_1(t) &= -2t^3 + 3t^2, \qquad K_0(t) = t^3 - 2t^2 + t, \qquad K_1(t) = t^3 - t^2. \tag{C.6} \end{align}

Let f : [-1,1] \to \mathbb R be some function and X = \{x_0 < \dots < x_n\} be equispaced points with x_0 = -1 and x_n = 1. We define h := \frac2n = x_{i+1}-x_i.

In this question we will prove that there exists an S:[-1,1]\to\mathbb R such that

(i) S is two times continuously differentiable,
(ii) S(x_i) = f(x_i) for i = 0,\dots,n,
(iii) S is a cubic polynomial on each subinterval [x_{i},x_{i+1}], and
(iv) S'(x_0) = f'(x_0) and S'(x_n) = f'(x_n).

Such a function is called a cubic spline (with clamped boundary conditions).

  1. [6 points]
    Fix an interval [x_i, x_{i+1}] for some i=0,\dots,n-1 and consider the change of variables t = \frac{x-x_i}{h}. Show that the function

\begin{align} S_i(x) := f(x_i) H_0(t) + f(x_{i+1}) H_1(t) + h \big( m_i K_0(t) + m_{i+1} K_1(t) \big) \tag{C.7} \end{align}

is a cubic polynomial with S_i(x_i) = f(x_i), S_i(x_{i+1}) = f(x_{i+1}), S'_i(x_i) = m_i, and S'_i(x_{i+1}) = m_{i+1}. (Recall that we defined H_0, H_1, K_0, K_1 in \textrm{(C.1)}-\textrm{(C.4)}, above.)

To obtain a cubic spline S : [-1,1] \to \mathbb R, we define S(x) = S_i(x) if x \in [x_i, x_{i+1}].

  1. [2 points] Explain why m_0 = f'(x_0) and m_{n} = f'(x_n).
  2. [2 points] Explain why S and S' are continuous on [-1,1],

Section C continued.

One can show that

\begin{align} h^2 S''_i(x_i) &= 6 ( f(x_{i+1}) - f(x_i) ) - 2 h \big( m_{i+1} + 2m_i \big) \nonumber\\ h^2 S''_i(x_{i+1}) &= -6( f(x_{i+1}) - f(x_i) ) + 2 h \big( 2m_{i+1} + m_{i} \big). \tag{C.8} \end{align}

In order to enforce continuity of the second derivatives of S, we require S_{i-1}''(x_{i}) = S_{i}''(x_{i}) for each i = 1,\dots,n-1 which, by \textrm{(C.8)}, is equivalent to

\begin{align} m_{i-1} + 4 m_i + m_{i+1} = 3 \frac{ f(x_{i+1}) - f(x_{i-1}) }{h} \tag{C.9} \end{align}

(you do not need to prove this).

  1. [4 points] Show that this linear system of equations can be written as A \bm m = \bm b where

\begin{align} A &= \begin{pmatrix} 4 & 1 & & \\ 1 & 4 & \ddots & \\ & \ddots & \ddots & 1 \\ & & 1 & 4 \end{pmatrix}. % \tag{C.10} %\\ % %b_i &= 3 \frac{f(x_{i+1}) - f(x_{i-1})}{h} - %\begin{cases} % f'(x_0) & \text{if } i=1\nonumber\\ % 0 & \text{if } 2 \leq i \leq n-2\\ % f'(x_n) & \text{if } i = n-1. %\end{cases} %\tag{C.12} \end{align}

It turns out that A is invertible and so there exists a unique \bm m which satisfies A\bm m = \bm b and so there exists a unique cubic spline satisfying the conditions (i)-(iv), above. Moreover, we can compute the cubic spline by using m = A \ b in Julia and evaluating \textrm{(C.7)}.

  1. [4 points] For smooth functions f, this cubic spline interpolation converges with rate O(n^{-4}). In the following we compute the cubic spline interpolation of f(x) = \cos(10 x) on [-1,1] and plot the approximation errors as a function of n. What scaling do we use for the graph (i.e. linear, log-y, or log-log)?

Choose TWO sections from B, C, D (each worth 30%)

29.4 Section D: Linear systems of equations

29.5 LU decomposition

  1. [6 points] Find the LU decomposition of

\begin{align} \begin{pmatrix} 2 & 1 & 0 & 0 \\ 1 & 2 & 1 & 0 \\ 0 & 1 & 2 & 1 \\ 0 & 0 & 1 & 2 \end{pmatrix} \tag{D.1} % % = \begin{pmatrix} % 1 & 0 & 0 & 0 \\ % 1/2 & 1 & 0 & 0 \\ % 0 & 2/3 & 1 & 0 \\ % 0 & 0 & 3/4 & 1 % \end{pmatrix} % \begin{pmatrix} % 2 & 1 & 0 & 0 \\ % 0 & 3/2 & 1 & 0 \\ % 0 & 0 & 4/3 & 1 \\ % 0 & 0 & 0 & 5/4 % \end{pmatrix} \end{align}

  1. [6 points] Hence solve

\begin{align} \begin{pmatrix} 2 & 1 & 0 & 0 \\ 1 & 2 & 1 & 0 \\ 0 & 1 & 2 & 1 \\ 0 & 0 & 1 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix} % = \begin{pmatrix} 12 \\ -24 \\ 8 \\ 26 \end{pmatrix} \tag{D.2} \end{align}

  1. [3 points] You may have noticed that when A is tridiagonal (as in the case of Question 1, above), then L and U are bidiagonal with L only having nonzero entries on the diagonal and lower diagonal and U only having nonzero entries on the diagonal and upper diagonal. Explain why this is true in general.

  2. [3 points] If A is a tridiagonal matrix, what is the computational cost in solving A \bm x = \bm b using LU decomposition and forwards and backwards substitution?

O(n)

O(n^2)

O(n^3)

Section D continued.

29.6 Singular value decomposition (SVD)

Recall that the SVD of A \in \mathbb R^{m\times n} is A = P\Sigma Q^T where P\in \mathbb R^{m\times m} is orthogonal (PP^T = P^TP = I), \Sigma\in \mathbb R^{m\times n} is diagonal with singular values \sigma_i := \Sigma_{ii}, \sigma_1 \geq \sigma_2 \geq \dots \geq \sigma_r \geq 0, and r := \min\{m,n\}, and Q \in \mathbb R^{n\times n} is orthogonal.

  1. [4 points] What are the singular values of the following matrix:

\begin{align} \begin{pmatrix} 1 & 0 \\ 1 & 0 \\ 1 & 0 \end{pmatrix}.\tag{D.3} \end{align}

Recall that \|A\| := \max_{\bm x\not= 0} \frac{|A\bm x|}{|\bm x|} where |\bm x|:= \sqrt{\sum_{i} x_i^2}.

  1. [4 points] Use the SVD of A, together with the fact that |u|^2 = \sum_{i=1}^m u_i^2 = u^T u, to show that

\begin{align} \|A\|^2 % = \max_{\bm y \not= 0} \frac {|P\Sigma y|^2} {|y|^2} % = \max_{\bm y \not= 0} \frac{\sum_{i} \sigma_i^2 y_i^2}{\sum_i y_i^2}.\tag{D.4} \end{align}

  1. [2 points] Hence show that \| A \| = \sigma_1.

  2. [2 points] Two matrices A, B \in \mathbb R^{512\times512} represent the following two grayscale images. Which matrix has singular values which decay faster? Why?