30  Final - (my) solutions

---

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.

---

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

Answer.

The relative error is

\begin{align} \left|\frac{c - 3\times10^8}{c} \right| % &= \frac{3 - 2.9979}{2.9979} \end{align}

2. [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?

Answer.

At least quadratic: there exists \eta_n such that

\begin{align} x_{n+1} -\xi = g(x_n) - g(\xi) = \frac{g''(\eta_n)}{2}(x_n-\xi)^2. \end{align}

3. [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?

Answer.

y because we are doing fewer floating point operations. Each iteration in the formula for x introduces errors.

Remark.

Here, we compute the value in Float64 using the built-in method (which doesn’t do 2025 repeated multiplications):

1.001^2025

7.568449119215128

4. [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]?

Answer.

(i) f is a continuous function on [0,\frac\pi2] that changes sign: f(0)= 1 and f(\frac\pi2) = -\pi. Therefore, there exists a root \xi \in (0, \frac\pi2).

(ii) We have \xi = \frac{1}2 \cos \xi and g:[0,\frac\pi2]\to[0,\frac12]\subset [0,\frac\pi2] defined by g(x) = \frac12\cos x is a contraction: |g'(x)| = |-\frac12 \sin x| \leq \frac12 < 1.

(iii) Yes, by the contraction mapping theorem.

5. [2 points] True or false: there exists a continuous function f with a root in [a,b] and f(a) f(b) > 0? Why?

Answer.

True. Any continuous function with f(a), f(b) \not= 0 and exactly two roots in [a,b] satisfies these conditions.

6. [2 points] True or false: all continuous functions f: [a,b] \to \mathbb R have fixed points \xi = f(\xi)? Why?

Answer.

False. For example f(x) = x + 1 has no fixed points.

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

Answer.

(i) Because x \mapsto 2x - \alpha x^2 is continuous, we have \xi = 2 \xi - \alpha \xi^2 and so \alpha\xi^2 = \xi. Thus \xi = 0 or \xi = \alpha^{-1}.

(ii) Define g(x) = 2x - \alpha x^2. Since g'(\alpha^{-1}) = 2 - 2\alpha \alpha^{-1} = 0, there exists an interval around \xi on which g is a contraction. Therefore, by the contraction mapping theorem, (x_{n}) \to \alpha^{-1} as n \to \infty,

8. [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?

Answer.

Let’s consider f as the interpolation of \{ 0, -40 \} on X = \{32, -40\}. In Newton form, we have

\begin{align} f(x) &= f(32) + \frac{f(32) - f(-40)}{32-(-40)} (x-32) % = \frac{40}{72}(x-32). \end{align}

When x = 68, we have f(68) = \frac{40}{72} \cdot 36 = 20.

9. [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}

Answer.

We have

\begin{align} I &= \sum_{k=1}^{n-1} \int_{k}^{k+1} \log(x) \mathrm{d}x\nonumber\\ % &\approx \sum_{k=1}^{n-1} \frac{(k+1)-k}{2}\left[ \log k + \log(k+1) \right] \nonumber\\ % &= \frac{1}{2}\left[ (\log 1 + \log 2) + (\log 2 + \log 3) + \dots + (\log (n-1) + log n) \right] \nonumber\\ % &= \log 1 +\log 2 \dots + \log n - \frac12 \log n \nonumber\\ % &= \log n! - \frac12 \log n. \end{align}

Remark.

Stirling’s approximation

10. [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]?

Answer.

(i) We have

\begin{align} \frac2{\sqrt\pi} \int_0^x e^{-t^2} \mathrm{d}t &\approx \frac2{\sqrt\pi} \left[ (x-0) e^{-0^2} \right] \nonumber\\ % &= \frac{2x}{\sqrt\pi} \end{align}

(ii) In this case, we get

\begin{align} \frac2{\sqrt\pi} \int_0^x e^{-t^2} \mathrm{d}t &\approx \frac2{\sqrt\pi} \frac{x-0}{2} \left[ e^{-0^2} + e^{-x^2} \right] \nonumber\\ % &= \frac{x}{\sqrt\pi} \left( 1 + e^{-x^2} \right). \end{align}

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.

11. [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]).

Answer.

Notice that, if x_0 \not\in [-1,1], then x\mapsto x-x_0 is either positive or negative on the whole interval [-1,1] and so we have

\begin{align} \left| \int_{-1}^1 (x-x_0) \prod_{j=1}^n (x-x_j)^2 \mathrm{d}x \right| > 0, \end{align}

contradicting the fact this is zero. As a result, we have x_0 \in [-1,1].

12. [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.

Answer.

First note that \ell_X(x)^2 is a polynomial of degree 2(n+1) = 2n+2. If the Gauss quadrature rule is exact for \ell_X(x)^2, then

\begin{align} \int_{-1}^1 \ell_X(x)^2 \mathrm{d}x % &= \sum_{j=0}^n w_j \ell_X(x_j)^2 = 0. \end{align}

However, the integral of a non-zero positive polynomial on a finite interval is strictly postive, a contradiction.

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

30.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,

Answer.

I_Xf(x) = \sum_{j=0}^n f(x_j) \ell_j(x) where \ell_j are the Lebesgue polynomials: \ell_j(x) = \prod_{k\not=j} \frac{x-x_k}{x_j-x_k}.

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}

2. [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)|.

Answer.

We have

\begin{align} \| I_Xf \|_{L^\infty} &= \max_{x\in[a,b]} \left|\sum_{j=0}^n f(x_j) \ell_j(x)\right| \nonumber\\ % &\leq \max_{x\in[a,b]} \sum_{j=0}^n \left|\ell_j(x)\right| \left(\max_j |f(x_j)|\right) \nonumber\\ % &\leq \Big( \max_{x\in[a,b]} \sum_{j=0}^n |\ell_j(x)| \Big) \| f \|_{L^\infty}. \end{align}

Since C := \Big( \max_{x\in[a,b]} \sum_{j=0}^n |\ell_j(x)| \Big) is a number for which \| I_Xf \|_{L^\infty} \leq C \| f \|_{L^\infty} and \Lambda_n is the smallest such number, we must have \Lambda_n \leq C.

Remark.

We actually have equality (this is the formula we use to plot the Lebesgue constants below).

3. [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.

Answer.

We have

\begin{align} I_X( f + \alpha g ) &= \sum_{j=0}^n (f+\alpha g)(x_j) \ell_j \nonumber\\ % &= \sum_{j=0}^n \left[ f(x_j) +\alpha g(x_j) \right] \ell_j \nonumber\\ % &= \sum_{j=0}^n f(x_j) \ell_j + \alpha \sum_{j=0}^n g(x_j) \ell_j \nonumber\\ % &= I_X f + \alpha I_X g. \end{align}

(this is equivalent to showing both I_X( f + g ) = I_Xf + I_Xg and I_X (\alpha f) = \alpha I_X f).

4. [5 points] Explain why I_X p_n = p_n for all polynomials p_n of degree at most n.

Answer.

I_X p_n is a polynomial of degree at most n for which I_Xp_n(x_j) = p_n(x_j) for all j=0,\dots,n-1. Since polynomials of degree n are determined by their values on n+1 distinct points (the polynomial I_X p_n - p_n has n+1 distinct zeros so must be the zero polynomial), we must have I_X p_n = p_n.

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}

5. [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}

Answer.

We have

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

We have used the fact that p_n^\star = I_X p_n^\star, I_X p_n^\star - I_Xf = I_X(p_n^\star-f), and the definition of \Lambda_n.

6. [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)}.

Answer.

The red dots correspond to the Chebyshev points and the black dots correspond to equispaced points. We saw in lectures that polynomial interpolation in equispaced points can blow up (even for smooth functions), whereas interpolation in Chebshev points converges for smooth functions. We would therefore expect the constant of proportionality between interpolation and best approximation to be much larger for equispaced points.

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

Answer.

These polynomials are the solutions to Hermite interpolation problems of the form we saw in lectures. That is, the interpolation set is X = \{\{0,0,1,1\}\} and we are specifying the values of the function and its derivative at these points. In class we saw that there there is a unique solution to this problem.

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

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

Answer.

\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: & 1 & & H_0'(0) = 0 & & \frac{-1 - 0}{1-0} =-1 & & \frac{1-(-1)}{1-0} = 2 \\ \\ t_1 = 0: & 1 & & \frac{0-1}{1-0} = -1 & & \frac{0-(-1)}{1-0} = 1 & & \\ \\ t_2 = 1: & 0 & & H_0'(1) = 0 & & \\ \\ t_3 = 1: & 0 & & & & \end{matrix} \tag{C.5} \end{align*}

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

Answer.

We have

\begin{align} H_0(t) &= H_0[0] + H_0[0,0] (t-0) + H_{0}[0,0,1] (t-0)^2 + H_0[0,0,1,1] (t-0)^2 (t-1) \nonumber\\ % &= 1 + 0 t - t^2 + 2 t^2 (t-1) \nonumber\\ % &= 2t^3 - 3t^2 + 1 \end{align}

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

4. [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.)

Answer.

Notice that, when x = x_i we have t = 0 and, when x = x_{i+1}, we have t = 1. Moreover, we have \mathrm{d}x = h^{-1} \mathrm{d}t. As a result, we have

\begin{align} S_i(x_i) &= f(x_i) H_0(0) + f(x_{i+1}) H_1(0) + h \big( m_i K_0(0) + m_{i+1} K_1(0) \big) \\ S_i(x_{i+1}) &= f(x_i) H_0(1) + f(x_{i+1}) H_1(1) + h \big( m_i K_0(1) + m_{i+1} K_1(1) \big) \\ S_i'(x_i) &= f(x_i) h^{-1} H_0'(0) + f(x_{i+1}) h^{-1} H_1'(0) + m_i K_0'(0) + m_{i+1} K_1'(0) \\ S_i'(x_{i+1}) &= f(x_i) h^{-1} H_0'(1) + f(x_{i+1}) h^{-1} H_1'(1) + m_i K_0'(1) + m_{i+1} K_1'(1) \end{align}

Using the properties of H_0, H_1, K_0, K_1, we find that 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}.

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

5. [2 points] Explain why m_0 = f'(x_0) and m_{n} = f'(x_n).

Answer.

We have S'(x_0) = S'_0(x_0) = m_i and S'(x_{n}) = S'_{n-1}(x_n) = m_{n} and we require S'(x_0) = f'(x_0) and S'(x_n) = f'(x_n).

6. [2 points] Explain why S and S' are continuous on [-1,1],

Answer.

Since S_i is a polynomial, we have S, S' continuous on each interval [x_i, x_{i+1}]. We just need to show that they are continuous across the boundaries. We have S_i(x_{i+1}) = f(x_{i+1}) = S_{i+1}(x_{i+1}) and so S is continuous. Moreover, we have S_i'(x_{i+1}) = m_{i+1} = S_{i+1}'(x_{i+1}) and so S' is continuous.

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

7. [4 points] Show that \mathrm{(C.9)} for i=1,\dots,n-1, 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}

Answer.

Since m_0 = f'(x_0), we have

\begin{align} \begin{pmatrix} 4 & 1 & 0 & \dots & 0 \end{pmatrix} \begin{pmatrix} m_1 \\ m_2 \\ \vdots \\m_{n-1} \end{pmatrix} = b_1 := 3 \frac{ f(x_2) - f(x_{0}) }{h} - f'(x_0). \end{align}

For i = 2,\dots, n-1, \mathrm{C.9} is equivalent to

\begin{align} \begin{pmatrix} 0 & \dots & 1 & 4 & 1 & \dots & 0 \end{pmatrix} \begin{pmatrix} \vdots \\ m_{i-1} \\ m_i \\ m_{i+1}\\ \vdots \end{pmatrix} = b_i := 3 \frac{ f(x_{i+1}) - f(x_{i-1}) }{h} . \end{align}

Finally, since m_n = f'(x_n), we have

\begin{align} \begin{pmatrix} 0 & \dots & 0 & 1 & 4 \end{pmatrix} \begin{pmatrix} m_1 \\ m_2 \\ \vdots \\m_{n-1} \end{pmatrix} = b_{n-1} := 3 \frac{ f(x_n) - f(x_{n-2}) }{h} - f'(x_n). \end{align}

Putting this together, we get A \bm m = \bm b.

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

8. [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)?

Answer.

It’s a log-log plot. Since y = n^{-4}, we have \log y = - 4 \log n and so we can see a straight line in a log-log plot.

30.4 Section D: Linear systems of equations

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

2. [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}

3. [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.

4. [3 points] If A \in \mathbb R^{n\times n} 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.

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

5. [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}.

6. [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}

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

8. [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?