📢 Updates

Lecture 9: Hermite interpolation available Assignment 3 available Lecture 8 online Lecture 7 online Schedule now available in "Admin" Problem class exercises now online Lecture 6 online Quiz 1 solutions online Lectures 2-3 exercises have solutions! Lecture 1 exercises have solutions!

Exercise bank

Lecture 1

Exercise 1.1

Add a new markdown cell below this one and write “Hello world!”

Exercise 1.2

Compute \(\eta\) for approximating \(e\) with \(3\), \(2.7\), and \(2.718\).

We have

Code

η( x, x_tilde ) = - log10( abs(x-x_tilde)/abs(x) )
η( exp(1), 3 ), η( exp(1), 2.7 ), η( exp(1), 2.718 )

(0.9844796209116626, 2.172274852323722, 3.984309635847452)

Exercise 1.3

Show that \(\frac{n + \cos(n)}{n^2} = O(\frac1n)\) and \(\frac{n^2+3}{n^{42}} = O(\frac1{n^{40}})\).

Since \(\cos(n) \leq n\) for all \(n \in \mathbb{N}\), we have

\[\begin{align} \frac{n + \cos n}{n^2} % \leq \frac{2n}{n^2} = \frac{2}{n}. \end{align}\]

That is, \(\frac{n + \cos n}{n^2}\) is bounded above by a constant multiple of \(n^{-1}\) so \(\frac{n + \cos n}{n^2} = O(n^{-1})\).

Similarly, \(3 \leq 3n^2\) so \(\frac{n^2 + 3}{n^{42}} \leq \frac{4}{n^{40}}\).

Exercise 1.4

Let \(e\) denote Euler’s constant.

1. \(e_N := \big( 1 + \frac{1}{N} \big)^N\) is an approximation \(e\). What is the error, absolute error, and relative error in approximating \(e\) with \(e_N\) for \(N = 1,2,...,10\)? What are the number of accurate digits in this approximation?
2. What is the rate of convergence of \(e_N\) to \(e\)?

Let’s define \(e\) and the sequence of errors \(|e - e_N|\):

Code

e = exp(1)

N = 1:100;
err = @. abs( e - (1 + 1/N)^N );

The error is equal to the absolute error in this case and is given by

Code

err[1:10]

10-element Vector{Float64}:
 0.7182818284590451
 0.4682818284590451
 0.34791145808867485
 0.2768755784590451
 0.22996182845904567
 0.19665545671693163
 0.17178213141833298
 0.1524973145086972
 0.13710703674584712
 0.12453936835904278

Computing the relative errors give

Code

err[1:10] ./ e

10-element Vector{Float64}:
 0.26424111765711533
 0.17227125736425472
 0.12798947277880338
 0.10185683307753335
 0.08459822894427681
 0.07234549952033957
 0.0631951145094156
 0.05610062684160521
 0.05043886005873465
 0.045815473235769066

Since \(0.0458... = 4.58... \times 10^{-2} \leq 5 \times 10^{-2}\), \(e_{10}\) is correct to \(2\) significant (decimal) digits.

One can also compute the \(\eta = - \log_{10} |\frac{e - e_N}{e}|\):

Code

- log10.( err[1:10]/e )

10-element Vector{Float64}:
 0.5779996023832055
 0.7637871764659977
 0.8928257498997241
 0.992009830990522
 1.072638728778693
 1.1405884803607715
 1.199316494876063
 1.2510322861176486
 1.2972347372416126
 1.338987823174671

Looking at the following plot, we can see that \(| e - e_N | = O(N^{-1})\):

Code

plot( err ;
    xaxis=:log, yaxis=:log, xlabel=L"N", title="absolute error", legend=false, lw=2, m=2 )
plot!( N.^(-1), color=:black, linestyle=:dash, label = L"n^{-1}" )

Exercise 1.5

Euler-Mascheroni constant, \(\gamma\), is defined as the limit of \(\gamma_n = \sum_{k=1}^n \frac{1}{k} - \log n\) as \(n \to \infty\).

1. What is the absolute error in approximating \(\gamma\) with \(\gamma_{1000}\)?
2. What is the rate of convergence of \(\gamma_n\) to \(\gamma\)?
3. [This was just so that we can remember some calculus and also rigorously prove the rate of convergence that we have seen above using a rigorous error estimate]
   Prove that \(|\gamma - \gamma_n| \leq \frac1n\).

First, we define the Euler-Mascheroni constant:

Code

γ = Base.MathConstants.eulergamma

γ = 0.5772156649015...

… and compute the absolute error in approximating \(\gamma\) with \(\gamma_{1000}\):

Code

Nmax = 1000;

N = 1:Nmax;
γN = zeros( Nmax );
γN[1] = 1;
for n ∈ 2:Nmax
    γN[n] = (1/n) + γN[n-1] + log( (n-1)/n );
end

println( "γ₁₀₀₀ = ", γN[1000] )
println( "Abs. error = " , abs( γ - γN[1000] ) )

γ₁₀₀₀ = 0.5777155815682078
Abs. error = 0.0004999166666749266

Here’s a plot of the absolute errors on a log-log; the gradient matches \(O(N^{-1})\):

Code

plot( γN .- γ , xaxis=:log, yaxis=:log , legend = false, lw = 2, m=2)   
plot!( N.^(-1) , linestyle=:dash)

Proof. First, notice that

\[\begin{align} \gamma_n - \gamma_N % &= \log \frac Nn - \sum_{k=n+1}^N \frac{1}{k}\\ % &= \sum_{k=n+1}^N \left( \log \frac{k}{k-1} - \frac{1}{k} \right). \end{align}\]

Here, \(\frac1k\) is the area under the rectangle with base \([k-1,k]\) and height \(\frac1k\) and \(\log \frac{k}{k-1}\) is the integral of \(\frac1x\) between \(x = k-1\) and \(x = k\):

Code

k=3; x = range(k-1,k, 100);
plot( x, [x.^(-1) ones(100)./k (x .+ 1).^(-1)], 
    ylims=(0,1/(k-1)), 
    label=[L"y=\frac{1}{x}" L"\frac{1}{k}" L"y=\frac{1}{x+1}"] )

In the graph above, \(\frac1k\) is the area under the red line, \(\log \frac{k}{k-1}\) is the area under the blue line. We upper bound the absolute value of the error between these areas, by computing the area between \(y = \frac1x\) and \(y=\frac1{x+1}\):

\[\begin{align} 0 \leq \gamma_n - \gamma_N % &\leq \sum_{k=n+1}^N \left( \log \frac{k}{k-1} - \frac{1}{k} \right) \\ % &\leq \sum_{k=n+1}^N \int_{k-1}^k \left( \frac{1}{x} - \frac{1}{x + 1} \right) \mathrm{d}x % = \int_{n}^\infty \frac{1}{x(x+1)} \mathrm{d}x \\ % &\leq \int_{n}^\infty \frac{1}{x^2} \mathrm{d}x = \frac1n \end{align}\]

Here, we have shown that \((\gamma_n)\) is a Cauchy sequence and thus it converges, and we have the same estimate for \(\gamma_n - \gamma\) by taking the limit in \(N\).

Exercise 1.6

Assignment 1, Fall 2025, 30 pts The quantities

\[\begin{align} \pi_N := 4 \sum_{n=0}^N \frac{(-1)^n}{2n+1} \end{align}\]

approximate \(\pi\) by using the Taylor expansion of \(\arctan\) together with the fact that \(\pi = 4 \arctan(1)\).

Code

using Plots, LaTeXStrings

Nmax = 100;
Ninter = 1:Nmax;

# vectors of length Nmax
p = zeros( Nmax ); q = zeros( Nmax );
p[1] = 4 * ( 1 - 1/3 ) ;
q[1] = abs( p[1] - π );
for N in 2:Nmax
    p[N] = p[N-1] + 4 * (-1)^N/(2*N + 1);
    q[N] = abs( p[N] - π );
end

plot(p, legend=false, xlabel="N", title=L"\pi_N")

Based on the following graph, what is the rate of convergence of this particular approximation to \(\pi\)? Why do we use a log-log scale for this graph?

Code

plot( q , xaxis=:log, yaxis=:log, xlabel=L"N", title="absolute error", legend=false)
plot!( Ninter.^(-1) , ls=:dash)

We have \(|\pi - \pi_N| = O(N^{-1})\). The log-log plot makes this behaviour easy to see.

In fact, one can show that

\[\begin{align} |\pi - \pi_N| &= 4 \int_{0}^1 \frac{t^{2(N+1)}}{1 + t^2} \mathrm{d}t\nonumber\\ % &\leq \frac{4}{2N+3}. \end{align}\]

Lecture 2

Exercise 2.1

In England, when we need to convert between lb and kg some people use the phrase “2.2 pounds of jam is the same as a kilogram”.

• What is the error, absolute error, and relative error of approximating 1 kg with 2.2 lb?
• To how many significant digits is this approximation correct to?

You may use 1 kg = 2.20462262185 lb.

The error is simply “exact minus approx.” \(= 2.20462262185 - 2.2 \approx 0.00462262185\). This is also the absolute error (because it is positive).

For the relative error, we divide by the exact value:

Code

x = 2.20462262185 
y = 2.2 

println("relative error = $( (x-y)/x )")

relative error = 0.0020967860005540483

The error is \(0.002096... \leq 5 \times 10^{-3}\) and so \(2.2\) is correct to 3 significant digits. (This is sufficient for all recipes you will likely use).

Exercise 2.2

Most people say a marathon is 26 miles when it is actually 26 miles 385 yards. In Europe, people say a marathon is 42 km when it is actually 42.194988 km.

• Which approximation has a smaller error?
• What is the advantage of computing the relative errors in this case?

Hint: You may use the fact that 1 mile \(=\) 1760 yards.

A marathon is \(26 \times 1760 + 385\) yards. The error and relative error in this approximation is

Code

x = 26*1760 + 385
y = 26*1760

err, r_err = x-y, (x-y)/x

(385, 0.00834326579261025)

The same calculation for the approximation in km is

Code

x′ = 42.194988
y′ = 42 

err′, r_err′ = x′-y′, (x′-y′)/x′

(0.19498800000000216, 0.004621117560218341)

Code

r_err / r_err′

1.8054649516893149

The approximation in miles has a \(1.8\times\) larger relative error than the approximation in km.

Why is the error in the approximation in miles much larger than \(1.8\times\) bigger than the error in the approximation in km?

Exercise 2.3

The astronomical unit is defined as \(1 \text{AU} = 149,597,870,700 \text{m}\). Write this in standard form.

\(1.495978707 \times 10^{11} \text{m}\)

Exercise 2.4

Write down all the floating point numbers with \(k = 2\) in the interval \([1,2)\). What about the interval \([2,4)\)? \([2^n,2^{n+1})\)?

There are \(2^k\) equally spaced floating point numbers in any interval \([2^n,2^{n+1})\) (as long as \(\ell\) is sufficiently large).

Floats between \(1\) and \(2\) are \((1.b_1b_2)_2 = 1 + b_1 2^{-1} + b_2 2^{-2}\) and so are \(1\) (\(b_1 = b_2 = 0\)), \(1.25\) (\(b_1 = 0, b_2 = 1\)), \(1.5\) (\(b_1 = 1, b_2 = 0\)), \(1.75\) (\(b_1 = b_2 = 1\)).

Similarly, floats in the interval \([2^n, 2^{n+1})\) are of the form \((1 + b_1 2^{-1} + b_2 2^{-2} ) 2^n\). That is, \(2^n\) (\(b_1 = b_2 = 0\)), \(1.25 \times 2^n\) (\(b_1 = 0, b_2 = 1\)), \(1.5 \times 2^n\) (\(b_1 = 1, b_2 = 0\)), \(1.75 \times 2^n\) (\(b_1 = b_2 = 1\)).

Exercise 2.5

Take \(k=\ell=2\).

• What number does \(0 0 1 1 0\) encode?
• Write the bitstring of \(-3.5\).
• What are the maximum and minimum positive (normal) numbers and the machine precision in this system?

Recall that floats are stored as \(s | c | f\) where \(s\) is the sign bit, \(c\) the characteristic, and \(f\) the fraction or mantissa. Recall that the bias is \(b = 2^{\ell-1}-1 = 2 - 1 = 1\) and \(c \in \{1, \dots, 2^\ell-2 \} = \{1,2\}\).

In this case, we have \(s = 0\), \(c = (01)_2 = 1\) and \(f = (0.10)_2 = \tfrac12\) and so this number represents \((-1)^s (1 + f) 2^{c - b} = (-1)^0 (1 + \tfrac12) 2^{1-1} = \tfrac32\).

\(-3.5 = (-1)^1 (1 + 0.75) 2 = (-1)^1 (1 + (0.11)_2 ) 2^{2-1}\) and so the bitstring is \(1 1 0 1 1\).

The maximum is \((-1)^0 (1 + (0.11)_2 ) 2^{ 2 - 1} = \tfrac72\),

The minimum positive float is \((-1)^0 (1 + (0.00)_2) 2^{ 1 - 1 } = 1\),

Machine precision is \(2^{-k} = \tfrac14\).

Exercise 2.6

What are the maximum and minimum positive numbers and the machine precision for each of Float16, Float32, and Float64?

Recall, the bitstring is stored as \(s | c | f\) with \(1 | \ell | k\) bits for each where

• Float16 = half precision = \(1|5|10\),
  
• Float32 = single precision = \(1|8|23\), and
  
• Float64 = double precision = \(1|11|52\).

In general, the maximum number is \(2^{2^{\ell-1} -1}( 2 - 2^{-k} ) \approx 2^{2^{\ell-1}}\), the minimum positive (normal) number is \(2^{2-2^{\ell-1}}\), and the machine precision is \(\epsilon_{\text{mach}} = 2^{-k}\). Therefore:

• Float16 (\(\ell = 5\), \(k = 10\)): the maximum is \(2^{15}(2 - 2^{-10}) = 65504\), the minimum positive (normal) number is \(2^{-14} \approx 6.10 \times 10^{-5}\), and the machine precision is \(2^{-10} \approx 9.77 \times 10^{-4}\).
• Float32 (\(\ell = 8\), \(k = 23\)): the maximum is \(2^{127}(2 - 2^{-23}) \approx 3.40 \times 10^{38}\), the minimum positive (normal) number is \(2^{-126} \approx 1.18 \times 10^{-38}\), and the machine precision is \(2^{-23} \approx 1.19 \times 10^{-7}\).
• Float64 (\(\ell = 11\), \(k = 52\)): the maximum is \(2^{1023}(2 - 2^{-52}) \approx 1.80 \times 10^{308}\), the minimum positive (normal) number is \(2^{-1022} \approx 2.23 \times 10^{-308}\), and the machine precision is \(2^{-52} \approx 2.22 \times 10^{-16}\).

You can check these in Julia with floatmax, floatmin, and eps:

Code

for T ∈ (Float16, Float32, Float64)
    label = "$T:  "
    pad = " "^length(label)
    println(label, "max = ", floatmax(T), "\n", pad, "min = ", floatmin(T), "\n", pad, "eps = ", eps(T))
end

Float16:  max = 6.55e4
          min = 6.104e-5
          eps = 0.000977
Float32:  max = 3.4028235e38
          min = 1.1754944e-38
          eps = 1.1920929e-7
Float64:  max = 1.7976931348623157e308
          min = 2.2250738585072014e-308
          eps = 2.220446049250313e-16

Exercise 2.7

Write down the bitstring for \(\pi\) in Float16.

\(\pi\) in Float16 rounds to

Code

x = Float16( π )
x == 3.140625

true

and, in Float16, we have \(1|\ell|k = 1|5|10\). Therefore, the bias is \(2^{\ell-1}-1 = 2^{4}-1 = 15\)

Since \(3.140625 = (-1)^0 (1 + 1/2 + 0/4 + 0/8 + 1/16 + 0/32 + 0/64 + 1/128)\times 2^{16-15}\), the bitstring is \(0 | 1 0 0 0 0 | 1 0 0 1 0 0 1 0 0 0\)

Code

bitstring(x)

"0100001001001000"

Exercise 2.8

Compute the relative condition numbers of (i) \(x \mapsto cx\) for \(c \not= 0\), (ii) \(\log(x)\), (iii) \(\exp\), (iv) \(x^p\), and (v) \(f\circ g\) in terms of \(\kappa_f\) and \(\kappa_g\)

Recall that \(\kappa_f(x) = \left| \frac{x f'(x)}{f(x)} \right|\).

Since \((x \mapsto cx)' = c\), \(\log'(x) = \tfrac1x\), \(\exp'(x) = \exp(x)\), \((x^p)' = p x^{p-1}\), and \((f \circ g)'(x) = f'(g(x)) g'(x)\), we have

\[\begin{align} \kappa_{x \mapsto c x}(x) &= \left| \frac{x c}{cx} \right| = 1 \nonumber\\ % \kappa_{\log}(x) &= \left| \frac{x \cdot \tfrac1x}{\log(x)} \right| = \frac{1}{|\log x|} \qquad \text{for } x \not= 1, \nonumber\\ % \kappa_{\exp}(x) &= \left| \frac{x e^x}{e^x} \right| = |x| \nonumber\\ % \kappa_{x\mapsto x^p}(x) &= \left| \frac{x \cdot p x^{p-1}}{x^p} \right| = p \nonumber \end{align}\]

and

\[\begin{align} \kappa_{f\circ g}(x) &= \left| \frac{x f'(g(x)) g'(x)}{f(g(x))} \right| \nonumber\\ % &=\left| \frac{g(x) f'(g(x))}{ f(g(x))} \right| \left| \frac{x g'(x)}{ g(x) } \right| \nonumber\\ % &= \kappa_{f}(g(x)) \kappa_g(x) \nonumber \end{align}\]

Exercise 2.9

Midterm, Fall 2025, Compulsory Section, 3 pts. Let us define

\[\begin{align*} f(x) = \frac{\sqrt{1+x}-1}{x}, \qquad g(x) = \frac{1}{\sqrt{1+x} + 1}. \end{align*}\]

One may verify that \(f\) and \(g\) are mathematically identical (but you don’t have to). The following table is the output of evaluating \(f\) and \(g\) on \(\{10^{-14}, 10^{-15}, 10^{-16}\}\) in double precision floating point. Why are the function values different?

Code

f = x -> ( sqrt(1 + x) - 1 ) / x
g = x-> 1/( 1 + sqrt(1 + x) )

t = [1e-14; 1e-15; 1e-16];
   
println( "                         f(x)                  g(x)" )
println( "x = 10^(-14):     ", f(t[1]), "    ",    g(t[1]) )
println( "x = 10^(-15):     ", f(t[2]), "   ",    g(t[2]) )
println( "x = 10^(-16):     ", f(t[3]), "                   ",    g(t[3]) )

                         f(x)                  g(x)
x = 10^(-14):     0.4884981308350689    0.4999999999999988
x = 10^(-15):     0.44408920985006256   0.4999999999999999
x = 10^(-16):     0.0                   0.5

Subtractive cancellation! The formula \(f\) suffers from subtractive cancellation, but the formula \(g\) does not.

(The formulae are the same since

\[\begin{align} f(x) &= \frac{\sqrt{1+x}-1}{x} \frac{\sqrt{1+x}+1}{\sqrt{1+x}+1} \nonumber\\ % &= \frac{1+x-1}{\sqrt{1+x}+1} = g(x). \nonumber \end{align}\]

Moreover, we can see that \(\lim_{x\to 0} f(x) = \lim_{x\to 0} g(x) = \frac12\).)

Quiz 0

Quiz 0: Exercise 1

Let \(x_n = \frac{1 + n^2}{ 2 n^3 } \cos n\). For which of the following sequences \((\varepsilon_n)\) do we have \(x_n = O(\varepsilon_n)\) as \(n\to\infty\)? Select all that apply.

(a). \(\varepsilon_n = e^{-n}\), (b). \(\varepsilon_n = n^{-2}\), (c). \(\varepsilon_n = n^{-1}\), (d). \(\varepsilon_n = 1\), (e). \(\varepsilon_n = n\).

(c), (d), and (e). Since \(|\cos n| \leq 1\) and \(1 \leq n^2\) for all \(n \in \mathbb N\), we have

\[\begin{align} |x_n| &= \left|\frac{1 + n^2}{ 2 n^3 } \cos n\right| \nonumber\\ % &\leq \frac{n^2 + n^2}{ 2 n^3 } \nonumber\\ % &= n^{-1} \leq 1 \leq n. \nonumber \end{align}\]

Therefore, \(x_n = O(n^{-1})\), \(O(1)\), and \(O(n)\).

However, \(x_n \not= O(n^{-2})\) (and so \(x_n \not= O(e^{-n})\) either): \(|\cos x| < \frac12\) only when \(x\) lies in an interval of the form \((\frac\pi3, \frac{2\pi}3) + k\pi\) (of length \(\frac\pi3 < 2\)), so at least one of any three consecutive integers \(n\) satisfies \(|\cos n| \geq \frac12\). For these \(n\), we have \(|x_n| \geq \frac{n^2}{2n^3} \cdot \frac12 = \frac1{4n}\), and so \(n^2 |x_n| \geq \frac{n}{4}\) is unbounded.

Quiz 0: Exercise 2

Let us define

\[\begin{align*} f(x) = \frac{\sqrt{1+x}-1}{x}, \qquad g(x) = \frac{1}{\sqrt{1+x} + 1}. \end{align*}\]

One may verify that \(f\) and \(g\) are mathematically identical (but you don’t have to). The following table is the output of evaluating \(f\) and \(g\) on \(\{10^{-14}, 10^{-15}, 10^{-16}\}\) in double precision floating point. Why are the function values different?

Code

f = x -> ( sqrt(1 + x) - 1 ) / x
g = x-> 1/( 1 + sqrt(1 + x) )

t = [1e-14; 1e-15; 1e-16];

println( "                         f(x)                  g(x)" )
println( "x = 10^(-14):     ", f(t[1]), "    ",    g(t[1]) )
println( "x = 10^(-15):     ", f(t[2]), "   ",    g(t[2]) )
println( "x = 10^(-16):     ", f(t[3]), "                   ",    g(t[3]) )

                         f(x)                  g(x)
x = 10^(-14):     0.4884981308350689    0.4999999999999988
x = 10^(-15):     0.44408920985006256   0.4999999999999999
x = 10^(-16):     0.0                   0.5

Subtractive cancellation! For small \(x\), \(\sqrt{1+x} \approx 1 + \frac{x}{2}\), so the numerator of \(f\) is the difference of two numbers that agree in almost all of their digits: the result has very few correct significant digits. At \(x = 10^{-16} < \frac12 \epsilon_{\rm mach}\) we even have 1 + 1e-16 == 1, so the numerator is exactly \(0\). The formula \(g\) only adds positive numbers, so there is no cancellation and \(g(x) \approx \frac12\) is accurate.

(The formulae are the same since

\[\begin{align} f(x) &= \frac{\sqrt{1+x}-1}{x} \frac{\sqrt{1+x}+1}{\sqrt{1+x}+1} \nonumber\\ % &= \frac{1+x-1}{x(\sqrt{1+x}+1)} = g(x). \nonumber \end{align}\]

Moreover, \(\lim_{x\to 0} f(x) = \lim_{x\to 0} g(x) = \frac12\).)

Lecture 3

Exercise 3.1

Let \(\epsilon\) be the machine precision. Is \((1 + \epsilon/2) - 1 = 1 + (\epsilon/2-1)\) in Float64?

Code

ϵ = eps( Float64 )
(1.0 + ϵ/2) - 1.0 == 1.0 + (ϵ/2 - 1.0)

false

This is false, since \(1 + \epsilon\) is the smallest float larger than \(1\), Julia has to round or truncate \(1 + \epsilon/2\):

Code

1 + ϵ/2

1.0

Therefore, \((1 + \epsilon/2) -1 = 0\). On the other hand \(\epsilon/2 - 1\) is a float and so \(1 + (\epsilon/2 - 1)\) is not equal to \(0\).

Exercise 3.2

Let \(X\) be an array containing \(10^{16}, -10^{16}\), and \(N\) copies of \(1\). What is the sum of all the entries of the array \(X\)? Does it depend on the order of summation? What about in Float64? Why?

Code

using Random

N = 10^2
X = vcat( [1e16, -1e16], fill(1., N) )

for _ ∈ 1:10 
    shuffle!(X)
    s = 0. 
    for x ∈ X
        s += x 
    end
    println(s)
end 

72.0
32.0
92.0
25.0
44.0
94.0
77.0
96.0
55.0
26.0

This happens because \(2^{53} \leq 10^{16} < 2^{54}\) and so the spacing between floats in this interval is \(2^{n-k} = 2^{53-52} = 2\) (in Float64, the binary precision is \(k=52\)). As a result \(10^{16} + 1\) is not a float; Julia rounds this to \(10^{16}\):

Code

1e16 + 1

1.0e16

Exercise 3.3

Exercise 1.4.1 from Driscoll and Braun (2017), [Fundamentals of Numerical Computation by Driscoll, Braun]

Code

f = x -> (1 - cos(x) )/sin(x)
g = x -> 2*sin(x/2)^2/sin(x)

(f(1e-6), g(1e-6))

(5.000444502912538e-7, 5.000000000000416e-7)

Since \(1 - \cos x = 2 \sin^2 \frac{x}{2}\), the formulas \(f(x) = \frac{1 - \cos x}{\sin x}\) and \(g(x) = \frac{2 \sin^2(x/2)}{\sin x}\) are mathematically identical (both equal \(\tan \frac{x}{2}\)).

In floating point, however, they are different for small \(x\). For \(x = 10^{-6}\) (see the code above):

• \(f(10^{-6}) \approx 5.000444 \times 10^{-7}\), with relative error \(\approx 9 \times 10^{-5}\) (only about \(4\) correct significant digits),
• \(g(10^{-6}) \approx 5.000000000000416 \times 10^{-7}\), which is correct to machine precision (relative error \(\approx 10^{-16}\)).

The problem itself is well-conditioned: with \(F(x) = \tan\frac{x}{2}\), we have \(\kappa_F(x) = \left|\frac{x F'(x)}{F(x)}\right| = \frac{x}{\sin x} \approx 1\) for small \(x\). The large error in \(f\) is therefore due to the algorithm: \(\cos(10^{-6}) \approx 1 - 5 \times 10^{-13}\), so computing \(1 - \cos x\) subtracts two nearly equal numbers (subtractive cancellation) and loses about \(12\) significant digits. The formula \(g\) involves no such subtraction, and so is stable.

Exercise 3.4

Exercise 1.4.2 from Driscoll and Braun (2017), [Fundamentals of Numerical Computation by Driscoll, Braun]

Code

x = @. 10.0^(-(1:11));
f = @. (exp(x) - 1)/x;

p=1;
for n ∈ 1:8
    p = @. p + x^n/( factorial(n+1) );
end 

abs.(f-p) ./ x

11-element Vector{Float64}:
    1.1102230246251565e-14
    1.0880185641326532e-12
    4.285460875053104e-11
    4.32542890393961e-9
    9.701794922989393e-7
    3.798317216308078e-5
    0.005663214341922184
    1.1077470940534795
   82.24037095772019
  826.9037099495335
 8273.537099867668

Code

# upper bound on the error between f and p as functions on x
err = @. (1/factorial(10)) * (x^9/(1-x))

11-element Vector{Float64}:
 3.061924358220656e-16
 2.7835675983824175e-25
 2.7584904128114e-34
 2.7560075231509052e-43
 2.7557594799933906e-52
 2.755734678133266e-61
 2.7557321979718108e-70
 2.755731949955909e-79
 2.755731925154322e-88
 2.7557319226741635e-97
 2.755731922426145e-106

Using the Taylor expansion of \(e^x\), we have

\[\begin{align} f(x) = \frac{e^x - 1}{x} = \sum_{n=0}^\infty \frac{x^n}{(n+1)!}, \qquad p(x) := \sum_{n=0}^{8} \frac{x^n}{(n+1)!}. \nonumber \end{align}\]

For \(0 < x < 1\), the truncation error satisfies

\[\begin{align} |f(x) - p(x)| = \sum_{n = 9}^\infty \frac{x^n}{(n+1)!} \leq \frac{1}{10!} \sum_{n=9}^\infty x^n = \frac{x^9}{10!\,(1-x)}, \nonumber \end{align}\]

which is less than \(3.1 \times 10^{-16}\) for \(x \leq 0.1\) (see the second code cell above). Since \(f(x) \geq 1\), the polynomial \(p\) is accurate to machine precision for all \(x \in (0, 0.1]\).

On the other hand, the formula \(\frac{e^x - 1}{x}\) suffers from subtractive cancellation for small \(x\) (\(e^x \approx 1\)): its relative error grows roughly like \(\epsilon_{\rm mach}/x\), from \(\approx 10^{-14}\) at \(x = 10^{-2}\) to \(\approx 10^{-8}\) at \(x = 10^{-8}\). The differences \(|f - p|\) computed above are therefore almost entirely the error of the naive formula, not the truncation error of \(p\). (Julia’s expm1 avoids this problem: expm1(x)/x is accurate for all small \(x\).)

Exercise 3.5

Let \(f(x) = e^x - 2x - 1\). Apply change of sign theorem to \(f\) to conclude there exists \(\xi \in [1,2]\) such that \(f(\xi) = 0\).

The function \(f(x) = e^x - 2x - 1\) is continuous with \(f(1) = e - 3 \approx -0.28 < 0\) and \(f(2) = e^2 - 5 \approx 2.39 > 0\). Therefore, \(f(1) f(2) < 0\) and, by the change of sign theorem, there exists \(\xi \in [1,2]\) with \(f(\xi) = 0\).

Exercise 3.6

Let’s return to Exercise 3.5. You showed that \(f(x) = e^x - 2x - 1\) has a root \(\xi \in [1,2]\).

Show that \(f(x) = 0\) can be solved by instead considering the fixed point problems \(g(x) = x\) or \(h(x) = x\) where

\[ \begin{align} g(x) &= \log( 2x + 1 ), \\ h(x) &= \frac{e^x - 1}{2}. \end{align} \]

Show that Brouwer’s fixed point theorem applies to \(g\) but not \(h\) on \([1,2]\).

For \(x \in [1,2]\) we have \(2x + 1 > 0\), and \(\log\) is one-to-one, so

\[\begin{align} g(x) = x &\iff \log(2x+1) = x \nonumber\\ &\iff 2x + 1 = e^x \nonumber\\ &\iff f(x) = 0, \nonumber \end{align}\]

and

\[\begin{align} h(x) = x &\iff \frac{e^x - 1}{2} = x \nonumber\\ &\iff e^x - 2x - 1 = 0 \nonumber\\ &\iff f(x) = 0. \nonumber \end{align}\]

That is, the roots of \(f\) are exactly the fixed points of \(g\) and of \(h\).

\(g\): \(g\) is continuous and increasing on \([1,2]\), with \(g(1) = \log 3 \approx 1.10\) and \(g(2) = \log 5 \approx 1.61\). So, for all \(x \in [1,2]\),

\[\begin{align} 1 < g(1) \leq g(x) \leq g(2) < 2. \nonumber \end{align}\]

Therefore \(g : [1,2] \to [1,2]\) and Brouwer’s fixed point theorem gives a fixed point \(\xi = g(\xi) \in [1,2]\).

\(h\): \(h(1) = \frac{e - 1}{2} \approx 0.86 < 1\) (and \(h(2) = \frac{e^2 - 1}{2} \approx 3.19 > 2\)), so \(h\) does not map \([1,2]\) into \([1,2]\) and Brouwer’s fixed point theorem does not apply. (Of course, \(h\) still has the fixed point \(\xi \in [1,2]\) found above: the theorem just cannot be used to show it.)

Exercise 3.7

Assignment 1, Fall 2025, [30 pts]. Prove that, for \(x \geq 0\), there exists a unique \(W(x) \geq 0\) for which \(W(x) e^{W(x)} = x\). Show that \(\Omega := W(1)\) satisfies \(0.56 < \Omega < 0.58\).

Intermediate value theorem is useful!

The function \(F(w) := w e^w\) is continuous on \([0, \infty)\) with \(F(0) = 0\) and \(F(w) \to \infty\) as \(w \to \infty\). Therefore, for \(x \geq 0\), there exists \(w_x \geq 0\) with \(F(0) \leq x \leq F(w_x)\) and, by the intermediate value theorem, there exists \(W \in [0, w_x]\) with \(F(W) = x\). Moreover, \(F'(w) = (1 + w) e^w > 0\), so \(F\) is strictly increasing on \([0,\infty)\): if \(W' < W < W''\) then \(F(W') < x < F(W'')\). Hence \(W(x) := W\) is the unique solution in \([0,\infty)\) of \(W e^W = x\).

Finally, \(F(0.56) = 0.56\, e^{0.56} \approx 0.980 < 1\) and \(F(0.58) = 0.58\, e^{0.58} \approx 1.035 > 1\), and so (again since \(F\) is continuous and strictly increasing) \(\Omega = W(1) \in (0.56, 0.58)\). In the plot below, \(\Omega\) is where the curve \(y = w e^w\) crosses the line \(y = 1\):

Code

using Plots, LaTeXStrings
F = w -> w * exp(w)

plot( F, 0.3, 0.75; label = L"y = w e^w", xlabel = L"w", lw = 2, legend = :topleft )
hline!( [1]; label = L"y = 1", linestyle = :dash )
scatter!( [0.56, 0.58], F.([0.56, 0.58]); label = L"w = 0.56,\ 0.58", color = :red )

Exercise 3.8

Assignment 2, Fall 2025, [10 pts]. For fixed \(f\), define

\[\begin{align} g(x) &= f(x - a), \\ h(x) &= b f(x), \\ j(x) &= f(c x), \end{align}\]

1. What is \(\kappa_g(x)\)?
2. Write \(\kappa_{h}(x)\) in terms of \(\kappa_f\),
3. Write \(\kappa_{j}(x)\) in terms of \(\kappa_f\),
4. What is \(\kappa_f(x)\) when \(f(x) = e^{-x^2}\)?

Recall that \(\kappa_f(x) = \left| \frac{x f'(x)}{f(x)} \right|\).

1. \(g'(x) = f'(x - a)\) and so \(\kappa_g(x) = \left| \frac{x f'(x-a)}{f(x-a)} \right| = \left| \frac{x}{x - a} \right| \kappa_f(x - a)\).
2. \(h'(x) = b f'(x)\) and so \(\kappa_h(x) = \left| \frac{x\, b f'(x)}{b f(x)} \right| = \kappa_f(x)\) (for \(b \not= 0\)): scaling the output doesn’t change the conditioning.
3. \(j'(x) = c f'(cx)\) and so \(\kappa_j(x) = \left| \frac{x\, c f'(cx)}{f(cx)} \right| = \kappa_f(cx)\). (Common mistakes: \(|c| \kappa_f(cx)\), \(\frac{1}{|c|}\kappa_f(x)\), or \(\frac{1}{|c|}\kappa_f(cx)\).)
4. \(f'(x) = -2x e^{-x^2}\) and so \(\kappa_f(x) = \left| \frac{x \cdot (-2x) e^{-x^2}}{e^{-x^2}} \right| = 2x^2\).

Exercise 3.9

Assignment 2, Fall 2025, [25 pts]. Consider the differential equation (initial value problem)

\[ \begin{align} \begin{split} u'(x) &= u(x) +e^{-x} \nonumber \\ u(0) &= -\frac12 . \end{split} \end{align} \]

1. Show that \(u(x) = -\frac12 e^{-x}\) solves (3.2).

Now we perturb the initial data, by considering

\[ \begin{align} \begin{split} u'(x) &= u(x) +e^{-x} \nonumber \\ u(0) &= \alpha \end{split} \end{align} \]

for \(\alpha := \epsilon - \frac12\).

2. Show that \(u_{\alpha}(x) = (\alpha + \frac12) e^x -\frac12 e^{-x}\) solves (3.3).
3. Plot the functions \(u\) and \(u_{\alpha}\) and comment on the difference between the solutions for large \(x\).
4. Let \(f( \alpha ) = u_\alpha( 100 )\). What is the condition number \(\kappa_f( \alpha )\)?
5. Explain what happens for \(\alpha \approx -\frac12\). Link your observations to your answer for question 3.

1. \(u'(x) = \frac12 e^{-x} = -\frac12 e^{-x} + e^{-x} = u(x) + e^{-x}\) and \(u(0) = -\frac12\).

2. \(u_\alpha'(x) = (\alpha + \frac12) e^x + \frac12 e^{-x} = u_\alpha(x) + e^{-x}\) and \(u_\alpha(0) = (\alpha + \frac12) - \frac12 = \alpha\).

3. The difference between the two solutions is \(u_\alpha(x) - u(x) = (\alpha + \frac12) e^{x} = \epsilon e^x\). Here are \(u\) and \(u_\alpha\) for \(\epsilon = 10^{-10}\):

Code

α = 1e-10 - 1/2;

u = x -> -exp(-x)/2;
uα = x -> (α + 1/2)*exp(x) + u(x);

plt1 = plot([u, uα], 0, 10, label=[L"u" L"u_\alpha"], title="Small x")
plt2 = plot([u, uα], 0, 100, label=[L"u" L"u_\alpha"], title="Large x")

plot( plt1, plt2, size=(1000,500))

For small \(x\) (left) the two solutions are indistinguishable: \(u_\alpha(x) \approx u(x)\). For large \(x\) (right) they are completely different: \(u(x) = -\frac12 e^{-x}\) decays to \(0\), whereas \(u_\alpha(x) \approx \epsilon e^x\) grows exponentially. For example, \(u(20) \approx -10^{-9}\) but \(u_\alpha(20) \approx 0.05\), and \(u(30) \approx -5 \times 10^{-14}\) but \(u_\alpha(30) \approx 1069\).

The difference is easiest to see on a logarithmic scale:

Code

err = x -> abs( u(x) - uα(x) ) 

plot( err, 0, 100 , yaxis=:log, legend = false, title="absolute error between u and uα" )

On a logarithmic \(y\)-axis, the error \(|u(x) - u_\alpha(x)| = \epsilon e^x\) is a straight line with slope \(1\) (in natural-log units): a perturbation of size \(\epsilon = 10^{-10}\) in the initial condition grows by a factor \(e\) for every unit increase in \(x\), and is about \(3 \times 10^{33}\) at \(x = 100\).

4. \(f(\alpha) = u_\alpha(100) = (\alpha + \frac12) e^{100} - \frac12 e^{-100}\) and \(f'(\alpha) = e^{100}\), so \[\begin{align} \kappa_f(\alpha) = \left| \frac{\alpha f'(\alpha)}{f(\alpha)} \right| = \left| \frac{\alpha e^{100}}{(\alpha + \frac12) e^{100} - \frac12 e^{-100}} \right|. \nonumber \end{align}\]

5. As \(\alpha \to -\frac12\), \(\kappa_f(\alpha) \to \left| \frac{-\frac12 e^{100}}{-\frac12 e^{-100}} \right| = e^{200} \approx 7 \times 10^{86}\): the problem is extremely ill-conditioned. A tiny relative change in the initial condition (e.g. storing \(-\frac12 + \epsilon\) in floating point) leads to a huge relative change in \(u(100)\), which is exactly what the plots in part 3 show.

   For example, with \(\epsilon = 10^{-10}\) as in the plots: the relative change in the input is \(\left|\frac{\epsilon}{-1/2}\right| = 2 \times 10^{-10}\), while \(u(100) \approx -1.9 \times 10^{-44}\) and \(u_\alpha(100) \approx 2.7 \times 10^{33}\), a relative change of about \(1.4 \times 10^{77}\) in the output. Their ratio is \(\approx 7.2 \times 10^{86} \approx e^{200} \approx \kappa_f(\alpha)\), as predicted.

Exercise 3.10

Assignment 2, Fall 2025, [10 pts]. Let \(p_{\alpha}(x) = \alpha x^n + \sum_{j=0}^{n-1} a_{j} x^j\) be a polynomial of degree \(n\) and let \(r(\alpha)\) be a root of \(p_{\alpha}\). What is the condition number \(\kappa_r(\alpha)\)?

Since \(p_\alpha(r(\alpha)) = 0\) for all \(\alpha\), differentiating with respect to \(\alpha\) gives

\[\begin{align} r(\alpha)^n + p_\alpha'\big(r(\alpha)\big)\, r'(\alpha) = 0 \qquad \text{and so} \qquad r'(\alpha) = - \frac{r(\alpha)^n}{p_\alpha'\big(r(\alpha)\big)} \nonumber \end{align}\]

(assuming \(r(\alpha)\) is a simple root, so that \(p_\alpha'(r(\alpha)) \not= 0\)). Therefore,

\[\begin{align} \kappa_r(\alpha) = \left| \frac{\alpha\, r'(\alpha)}{r(\alpha)} \right| = \left| \frac{\alpha\, r(\alpha)^{n-1}}{p_\alpha'\big(r(\alpha)\big)} \right|. \nonumber \end{align}\]

In particular, the root is ill-conditioned when \(p_\alpha'(r(\alpha))\) is small, i.e. when \(r(\alpha)\) is close to a repeated root.

Lecture 4

Exercise 4.1

Suppose \(g :[a,b] \to \mathbb R\) is continuously differentiable with \(|g'(x)| \leq L < 1\) for all \(x \in [a,b]\). Show that \(g\) is a contraction with constant \(L\).

Use Taylor/mean value theorem.

By the mean value theorem, for \(x, y \in [a,b]\), there exists \(c \in [a,b]\) such that

\[\begin{align} |g(x) - g(y)|&= |g'(c)| |x-y| \nonumber\\ % &\leq L |x-y|. \end{align}\]

As a result, \(g\) is a contraction.

Exercise 4.2

Take the simple iteration \(x_{n+1} = g(x_n)\) for some contraction \(g : [a,b] \to [a,b]\). Show that \((x_n)\) converges to the fixed point of \(g\) in \([a,b]\) at least linearly with asymptotic error constant \(\mu = |g'(\xi)|\).

Notice that, by the mean value theorem, there exists \(\eta_n\) between \(x_{n}\) and \(\xi\) for which

\[\begin{align} \frac {|x_{n+1} - \xi|} {|x_n - \xi|} % &= \frac {|g(x_n) - g(\xi)|} {|x_n - \xi|} \nonumber\\ % &= \frac {|g'(\eta_n)| |x_n - \xi|} {|x_n - \xi|} \nonumber\\ % &= |g'(\eta_n)| \nonumber \end{align}\]

Following the hint: since \((x_n) \to \xi\), and \(\eta_n\) is between \(x_n\) and \(\xi\), we have

\[\begin{align} \lim_{n\to\infty} \frac {|x_{n+1} - \xi|} {|x_n - \xi|} % &= |g'(\xi)| \nonumber \end{align}\]

and the order of convergence is at least linear with asymptotic error constant \(\mu = |g'(\xi)|\).

(Recal that, if \(\mu = 0\), the order of convergence is superlinear)

Exercise 4.3

Show that \(g(x) = \log( 2x + 1)\) is a contraction on \([1,2]\) with \(L = \frac23\). Hence conclude that \(x_{n+1} := g(x_n)\) converges for all \(x_1 \in [1,2]\).

Notice that \(g'(x) = \frac{2}{2x+1}\) which is positive and decreasing on \([1,2]\). We therefore have \(\frac25 = g'(2) \leq g'(x) \leq g'(1) = \frac{2}{3}\) and so \(g\) is a contraction with \(L = \frac23\) by Exercise 4.1. Applying the contraction mapping theorem shows that \((x_n)\) converges linearly to the fixed point \(\xi \in [1,2]\) of \(g\).

Exercise 4.4

Assignment 2, Fall 2025, [25 pts]. Fix \(a, x_1 > 0\) and define the iteration \(x_{n+1} = \frac{1}{2}\big( x_n + \frac{a}{x_n} \big)\).

1. ✍ Show that \(x_{n+1}\) has the same sign as \(x_n\),
2. ✍ Suppose that \((x_n) \to \xi\). What is the value of \(\xi\)?
3. 💻 Let \(x\) be as in question 2. Show numerically that \((x_n) \to \xi\) for a few different choices of \(a\) and \(x_1\),
4. 💻 This iteration can be seen as a way of approximating \(\xi\). What is the correseponding order of convergence?

1. Note \(\frac12 x_n\) and \(\frac1{2 x_n}\) have the same sign as \(x_n\). Moreover, since \(a> 0\), \(\frac{a}{x_n}\) has the same sign as \(x_n\). As a result \(x_{n+1}\) has the same sign as \(x_n\). Since \(x_1 > 0\), we find that \(x_n > 0\) for all \(n\).

2. Since \(g(x) := \frac12\big( x + \frac{a}{x} \big)\) is continuous on \((0,\infty)\) and \((x_n)\to\xi\), we have \(\xi = g(\xi)\). Therefore, we have \(\frac12\xi = \frac12 \frac{a}{\xi}\) and so \(\xi^2 = a\). Since \(\xi > 0\) since \(x_n > 0\) for all \(n\), we have \(\xi = + \sqrt{a}\).

Code

A = [2, 4, 9, 81]
X = [1., 2., 1e-2]
for a ∈ A
    for x1 ∈ X 
        g = x -> (1/2) * (x + a/x)
        print( "(a, x₁) = (", a, ", ", x1, "); " )

        x, boo = simple_iteration( g, x1 )
        if boo 
            println( "converged to ", x[end] )
        else
            println( " [didn't converge]" )
        end
    end
    println("")
end

(a, x₁) = (2, 1.0); converged to 1.414213562373095
(a, x₁) = (2, 2.0); converged to 1.414213562373095
(a, x₁) = (2, 0.01); converged to 1.414213562373095

(a, x₁) = (4, 1.0); converged to 2.0
(a, x₁) = (4, 2.0); converged to 2.0
(a, x₁) = (4, 0.01); converged to 2.0

(a, x₁) = (9, 1.0); converged to 3.0
(a, x₁) = (9, 2.0); converged to 3.0
(a, x₁) = (9, 0.01); converged to 3.0

(a, x₁) = (81, 1.0); converged to 9.0
(a, x₁) = (81, 2.0); converged to 9.0
(a, x₁) = (81, 0.01); converged to 9.0

This method converges to \(\sqrt{a}\) for all \(a\) and \(x_1\) that we chose.

4. Here, we compute \(x_{n+1} = g(x_n)\) for \(a = 2\) and \(x_1 = 1\):

Code

function μ( x, ξ; α=1 )
    return @. abs( x[2:end] - ξ ) / ( abs(x[1:end-1] - ξ )^α );
end 

g = x -> (1/2) * (x + 2/x)
x, _ = simple_iteration(g, 1; tol=1e-5)

x

5-element Vector{Float64}:
 1.0
 1.5
 1.4166666666666665
 1.4142156862745097
 1.4142135623746899

Next, we compute \(\frac{|x_{n+1} - \xi|}{|x_n-\xi|^\alpha}\) for \(\alpha = 1\):

Code

display( μ( x, √2 ) ) 

4-element Vector{Float64}:
 0.20710678118654724
 0.028595479208965498
 0.0008658015152821018
 7.508466926336308e-7

This appears to converge to \(0\); suggesting the convergence is superlinear.

For \(\alpha = 2\), we have:

Code

display( μ( x, √2; α=2) ) 

4-element Vector{Float64}:
 0.49999999999999917
 0.33333333333330084
 0.35294117643144196
 0.35352238456116425

This is a number \(>0\), suggesting the convergence is quadratic.

In fact, in the Lecture 5, we will prove that the convergence is quadratic with \(\mu = \frac1{2\sqrt{a}}\). In the last case, this is

Code

1/( 2√2 )

0.35355339059327373

which matches the numerically computed value above.

We could have also numerically approximated \(\alpha\) by computing \(\frac{\log |x_{n+1} - \xi|}{\log |x_n - \xi|}\):

Code

@. log(abs(x[2:end]-√2)) / log(abs(x[1:end-1]-√2))

4-element Vector{Float64}:
 2.7864397013573967
 2.4473369494135158
 2.1732752731987097
 2.0796040502714805

Exercise 4.5

Assignment 3, Fall 2025. The Lambert W function (also called the product-log function) \(W : [0,\infty) \to [0,\infty)\) satisfies \(W(a) e^{W(a)} = a\). Define \(\Omega := W(1)\).

1. ✍ Show that \(\Omega \in [0.56, 0.58]\).
2. ✍ Show that \(\Omega\) is a fixed point of \(g\) and \(h\) as defined by

\[\begin{align} g(x) = e^{-x}, \qquad h(x) = \log \left|\frac{1}{x}\right| \end{align}\]

3. 💻 Show numerically that \(x_{n+1} = g(x_n)\) converges and \(y_{n+1} = h(y_n)\) diverges for a few different choices of \(x_1\) and \(y_1\).
4. ✍ Show that \(g:[\frac{1}{e}, 1]\to[\frac{1}{e},1]\) with \(|g'| < 1\),
5. ✍ Explain why this means that the iteration \(x_{n+1} = g(x_n)\) converges for all \(x_1 \in [\frac{1}{e}, 1]\),
6. ✍ Now fix \(x_1 \in \mathbb R\). Show that the sequence \((x_n)\) is eventually contained in \([0,1]\) and that \(g : [0,1] \to [\frac{1}{e}, 1]\). Explain why this means that \((x_n)\) converges for all \(x_1 \in \mathbb R\),
7. ✍ What is the order of convergence of this method?
8. 💻 Confirm this numerically.
9. ✍ Show that \((y_n)\) doesn’t converge to \(\Omega\) for any choice of \(y_1 \not= \Omega\).

Solution

Exercise 4.6

Midterm, Fall 2025, Optional section, 30 pts. Define \(g(x) = 1 + \frac{1}{1 + x}\) and \(I := [1,2]\).

1. [3 points] Show that \(g(x) \in I\) for all \(x \in I\).
2. [3 points] Which theorem guarantees the existence of a fixed point \(\xi = g(\xi) \in I\).
3. [3 points] Compute the value of \(\xi\).
4. [3 points] Show that \(L := \max_{x \in I} |g'(x)| = \frac14\).
5. [3 points] Explain why this means the iteration \(x_{n+1} = g(x_n)\) converges to \(\xi\) for all \(x_1 \in I\).
6. [5 points] What is the order of convergence?
7. [5 points] What is the asymptotic error constant?
8. [5 points] Fix \(x_1 \in I\). Find \(n\) for which the error bound \(|x_{n+1} - \xi| \leq \frac{1}{64}\) is guaranteed.

Solution

Quiz 1

Quiz 1: Exercise 1

Let \(x_n = \frac{1 + n^2}{ 3 n^4 } \cos n^5\). For which of the following sequences \((\varepsilon_n)\) do we have \(x_n = O(\varepsilon_n)\) as \(n\to\infty\)? Select all that apply.

• \(\varepsilon_n = e^{-n}\),
  
• \(\varepsilon_n = n^{-2}\),
  
• \(\varepsilon_n = n^{-1}\),
  
• \(\varepsilon_n = 1\),
  
• \(\varepsilon_n = n\).

(b), (c), (d), and (e). Since \(|\cos n^5| \leq 1\) and \(1 \leq n^2\) for all \(n \in \mathbb N\), we have

\[\begin{align} |x_n| &= \left|\frac{1 + n^2}{ 3 n^4 } \cos n^5\right| \nonumber\\ % &\leq \frac{n^2 + n^2}{ 3 n^4 } \nonumber\\ % &= \frac23 n^{-2} \leq n^{-1} \leq 1 \leq n. \nonumber \end{align}\]

Therefore, \(x_n = O(n^{-2})\), \(O(n^{-1})\), \(O(1)\), and \(O(n)\). However, \(x_n \not= O(e^{-n})\): e.g. along the subsequence of \(n\) for which \(|\cos n^5| \geq \frac12\) (there are infinitely many), \(|x_n| \geq \frac16 n^{-2}\), which decays much slower than \(e^{-n}\).

Quiz 1: Exercise 2

Note that \(x := 1.001^{2026}\) and \(y := e^{2026 \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):

Code

b, n, x = Float16( 1.001 ), 2026, Float16( 1 );
for i ∈ 1:n
    x = b*x
end
y = exp( n * log( b ) )
println( "x = ", x, "; y = ", y )

x = 5.914; y = 7.227

Which result do you trust more and why?

Trust \(y\) more. Computing \(x\) by repeated multiplication involves \(2026\) multiplications, each of which introduces a round-off error (of relative size up to \(\frac12 \epsilon_{\rm mach} = 2^{-11}\) in Float16) and these errors accumulate. On the other hand, \(y\) is computed with only three operations (log, *, and exp), so there are far fewer opportunities for round-off error.

(In fact, neither value equals \(1.001^{2026} \approx 7.576\): Float16(1.001) is \(1 + 2^{-10} \approx 1.000977\), and \(y \approx 7.227\) is an accurate value of \((1 + 2^{-10})^{2026} \approx 7.225\).)

Quiz 1: Exercise 3

True or false: there exists a continuous function \(f\) with a root in \([a,b]\) and \(f(a) f(b) > 0\)? Why or why not?

True. For example, \(f(x) = x^2\) has a root in \([-1,1]\) but \(f(-1) f(1) = 1 > 0\). That is, the converse of the change of sign theorem is false.

Quiz 1: Exercise 4

True or false: the function \(f(x) = x^2\) is a contraction on some interval \([a,b]\). Show your working.

True. For example, on \([0, \frac14]\) we have \(|x^2 - y^2| = |x+y|\,|x-y| \leq \frac12 |x-y|\), so \(f\) is a contraction with \(L = \frac12\).

More generally, \(|x^2 - y^2| = |x+y|\,|x-y| \leq 2\max(|a|,|b|)\,|x-y|\) for \(x, y \in [a,b]\) (or use the mean value theorem: \(|x^2 - y^2| = 2|\xi|\,|x-y|\) for some \(\xi\) between \(x\) and \(y\)). Therefore, \(f\) is a contraction on any \([a,b] \subset (-\frac12, \frac12)\) with \(L = 2\max(|a|,|b|) < 1\).

Lecture 5

Exercise 5.1

Suppose \(g : [a,b] \to [a,b]\) is a contraction with constant \(L < 1\) and consider the simple iteration \(x_{n+1} = g(x_n)\). Let \(\xi \in [a,b]\) be the unique fixed point of \(g\) in \([a,b]\).

1. Show that \(|x_1 - \xi| \leq \frac1{1-L} |x_1 - x_2|\),
2. Hence, write down a lower bound on the number of iterations \(n\) to ensure \(|x_{n+1} - \xi|\leq \varepsilon\).

1. Use the triangle inequaltiy,
2. Use the fact that

\[\begin{align} |x_{n+1} - \xi| \leq \frac{L^n}{1-L} |x_1 - x_2| \end{align}\]

Recall that \(|x_{n+1} - \xi| \leq L |x_{n} - \xi|\). In particular, we have \(|x_1 - \xi| \leq |x_1 - x_2| + |x_2 - \xi| \leq |x_1 - x_2| + L|x_1-\xi|\) and thus \(|x_1 - \xi| \leq \frac1{1-L} |x_1 - x_2|\). As a result, we have

\[\begin{align} |x_{n+1} - \xi| &\leq L^n |x_1 - \xi| \nonumber\\ &\leq \frac {L^n} {1-L} |x_1 - x_2| \end{align}\]

This number is less than \(\epsilon\) when

\[\begin{align} n \geq \frac{ \log \frac{|x_1-x_2|}{\epsilon(1-L)} }{\log\frac1L}. \end{align}\]

This gives an upper bound on the number of iterations needed in order to be within some fixed tolerance.

Let’s back to the function \(g(x) = \log(2x+1)\). For an error of \(\epsilon = 10^{-10}\), we require around \(54\) iterations. Actually, \(40\) will do (but we have an upper bound):

Code

g(x) = log( 2x + 1 )
x, _ = simple_iteration( g, 1.0; N=1000, tol=0.0 )    # fixed point ξ ≈ 1.256 
ξ = x[end]

ϵ = 1e-10; 
L = 2/3;
x1 = 1; x2 = g(x1);
@show N = log( abs(x1-x2) / ( ϵ * (1-L) ) ) / log( 1/L )

x, _ = simple_iteration( g, 1.0; N=54, tol=1e-20 )

abs( x[40] - ξ )

N = log(abs(x1 - x2) / (ϵ * (1 - L))) / log(1 / L) = 53.78490871074893

8.923461969345681e-11

Exercise 5.2

Exercise 4.4 continued…

This is an example of Newton’s method!

\[\begin{align} x_{n+1} = \frac12 \left( x_n + \frac{a}{x_n} \right) \end{align}\]

is the Newton iteration with \(f(x) = x^2 - a\). Notice that \(\xi = \sqrt{a}\) and \(f'(\xi) = 2\sqrt{a}\), \(f''(\xi) = 2\).

Compute the order of convergence and asymptotic error constant numerically, and compare to Theorem 5.3.

Here, we plot the order of convergence and asymptotic error constant:

Code

g(x) = x^2 - 2
g′(x) = 2*x 
orderOfConvergence( Newton(g,g′,1.), sqrt(2); α=2 )

println( "theoretical value of μ (for α=2) = ",  (1/2) * ( 2/ g′(sqrt(2)) ) )

iteration  abs error        log error    α          μ (α=2.0)     
1          4.142136e-01     -0.3828      NaN        NaN           
2          8.578644e-02     -1.0666      2.7864     0.5000        
3          2.453104e-03     -2.6103      2.4473     0.3333        
4          2.123901e-06     -5.6729      2.1733     0.3529        
5          1.594724e-12     -11.7973     2.0796     0.3535        
theoretical value of μ (for α=2) = 0.35355339059327373

Exercise 5.3

Assignment 5, Fall 2025, Section B

What happens when applying Newton to Kepler’s equation with \(\epsilon = 1\) (in this case \(f(2\pi) = f'(2\pi) = f''(2\pi) = 0\)).

Exercise 5.4

Suppose that \(g\) is twice continuously differentiable with \(g(\xi) = \xi\) and \(g'(\xi) = 0\).

1. ✍️ Show that there exists \(\delta > 0\) such that \(x_{n+1} = g(x_n)\) converges to \(\xi\) for all \(x_1 \in [\xi - \delta, \xi + \delta]\).
   
2. ✍️ Show that

\[\begin{align} \lim_{n\to\infty} \frac{|x_{n+1} - \xi|}{|x_n - \xi|^2} = \frac12 |g''(\xi)| \nonumber \end{align}\]

(whenever \(x_n \not= \xi\) for all \(n\)). What is the order of convergence when \(g''(\xi) \not= 0\)?

Exercise 5.5

Newton’s method \(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\) is the simple iteration with \(g(x) := x - \frac{f(x)}{f'(x)}\).

✍️ Suppose that \(f\) is three times continuously differentiable with \(f(\xi) = 0\) and \(f'(\xi) \not= 0\). Show that \(g'(\xi) = 0\) and \(g''(\xi) = \frac{f''(\xi)}{f'(\xi)}\). What does this tell you about Newton’s method?

Taylor’s theorem.

Exercise 5.6

Fix \(a > 0\).

1. ✍️ Show that Newton’s method applied to \(f(x) := \frac1x - a\) is given by \(x_{n+1} = x_n(2 - a x_n)\).
   
2. ✍️ Define \(e_n := 1 - a x_n\) and show that \(e_{n+1} = e_n^2\). Hence, find all \(x_1\) for which \((x_n) \to \frac1a\).
   
3. ✍️ Show that \(\left| x_{n+1} - \frac1a \right| = a \left| x_n - \frac1a \right|^2\). What are the order of convergence and the asymptotic error constant? Compare with the asymptotic error constant \(\mu = \frac12 \left| \frac{f''(\xi)}{f'(\xi)} \right|\) from lectures.
   
4. 💻 Take \(a = 7\). Verify your answers to (ii) and (iii) numerically using \(x_1 = 0.1, 0.28\), and \(0.3\).

1. We have \(f(x) = \frac1x - a\) and \(f'(x) = -\frac1{x^2}\) and so Newton is

\[\begin{align} x_{n+1} &= x_n - \frac{f(x_n)}{f'(x_n)} \nonumber\\ % &= x_n - \frac{\frac1{x_n} - a}{ - \frac{1}{x_n^2} } \nonumber\\ % &= x_n( 2 - a x_n ) \end{align}\]

If \(x_n \to \xi\), then \(\xi = \frac1a\).

2. Notice \(e_n = 1 - a x_n\) and so

\[\begin{align} e_{n+1} &= 1 - a x_{n+1} \nonumber\\ &= 1 - a x_n (2 - ax_n) \nonumber\\ &= 1 - 2 a x_n + (ax_n)^2 \nonumber\\ &= (1 - a x_n)^2 = e_n^2. \end{align}\]

We have \(x_n \to \frac1a\) if and only if \(e_n \to 0\). Since \(e_{n+1} = e_n^2 = e_1^{2^n}\), we have \(e_n \to 0\) for all \(x_1\) with \(|e_n| = |1 - a x_1| < 1\). That is, if \(0 < x_1 < 2/a\).

3. We have

\[\begin{align} |x_{n+1} - \frac1a| &= \frac{1}{a} | 1 - a x_{n+1} | \nonumber\\ &= \frac1a e_{n+1} \nonumber\\ &= \frac1a |1 - a x_n|^2 \nonumber\\ &= a |x_n - \frac1a|^2. \end{align}\]

Therefore, the order of convergence is \(\alpha = 2\) with asymptotic error constant \(\mu = a\). Here, we are applying Newton’s method to \(f(x) = \frac1x - a\) and so \(f'(x) = -\frac1{x^2}\), \(f''(x) = \frac2{x^3}\), and \(\frac12 \big|\frac{f''(1/a)}{f'(1/a)}\big| = \frac12 \frac2{a^{-3}} (a^{-2}) = a\) (which agrees with the above).

Problem Class 1

Exercise 1

Show that each of the following have fixed points on \([2,4]\):

\[\begin{align} g_1(x) &= 1 + x - \frac19 x^2 \nonumber\\ g_2(x) &= \pi + 0.99 \sin x \nonumber\\ g_3(x) &= x + 1 - \tan\frac{x}4 \end{align}\]

1. Do you expect the simple iterations \(x_{n+1} = g_j( x_n )\) to converge?
2. Why?
3. Test your hypothesis numerically.
4. Compute the order of convergence and asymptotic error constants.

Exercise 2

Suppose that \(g(\xi) = \xi\), \(g'(\xi) = 0\) and \(|g''|\) is bounded near \(\xi\).

1. Show that the simple iteration \(x_{n+1} = g(x_n)\) with \(x_1 \approx \xi\) converges to \(\xi\).
2. What is the order of convergence?
3. Can you generalise this statement to the case where \(g'(\xi) = \cdots = g^{(\ell-1)}(\xi) = 0\) and \(g^{(\ell)}\) bounded near \(\xi\)?

Exercise 3

For \(\alpha \in (0,1]\), define \(x_{n+1} := x_n + \alpha ( g(x_n) - x_n )\) where \(g\) is continuous.

1. Show that if \((x_n) \to \xi\) then \(\xi\) is a fixed point of \(g\)
2. Suppose that \(g'\) is continuous with \(g'(\xi) < 0\). Show that there exists \(\alpha_0 > 0\) such that, if \(\alpha \leq \alpha_0\), then \((x_n)\) converges to \(\xi\) for all \(x_1\) sufficiently close to \(\xi\)

Lecture 6

Exercise 6.1

Show that the degree \(m_0\) Taylor polynomial of \(f\) about \(x_0\) solves (6.1) when \(X = \{ x_0 \}\).

Recall that

\[\begin{align} p(x) = f(x_0) + f'(x_0)(x-x_0) + \cdots + \frac{f^{(m_0)}(x_0)}{m_0!}(x-x_0)^{m_0} % \nonumber \end{align}\]

is the degree \(m_0\) Taylor polynomial of \(f\) about \(x_0\).

Exercise 6.2

Write down the polynomial of least degree interpolating \(f(x) = \sin(x)\) at the points \(0, \frac{\pi}{2}, \pi\)

Exercise 6.3

Write down an error estimate for approximating \(f(x) = \sin(x)\) by \(I_{\{0, \frac\pi2, \pi\}}f(x) = - \frac4{\pi^2} x (x - \pi)\).

Exercise 6.4

Consider \(f(x) = e^{-x}\) on \([-1,1]\).

1. Find a Lagrange interpolating polynomial \(p\) of \(f\) on \(\{ -1, 0 , 1\}\),

2. Write down an error estimate for approximating \(f\) with \(p\) on \([-1,1]\).

Exercise 6.5

(\(\star\) more difficult)

Let \(X = \{x_0 < \cdots < x_n\}\) be \(n+1\) equispaced nodes on \([-1,1]\), i.e. \(x_j = -1 + jh\) with \(h = 2/n\), and let \(\ell_X(x) = (x-x_0)(x-x_1)\cdots(x-x_n)\) be the node polynomial.

1. Fix \(x \in [x_k, x_{k+1}]\) for some \(0 \leq k \leq n-1\). By bounding \(|x - x_j|\) for \(j < k\) and \(j > k+1\) and using \(|x-x_k|\,|x-x_{k+1}| \leq h^2/4\), show that

   \[ |\ell_X(x)| \leq \frac{h^{n+1}}{4} (k+1)!\,(n-k)!, \]

   and deduce that \(\max_{x\in[-1,1]} |\ell_X(x)| \leq \dfrac{n!}{4} h^{n+1}\).

2. Use Stirling’s formula, \(n! \sim \sqrt{2\pi n}\,(n/e)^n\), to show that there is a constant \(C > 0\), independent of \(n\), such that

   \[ \max_{x\in[-1,1]} |\ell_X(x)| \leq C\, \frac{(2/e)^n}{\sqrt n} \qquad \text{for all } n. \]

3. By evaluating \(\ell_X\) at the point \(x^* = x_0 + h/2\), show that there is a constant \(c > 0\), independent of \(n\), such that

   \[ \max_{x\in[-1,1]} |\ell_X(x)| \geq c\, \frac{(2/e)^n}{n} \qquad \text{for all } n. \]

Note that \(2/e \approx 0.736\): parts 2 and 3 together show that \(\max\limits_{x\in[-1,1]} |\ell_X(x)|\) decays exponentially, but at a rate no better than \((2/e)^n\).

Exercise 6.6

(\(\star\) more difficult)

Fix \(\alpha > 0\) and let \(f(x) = \dfrac{1}{1+\alpha x^2}\), as in Example 6.2 (with \(\alpha=25\)).

1. Using the geometric series \(f(x) = \sum_{k=0}^\infty (-\alpha)^k x^{2k}\) for \(|x| < 1/\sqrt\alpha\), show that

   \[\begin{align} f^{(2k)}(0) &= (2k)!\,(-\alpha)^k, \nonumber\\ f^{(2k+1)}(0) &= 0. \nonumber \end{align}\]

2. Using partial fractions (the poles of \(f\) are at \(x = \pm i/\sqrt\alpha\)), show more generally that

   \[ \max_{x \in [-1,1]} |f^{(n)}(x)| = n!\, \alpha^{n/2}, \]

   with the maximum attained at \(x=0\).

3. Combine part 2 with Exercise 6.5 and Theorem 6.2 to show that, for \(X\) the equispaced nodes on \([-1,1]\),

   \[ |f(x) - I_Xf(x)| \leq C\sqrt\alpha \, \frac{1}{\sqrt n} \left( \frac{2\sqrt\alpha}{e} \right)^n, \]

   and deduce that this bound converges to \(0\) whenever \(\alpha < \frac{e^2}{4} \approx 1.847\).

It turns out that we have convergence for all \(\alpha < \alpha_{\rm c} \approx 3.621\) Trefethen (2019), Epperson (1987).

Next week, we will see that for Chebyshev nodes, the node polynomial instead satisfies \(|\ell_X(x)| \leq 2^{1-n}\) and so

\[ |f(x) - I_Xf(x)| \leq 2\sqrt\alpha \left( \frac{\sqrt\alpha}{2} \right)^n, \]

which converges to \(0\) whenever \(\alpha < 4\). (A sharper argument, using that \(f\) is analytic throughout a neighbourhood of \([-1,1]\) in \(\mathbb C\) rather than just bounding derivatives along the real axis shows that Chebyshev interpolation converges for all \(\alpha > 0\).)

Lecture 7

Exercise 7.1

Show that \(T_0(x) = 1\) and \(T_1(x) = x\).

Exercise 7.2

Use \(\cos(x+y) = \cos x \cos y - \sin x \sin y\) to show that

• \(T_{n+1}(x) = 2x T_{n}(x) - T_{n-1}(x)\) for \(n =1,2,\dots\),
• Hence show that \(T_n\) is a polynomial of degree \(n\),
• Show that the coefficent of \(x^n\) in \(T_n\) is \(2^{n-1}\) for \(n \geq 1\).

Exercise 7.3

Show that

• \(\|T_n\|_{L^\infty([-1,1])} = 1\),
• \(T_{n}\big( \cos\tfrac{j\pi}{n} \big) = (-1)^j\) for \(j = 0,\dots,n\) (extreme points; Chebyshev nodes of the second kind),
• \(T_{n+1} = 0\) on \(\big\{ \cos \frac{(2j+1)\pi}{2(n+1)} \big\}_{j=0}^{n}\) (zeros; Chebyshev nodes of the first kind).

Exercise 7.4

Recall that \(T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x)\) for \(n=1,2,...\).

1. Use the recursion to write down \(T_2\), \(T_3\) and \(T_4\) in the monomial basis.
2. Write down the Chebyshev nodes of the first kind \(X_{\mathrm I}\) for \(n = 2\) (i.e. the zeros of \(T_3\)) and check directly that \(\ell_{X_{\mathrm I}}(x) = t_3(x) = x^3 - \tfrac34 x\).
3. Find \(\| \ell_{X_{\mathrm I}} \|_{L^\infty([-1,1])}\) for \(n=2\) by calculus, and compare with the equispaced nodes \(\{-1, 0, 1\}\).

Exercise 7.5

In the notes, we claimed that \(\ell_{X_{\mathrm{II}}}(x) = 2^{-n} \big( T_{n+1}(x) - T_{n-1}(x) \big)\) for \(n \geq 1\). Since \(\|T_n\|_{L^\infty} = 1\), this gives \(\| \ell_{X_{\rm II}} \|_{L^\infty} \leq 2^{1-n}\).

Prove the claim as follows:

1. Show that \(T_{n+1}(\cos\theta) - T_{n-1}(\cos\theta) = -2 \sin\theta \sin n\theta\).
2. Deduce that \(T_{n+1} - T_{n-1}\) vanishes on \(X_{\mathrm{II}} = \big\{ \cos\frac{j\pi}{n} \big\}_{j=0}^n\).
3. Show that \(2^{-n} \big( T_{n+1} - T_{n-1} \big)\) is monic of degree \(n+1\), and conclude.

Exercise 7.6

Suppose we want to interpolate on a general interval \([a,b]\) rather than \([-1,1]\).

1. Let \(\phi(t) = \frac{a+b}{2} + \frac{b-a}{2} t\), so that \(\phi\) maps \([-1,1]\) onto \([a,b]\). If \(Y = \{y_0,\dots,y_n\} \subset [-1,1]\) and \(X = \phi(Y) = \{\phi(y_0), \dots, \phi(y_n)\}\), show that \[ \ell_X\big(\phi(t)\big) = \left( \frac{b-a}{2} \right)^{n+1} \ell_Y(t). \]
2. Deduce that, if \(X\) is the image of \(X_{\mathrm I}\) under \(\phi\), then \(\max_{x \in [a,b]} |\ell_X(x)| = 2 \left( \frac{b-a}{4} \right)^{n+1}\).
3. Let \(f(x) = \sin x\) on \([0, \pi]\), interpolated at the image of \(X_{\mathrm I}\). Write down an error bound and find the smallest \(n\) that guarantees \(\max_{x\in[0,\pi]} |f(x) - I_Xf(x)| \leq 10^{-6}\).

Quiz 2

Quiz 2: Exercise 1

Recall that \(g : [a,b] \to \mathbb R\) is a contraction if there exists \(L < 1\) such that \(|g(x) - g(y)| \leq L |x-y|\) for all \(x,y \in [a,b]\).

Prove the uniqueness part of the contraction mapping theorem. That is, suppose that \(g :[a,b] \to [a,b]\) is a contraction with a fixed point \(\xi = g(\xi) \in [a,b]\) and show that \(\xi\) is the only fixed point of \(g\) in \([a,b]\).

Suppose \(\zeta = g(\zeta) \in [a,b]\) is a fixed point. Then

\[\begin{align} |\xi - \zeta| = |g(\xi) - g(\zeta)| \leq L |\xi-\zeta| \nonumber \end{align}\]

and so \((1 - L)|\xi - \zeta| \leq 0\). Since \(L < 1\), we must have \(|\xi - \zeta| = 0\) and thus \(\xi = \zeta\).

Quiz 2: Exercise 2

Show that \(g(x) = \cos x\) maps \([0,\tfrac\pi3]\) to \([0,\tfrac\pi3]\) and is a contraction on \([0,\tfrac\pi3]\). What does this tell you about the iteration \(x_{n+1} = \cos x_n\) with \(x_1 = 0\)?

Since \(\cos\) is decreasing on \([0,\tfrac\pi3]\), we have \(0 < \tfrac12 = \cos \tfrac\pi3 \leq g(x) \leq \cos 0 = 1 < \tfrac\pi3\) for all \(x \in [0, \tfrac\pi3]\) and so \(g : [0,\tfrac\pi3] \to [\tfrac12, 1] \subset [0,\tfrac\pi3]\).

Moreover, \(|g'(x)| = |\sin x| \leq \sin \tfrac\pi3 = \frac{\sqrt{3}}{2} < 1\) for all \(x \in [0,\tfrac\pi3]\) and so, by the mean value theorem, \(g\) is a contraction with \(L = \frac{\sqrt3}{2}\).

By the contraction mapping theorem, there exists a unique fixed point \(\xi = \cos \xi \in [0,\tfrac\pi3]\) and, since \(x_1 = 0 \in [0,\tfrac\pi3]\), the iteration \(x_{n+1} = \cos x_n\) converges to \(\xi\) (at least linearly, with asymptotic error constant \(\mu = |g'(\xi)| = \sin \xi\)).

Quiz 2: Exercise 3

Recall that if \(f\) is \(n+1\) times continuously differentiable and \(X = \{x_0 < \dots < x_n\}\), then, for all \(x\), there exists \(\xi_x\) (between the smallest and largest of \(x, x_0, \dots, x_n\)) such that

\[\begin{align} f(x) - I_Xf(x) = \frac{f^{(n+1)}(\xi_x)}{(n+1)!} (x - x_0)(x-x_1) \cdots (x-x_n). \nonumber \end{align}\]

Let \(f(x) = x^3\) and \(X = \{ x_0 < x_1 < x_2 \} := \{-1, 0, 1\}\). Find a polynomial interpolant \(p\) of \(f\) on \(X\). Write down the interpolation error formula (above) in this case.

We have \(f(-1) = -1\), \(f(0) = 0\), \(f(1) = 1\). The polynomial \(p(x) = x \in \mathcal P_2\) satisfies \(p(x_j) = f(x_j)\) for \(j = 0,1,2\) and so, by uniqueness, \(I_Xf(x) = x\). (Alternatively, use the Lagrange polynomials: \(I_Xf(x) = -\ell_0(x) + \ell_2(x) = -\frac{x(x-1)}{2} + \frac{x(x+1)}{2} = x\).)

Here \(n = 2\) and \(f'''(x) = 6\) for all \(x\), so the error formula gives

\[\begin{align} f(x) - I_Xf(x) = \frac{6}{3!}(x+1)x(x-1) = x^3 - x, \nonumber \end{align}\]

which agrees with \(x^3 - x\) (for any choice of \(\xi_x\)).

Quiz 2: Exercise 4

Let \(f(x) = \frac1x\) and \(X = \{1, 2\}\).

1. Find the linear interpolant \(I_Xf\).
2. Use the interpolation error formula to show that \(|f(x) - I_Xf(x)| \leq \frac14\) for all \(x \in [1,2]\). Hint: \(|(x-1)(x-2)|\) is largest on \([1,2]\) at \(x = \frac32\).
3. Compute the actual error at \(x = \frac32\).

1. Using the Lagrange polynomials \(\ell_0(x) = \frac{x-2}{1-2} = 2 - x\) and \(\ell_1(x) = \frac{x-1}{2-1} = x - 1\), we have

\[\begin{align} I_Xf(x) = f(1) \ell_0(x) + f(2) \ell_1(x) = (2 - x) + \tfrac12 (x-1) = \frac{3 - x}{2}. \nonumber \end{align}\]

2. Here \(n = 1\) and \(f''(x) = \frac{2}{x^3}\), so \(\max_{\xi \in [1,2]} |f''(\xi)| = 2\). Moreover, \(|(x-1)(x-2)|\) is maximised on \([1,2]\) at the midpoint \(x = \frac32\) where it equals \(\frac14\). Therefore,

\[\begin{align} |f(x) - I_Xf(x)| \leq \frac{\max_{[1,2]}|f''|}{2!} |(x-1)(x-2)| \leq \frac{2}{2} \cdot \frac14 = \frac14. \nonumber \end{align}\]

3. \(f(\frac32) = \frac23\) and \(I_Xf(\frac32) = \frac34\), so the error is \(|\frac23 - \frac34| = \frac1{12} \leq \frac14\).

(In fact, \(f(x) - I_Xf(x) = \frac{(x-1)(x-2)}{2x}\), which agrees with the error formula with \(\xi_x = \sqrt[3]{2x}\), e.g. \(\xi_{3/2} = \sqrt[3]{3} \in [1,2]\).)

Lecture 8

Exercise 8.1

Let’s use the formula

\[ \begin{align} I_Xf(x) &= \sum_{j=0}^n f(x_j) \ell_j(x), \\ % \ell_j(x) &:= \prod_{k\not=j} \frac{x - x_k}{x_j - x_k} \nonumber \end{align} \]

to compute the polynomial interpolation of \(f(x) = \frac1{1+25x^2}\) (Runge function, as above):

Code

Runge = t -> 1/(1 + 25t^2)

function ℓⱼ( j, t, x ) 
    if j == 0
        return prod( @. t - x[2:end] ) / prod( @. x[1] - x[2:end] )
    elseif j == length(x)-1 
        return prod( @. t - x[1:end-1] ) / prod( @. x[end] - x[1:end-1] )
    else
        return prod( @. t - x[1:j] ) * prod( @. t - x[(j+2):end] ) / ( prod( @. x[j+1] - x[1:j] ) * prod( @. x[j+1] - x[(j+2):end] ) )
    end
end 

function p_naive( t, x, y ) 
    r = 0
    for j in 0:(length(x)-1)
        r = r + ℓⱼ(j, t, x) * y[j+1] 
    end
    return r
end 

egrid = range(-1, 1, length=401)

anim = @animate for n ∈ [10, 20, 50, 100, 200, 400, 600, 700, 760, 780:20:1020...]
    x = @. cos( π*(0:n)/n )
    y = @. Runge( x )
    p = p_naive.( egrid, Ref(x), Ref(y) )
    err = @. max( abs( p - Runge(egrid) ), 1e-17 )

    plt1 = plot(Runge, -1, 1; label=L"y = f(x)", lw = 3, linestyle = :dash, 
            title = "Naive Lagrange formula: n = $n", ylims = (-0.5, 1.5))
    plot!(plt1, egrid, p; label=L"y = p(x)", lw = 3 )
    scatter!(plt1, x, y; primary = false, markersize = 2, markerstrokewidth = 0)

    plt2 = plot(egrid, err; yscale = :log10, ylims = (1e-17, 1e2), legend = false, lw = 2,
            title = @sprintf("max error = %.2e", maximum(err)))
    plot( plt1, plt2; size=(1000, 500) )
end

gif(anim, "../images/Naive Lagrange.gif", fps = 1, show_msg = false)

What is happening here?

Underflow! The products \(\prod_{k} (t - x_k)\) and \(\prod_{k\not=j} (x_j - x_k)\) are of size roughly \(2^{-n}\) (and the partial products computed in the code can be even smaller). Around \(n \approx 800\), these partial products fall below the smallest normal double (\(\approx 2.2 \times 10^{-308}\)) into the subnormal range, where precision is lost.

Exercise 8.2

Suppose that \(X = \{ x_0 < ... < x_n \}\) and let \(\ell_j\) be the corresponding Lagrange polynomials. Show that \(\sum_{j=0}^n \ell_j(x) = 1\) for all \(x\).

Exercise 8.3

Explain why one may instead consider any weights of the form \(w_j := C \lambda_j\) in the Barycentric formula (8.3). Why might this be a good idea? See (8.4) for the weights appearing in the Barycentric formula for Chebyshev nodes of the second kind.

Exercise 8.4

Recall that the degree \(2\) polynomial interpolation of \(f(x) = \sin x\) on \(X = \{ x_0 < x_1 < x_2 \} = \{ 0, \frac\pi2, \pi \}\) is \(I_X f(x) = -\frac4{\pi^2} x(x-\pi)\). Write the Barycentric version of this formula (8.3).

Explain how you would have to update this formula if we added an extra point \(x_{3} = -2\pi\) to the interpolation set.

The Barycentric formula (8.3) depends on the weights \(\lambda_j = \prod_{k\not=j} (x_j-x_k)^{-1}\) which in this case are given by

\[ \begin{align} \lambda_0 &= \frac{1}{(0-\frac\pi2)(0-\pi)} = \frac2{\pi^2}, \nonumber\\ \lambda_1 &= \frac{1}{(\frac\pi2-0)(\frac\pi2 - \pi)} = -\frac4{\pi^2}, \nonumber\\ \lambda_2 &= \frac{1}{(\pi - 0)(\pi-\frac\pi2)} = \frac2{\pi^2} \end{align} \]

(which makes sense because \(\{ 0, \frac\pi2, \pi \}\) are Chebyshev nodes of the second kind scaled to \([0,\pi]\). The scaling only multiplies all of the weights by a common constant, which cancels in the Barycentric formula. In general, the Chebyshev weights can be chosen to be \(\frac12, -1, 1, -1, \dots, (-1)^{n-1}, \frac{(-1)^n}{2}\).)

Since \(f(x_{j}) = 0\) for \(j=0,2\), we have

\[ \begin{align} I_Xf(x) &= \left.\left( \frac{\lambda_0 f(x_0)}{x - x_0} + \frac{\lambda_1 f(x_1)}{x - x_1} + \frac{\lambda_2 f(x_2)}{x - x_2} \right) \middle/ \left( \frac{\lambda_0}{x - x_0} + \frac{\lambda_1}{x - x_1} + \frac{\lambda_2}{x - x_2} \right) \right. \nonumber\\ % &= \left.\left( \frac{-1}{x - \frac\pi2} \right) \middle/ \left( \frac{1/2}{x} + \frac{-1}{x - \frac\pi2} + \frac{1/2}{x - \pi} \right) \right. \end{align} \]

where, in the second line, we have multiplied all of the weights by \(\frac{\pi^2}{4}\) (which does not change the formula) to obtain the weights \(\frac12, -1, \frac12\).

We would have to update \(\lambda_{j}\) by dividing by \(x_{j} - (-2\pi)\) for each \(j=0,1,2\) so we get \(\lambda_0,\lambda_1,\lambda_2 = \frac{1}{\pi^3}, -\frac{8}{5\pi^3}, \frac{2}{3\pi^3}\) and write down

\[ \begin{align} \lambda_3 = \frac{1}{(-2\pi-0)(-2\pi-\frac\pi2)(-2\pi-\pi)} = -\frac{1}{15\pi^3} \end{align} \]

Exercise 8.5

Let \(n = 2\) and let \(X = X_{\mathrm{II}} = \big\{ \cos\frac{j\pi}{2} \big\}_{j=0}^2\) be the Chebyshev nodes of the second kind, so that \(x_0 = 1\), \(x_1 = 0\), and \(x_2 = -1\). Compute \(\lambda_0, \lambda_1, \lambda_2\) directly from the definition \(\lambda_j = \big( \prod_{k\not=j} (x_j - x_k) \big)^{-1}\) and check that they agree with (8.4).

We have \[ \begin{align} \lambda_0 &= \frac{1}{(1 - 0)(1 - (-1))} = \frac12, \nonumber\\ \lambda_1 &= \frac{1}{(0 - 1)(0 - (-1))} = -1, \nonumber\\ \lambda_2 &= \frac{1}{(-1 - 1)(-1 - 0)} = \frac12. \nonumber \end{align} \] With \(n = 2\), (8.4) gives \(\lambda_0 = \frac{2^2}{8} = \frac12\), \(\lambda_1 = (-1)^1 \frac{2^2}{4} = -1\), and \(\lambda_2 = (-1)^2 \frac{2^2}{8} = \frac12\).

Exercise 8.6

In the notes, we claimed that for the Chebyshev nodes \(X_{\mathrm{II}} = \big\{ x_j = \cos\frac{j\pi}{n} \big\}_{j=0}^n\) the barycentric weights \(\lambda_j := \big( \prod_{k\not=j} (x_j - x_k) \big)^{-1}\) are given by

\[ \lambda_j = \begin{cases} \frac{2^{n}}{4n} & \text{if } j = 0 \\ (-1)^j \frac{2^{n}}{2n} & \text{if } j = 1,\dots,n-1 \\ (-1)^n \frac{2^{n}}{4n} & \text{if } j = n \end{cases} \]

Prove the claim as follows:

1. Show that \(\lambda_j = 1 / \ell_X'(x_j)\), where \(\ell_X(x) = \prod_{k=0}^n (x - x_k)\) is the node polynomial.
2. Show that \(T_n'(\cos\theta) = n \frac{\sin n\theta}{\sin\theta}\). Using Exercise 7.5, deduce that \[ \ell_X(x) = \frac{2^{1-n}}{n} (x^2 - 1) T_n'(x). \]
3. By differentiating \(T_n(\cos\theta) = \cos n\theta\) twice, show that \(T_n\) satisfies the Chebyshev differential equation \[ (1 - x^2) T_n''(x) - x T_n'(x) + n^2 T_n(x) = 0. \]
4. Let \(w(x) := (x^2-1) T_n'(x)\). Show that \[ w'(x) = x T_n'(x) + n^2 T_n(x). \]
5. Compute \(w'(x_j)\) for \(j = 0,\dots,n\) and conclude.
6. Explain why we may replace \(\lambda_j\) by \((-1)^j\) (halved for \(j = 0, n\)) in the barycentric formula. Why is this a good idea when \(n\) is large?

For 5, treat the interior nodes and the endpoints separately: \(x_1,\dots,x_{n-1}\) are the points where \(T_n\) attains \(\pm 1\) in \((-1,1)\). For the endpoints, evaluate the differential equation from 3 at \(x = \pm 1\) to find \(T_n'(\pm1)\).

1. Differentiating \(\ell_X(x) = (x - x_j) \prod_{k\not=j} (x - x_k)\) with the product rule and setting \(x = x_j\) gives \(\ell_X'(x_j) = \prod_{k\not=j} (x_j - x_k) = 1/\lambda_j\).

2. Differentiating \(T_n(\cos\theta) = \cos n\theta\) gives \(-\sin\theta \, T_n'(\cos\theta) = -n \sin n\theta\). Therefore, with \(x = \cos\theta\), \[ (x^2 - 1) T_n'(x) = -\sin^2\theta \cdot n \frac{\sin n\theta}{\sin \theta} = -n \sin\theta \sin n\theta = \frac{n}{2} \big( T_{n+1}(x) - T_{n-1}(x) \big) \] by Exercise 7.5. Since both sides are polynomials which agree on \([-1,1]\), they are equal, and so \(\ell_X(x) = 2^{-n}\big( T_{n+1}(x) - T_{n-1}(x) \big) = \frac{2^{1-n}}{n} (x^2-1) T_n'(x)\).

3. Differentiating \(-\sin\theta \, T_n'(\cos\theta) = -n \sin n\theta\) again gives \(\sin^2\theta \, T_n''(\cos\theta) - \cos\theta \, T_n'(\cos\theta) = -n^2 \cos n\theta = -n^2 T_n(\cos\theta)\). With \(x = \cos\theta\) and \(\sin^2\theta = 1 - x^2\), this is the Chebyshev differential equation on \([-1,1]\) and therefore everywhere (both sides are polynomials).

4. By the product rule, \(w' = 2x T_n' + (x^2-1) T_n''\). By 3, \((x^2-1) T_n'' = n^2 T_n - x T_n'\) and so \(w' = x T_n' + n^2 T_n\).

5. For \(j = 1,\dots,n-1\), \(T_n(x_j) = \cos j\pi = (-1)^j\) and, since \(x_j\) is an interior maximum or minimum of \(T_n\), \(T_n'(x_j) = 0\). Therefore \(w'(x_j) = n^2 (-1)^j\).

At \(x = 1\), the differential equation gives \(-T_n'(1) + n^2 T_n(1) = 0\) and so \(T_n'(1) = n^2\). Therefore, \(w'(x_0) = w'(1) = n^2 + n^2 = 2n^2\). Similarly, at \(x = -1\), we have \(T_n'(-1) = -n^2 T_n(-1) = (-1)^{n+1} n^2\) and so \(w'(x_n) = w'(-1) = (-1)^n n^2 + (-1)^n n^2 = 2 (-1)^n n^2\).

By 1 and 2, \(\lambda_j = 1 / \ell_X'(x_j) = \frac{n 2^{n-1}}{w'(x_j)}\), which gives \[ \begin{align} \lambda_0 = \frac{n 2^{n-1}}{2n^2} = \frac{2^n}{4n}, \qquad \lambda_j = \frac{n 2^{n-1}}{(-1)^j n^2} = (-1)^j \frac{2^n}{2n}, \qquad \lambda_n = \frac{n 2^{n-1}}{2(-1)^n n^2} = (-1)^n \frac{2^n}{4n}. \nonumber \end{align} \]

6. The barycentric formula is unchanged if every \(\lambda_j\) is multiplied by the same constant \(c \not= 0\) (the factor \(c\) cancels between the numerator and the denominator). Taking \(c = \frac{2n}{2^n}\) gives the weights \((-1)^j\), halved for \(j = 0, n\). For large \(n\), \(2^n\) overflows in floating point arithmetic (for example, 2.0^1024 == Inf) whereas the rescaled weights are all of size \(1\) or \(\frac12\).

Exercise 8.7

Let \(X = \{x_0 < \dots < x_n\}\) be the equispaced nodes on \([-1,1]\): \(x_j = -1 + jh\) for \(j = 0, \dots, n\), where \(h = \frac2n\).

1. Show that \(\prod_{k\not= j} (x_j - x_k) = (-1)^{n-j} h^n \, j!\, (n-j)!\), and deduce that the barycentric weights are \[ \lambda_j = \frac{(-1)^{n-j}}{h^n\, n!} \binom{n}{j}. \]
2. Using Exercise 8.3, explain why we may use the weights \(w_j := (-1)^j \binom{n}{j}\) in the Barycentric formula (8.3) instead. Compare the ratio \(\max_j |w_j| / \min_j |w_j|\) for these weights with the same ratio for the weights \((-1)^j \tau_j\) for Chebyshev nodes (\(\tau_j = \frac12\) for \(j = 0,n\) and \(\tau_j = 1\) otherwise) as \(n \to \infty\).
3. Write a function EquispacedInterpolant(f, n), analogous to ChebyshevInterpolant from the notes, which returns the polynomial interpolant of f on the equispaced nodes (using the weights \(w_j\) from 2; in Julia, binomial(n, j) computes \(\binom{n}{j}\)). Compute the maximum error \(\max_{t} |f(t) - I_Xf(t)|\) over the grid range(-1, 1, length=2001) for:
   • the Runge function \(f(x) = \frac{1}{1 + 25x^2}\) with \(n = 10, 20, 40\),
   • \(f(x) = e^x\) with \(n = 10, 20, 40, 60\).
   Compare your results with the error estimate from Theorem 6.2 together with the bound \(\| \ell_X \|_{L^\infty([-1,1])} \leq \frac{e}{2\sqrt n} \big( \frac2e \big)^n\) for equispaced nodes from Lecture 6. Comment on what you see.

For 1, split the product into the terms with \(k < j\) (where \(x_j - x_k = (j - k)h > 0\)) and the terms with \(k > j\) (where \(x_j - x_k = -(k - j)h < 0\)). For 3, think about how errors of size \(10^{-16}\) in the values \(f(x_j)\) (or in the arithmetic) are amplified when the weights vary by many orders of magnitude.

1. For \(k < j\), we have \(x_j - x_k = (j-k)h\) and so \(\prod_{k=0}^{j-1} (x_j - x_k) = h^j \, j!\). For \(k > j\), we have \(x_j - x_k = -(k-j)h\) and so \(\prod_{k=j+1}^{n} (x_j - x_k) = (-1)^{n-j} h^{n-j} (n-j)!\). Multiplying gives \(\prod_{k\not= j} (x_j - x_k) = (-1)^{n-j} h^n \, j! \, (n-j)!\) and therefore \[ \lambda_j = \frac{(-1)^{n-j}}{h^n \, j! \, (n-j)!} = \frac{(-1)^{n-j}}{h^n \, n!} \binom{n}{j}. \]

2. We have \(w_j = C \lambda_j\) with the constant \(C = (-1)^n h^n n!\) (which does not depend on \(j\)), and so the factor \(C\) cancels between the numerator and the denominator of (8.3). The smallest weights are \(|w_0| = |w_n| = 1\) and the largest is \(\binom{n}{\lfloor n/2 \rfloor} \approx 2^n \sqrt{\frac{2}{\pi n}}\). Therefore, the ratio grows exponentially: it is \(252\) for \(n = 10\), about \(1.8 \times 10^5\) for \(n = 20\), and about \(1.2 \times 10^{17}\) for \(n = 60\). For Chebyshev nodes, the ratio is \(2\) for all \(n\).

3. For example:

function EquispacedInterpolant( f, n )
    x = range(-1, 1, length=n+1)
    y = f.(x)
    w = [ (-1.0)^j * binomial(n, j) for j in 0:n ]

    return function ( t )
        num = 0.0
        den = 0.0
        for j ∈ 1:(n+1)
            if t == x[j]
                return y[j]
            end
            num += w[j] * y[j] / (t - x[j])
            den += w[j] / (t - x[j])
        end
        return num / den
    end
end

This gives the following maximum errors:

\(n\)	Runge	\(e^x\)	error bound for \(e^x\)
\(10\)	\(1.9\)	\(2.4 \times 10^{-10}\)	\(1.4 \times 10^{-9}\)
\(20\)	\(6.0 \times 10^{1}\)	\(1.4 \times 10^{-12}\)	\(3.5 \times 10^{-23}\)
\(40\)	\(1.0 \times 10^{5}\)	\(3.9 \times 10^{-7}\)	\(8.2 \times 10^{-56}\)
\(60\)		\(2.0 \times 10^{-1}\)	\(9.5 \times 10^{-93}\)

Here, the error bound is \(\frac{e}{(n+1)!} \cdot \frac{e}{2\sqrt n} \big(\frac2e\big)^n\) (using \(\max_{[-1,1]} |f^{(n+1)}| = e\)). For the Runge function, the error grows exponentially: this is the Runge phenomenon from Lecture 6 (the derivatives of the Runge function grow too quickly for the error estimate to give convergence). For \(e^x\), the exact interpolant converges extremely quickly, and the computed error decreases until \(n \approx 20\). After that, the computed error grows even though the exact interpolant is more and more accurate. This is caused by rounding errors: perturbing the data \(f(x_j)\) by \(\delta_j\) changes the interpolant by \(\sum_j \delta_j \ell_j(x)\), and for equispaced nodes \(\max_x \sum_j |\ell_j(x)|\) grows exponentially in \(n\) (like the ratio of the weights in 2). Rounding errors of size \(10^{-16}\) are therefore amplified by a factor of order \(2^n\). For Chebyshev nodes, this amplification factor only grows like \(\log n\) (see Berrut and Trefethen (2004)).

Lecture 9

Exercise 9.1

Let \(f(x) = e^x\) and \(X = \{\{0, 0, 1\}\}\) (so that \(z_0 = 0\) with \(m_0 = 2\), \(z_1 = 1\) with \(m_1 = 1\), and \(n = 2\)). Find the Hermite interpolant \(p \in \mathcal P_2\) of \(f\) on \(X\) by writing \(p(x) = a_0 + a_1 x + a_2 x^2\) and solving for the coefficients. Compare \(p\) with the Taylor polynomial of degree \(1\) about \(0\) and with the Lagrange interpolant of \(f\) on \(\{0, 1\}\).

Exercise 9.2

What is the least degree polynomial interpolant on \(X = \{\{x_0, x_0, \dots, x_0\}\}\) (where \(x_0\) is repeated \(n+1\) times)?

Exercise 9.3

Return to Exercise 9.1: compute the Hermite interpolant of \(f(x) = e^x\) on \(\{\{0, 0, 1\}\}\) again, this time using the divided difference table and Theorem 9.3.

Exercise 9.4

What is the Hermite interpolant of \(f(x) = \sin x\) on \(\{\{ 0,0,\pi, \pi\}\}\)? Hint: you can compute the divided difference tables in the same way as before.

Exercise 9.5

In this exercise, we prove Theorem 9.1. As in the notes, let \(z_0, \dots, z_k\) be distinct nodes with multiplicities \(m_0, \dots, m_k \geq 1\), and let \(n + 1 = \sum_{j=0}^k m_j\).

1. Let \(r \in \mathcal P_n\) and suppose that \(r^{(l)}(z) = 0\) for \(l = 0, \dots, m-1\). Show that \(r(x) = (x - z)^m s(x)\) for some polynomial \(s\).
2. Suppose that \(r \in \mathcal P_n\) satisfies \(r^{(l)}(z_j) = 0\) for all \(l = 0,\dots,m_j - 1\) and \(j = 0,\dots,k\). Show that \(r = 0\).
3. Deduce that the Hermite interpolation problem has at most one solution in \(\mathcal P_n\).
4. Consider the map \(L : \mathcal P_n \to \mathbb R^{n+1}\) defined by \(L p := \big( p^{(l)}(z_j) \big)_{j = 0,\dots,k, \; l = 0, \dots, m_j - 1}\). Show that \(L\) is linear and injective, and conclude that the Hermite interpolation problem has exactly one solution in \(\mathcal P_n\).

For 1, write \(r\) as its Taylor polynomial about \(z\) (which is exact because \(r\) is a polynomial). For 2, apply 1 at \(z_0\) and then show that \(s\) still has zeros of order \(m_j\) at \(z_j\) for \(j \geq 1\). For 4, recall that an injective linear map between vector spaces of the same finite dimension is also surjective.

1. Since \(r \in \mathcal P_n\), Taylor’s theorem is exact: \(r(x) = \sum_{l=0}^n \frac{r^{(l)}(z)}{l!} (x - z)^l = \sum_{l=m}^n \frac{r^{(l)}(z)}{l!} (x - z)^l = (x-z)^m s(x)\) where \(s(x) := \sum_{l=m}^n \frac{r^{(l)}(z)}{l!} (x - z)^{l-m}\).

2. By 1, \(r(x) = (x - z_0)^{m_0} s_0(x)\). Fix \(j \geq 1\) and write \(u(x) := (x - z_0)^{m_0}\) (so that \(u(z_j) \not= 0\) because the nodes are distinct). By the Leibniz rule, \[ r^{(l)}(z_j) = \sum_{i=0}^l \binom{l}{i} u^{(l-i)}(z_j) \, s_0^{(i)}(z_j) = u(z_j) \, s_0^{(l)}(z_j) + \sum_{i=0}^{l-1} \binom{l}{i} u^{(l-i)}(z_j) \, s_0^{(i)}(z_j). \] By induction on \(l = 0, \dots, m_j - 1\), we obtain \(s_0^{(l)}(z_j) = 0\). Therefore, \(s_0\) satisfies the same conditions at \(z_1, \dots, z_k\) and, repeating the argument, we have \[ r(x) = (x - z_0)^{m_0} (x - z_1)^{m_1} \cdots (x - z_k)^{m_k} s(x) \] for some polynomial \(s\). If \(s \not= 0\), then \(r\) has degree at least \(\sum_j m_j = n+1\), which contradicts \(r \in \mathcal P_n\). Therefore \(s = 0\) and \(r = 0\).

3. If \(p, q \in \mathcal P_n\) both solve the Hermite interpolation problem, then \(r := p - q \in \mathcal P_n\) satisfies \(r^{(l)}(z_j) = f^{(l)}(z_j) - f^{(l)}(z_j) = 0\) for all \(j\) and \(l = 0, \dots, m_j-1\). By 2, \(r = 0\) and so \(p = q\).

4. Differentiation and evaluation are linear and so \(L\) is linear. If \(L p = 0\), then \(p = 0\) by 2, and so \(L\) is injective. Since \(\dim \mathcal P_n = n + 1 = \dim \mathbb R^{n+1}\), the rank-nullity theorem shows that \(L\) is also surjective. In particular, there exists \(p \in \mathcal P_n\) with \(L p = \big( f^{(l)}(z_j) \big)_{j,l}\), which is exactly a solution to the Hermite interpolation problem. By 3, this solution is unique.

Exercise 9.6

In this exercise, we prove the error estimate for Hermite interpolation (compare with the Lagrange case in Lecture 6). We say that \(g\) has a zero of order \(m\) at \(z\) if \(g(z) = g'(z) = \dots = g^{(m-1)}(z) = 0\).

1. (Rolle’s theorem with multiplicities) Suppose that \(g \in C^1([a,b])\) has zeros at distinct points \(w_1 < \dots < w_r\) in \([a,b]\) of orders \(k_1, \dots, k_r\), and let \(N := k_1 + \dots + k_r\). Show that \(g'\) has zeros in \([w_1, w_r]\) whose orders add up to at least \(N - 1\).
2. Let \(f \in C^{n+1}([a,b])\), let \(X = \{\{x_0, \dots, x_n\}\} \subset [a,b]\) be a multiset, and let \(p \in \mathcal P_n\) be the Hermite interpolant of \(f\) on \(X\). Define the node polynomial \(\ell_X(x) := (x - x_0)(x - x_1) \cdots (x - x_n)\) (so repeated nodes give repeated factors). Show that, for every \(x \in [a,b]\), there exists \(\xi_x \in [\min (\{x\} \cup X), \max (\{x\} \cup X)]\) with \[ f(x) - p(x) = \frac{f^{(n+1)}(\xi_x)}{(n+1)!} \ell_X(x). \]

Exercise 9.7

Let \(X = \{\{0, 0, 1, 1\}\}\) and suppose that we know \(f(0) = y_0\), \(f'(0) = s_0\), \(f(1) = y_1\) and \(f'(1) = s_1\).

1. Compute the divided difference table and show that the Hermite interpolant of \(f\) on \(X\) is \[ p(x) = y_0 + s_0 x + (d - s_0) x^2 + (s_0 + s_1 - 2d) \, x^2 (x - 1), \qquad \text{where } d := y_1 - y_0. \]
2. Write down \(p\) for \(f(x) = e^x\) and compare it with Exercise 9.1.
3. Use Exercise 9.6 to show that \(\max_{x \in [0,1]} |e^x - p(x)| \leq \frac{e}{384}\).
4. More generally, suppose that \(p\) is the Hermite interpolant of \(f\) on \(\{\{a, a, b, b\}\}\). Show that \(\max_{x \in [a,b]} |f(x) - p(x)| \leq \frac{(b-a)^4}{384} \max_{[a,b]} |f^{(4)}|\). What does this tell you about approximating \(f\) on \([a,b]\) by a cubic Hermite interpolant on each of \(N\) subintervals of length \(h = \frac{b-a}{N}\)?

Exercise 9.8

1. Let \(f \in C^n([a,b])\) and let \(X = \{\{x_0, \dots, x_n\}\} \subset [a,b]\) be a multiset. Show that there exists \(\xi \in [\min X, \max X]\) with \[ f[x_0, \dots, x_n] = \frac{f^{(n)}(\xi)}{n!}. \]
2. Deduce that \(f[x_0, x_0 + h] \to f'(x_0)\) and, more generally, \(f[x_0, x_0 + h, \dots, x_0 + nh] \to \frac{f^{(n)}(x_0)}{n!}\) as \(h \to 0\). Explain why this makes the first case of Theorem 9.2 a natural definition.

Exercise 9.9

1. Write a Julia function hermite_dd(x, df) which takes a vector of nodes x (sorted so that equal nodes are next to each other) and a function df(l, z) returning \(f^{(l)}(z)\), and returns the Newton coefficients \(f[x_0], f[x_0,x_1], \dots, f[x_0,\dots,x_n]\) using the divided difference table.
2. Write a function newton_eval(c, x, t) which evaluates \(\sum_{k=0}^n c_k \pi_k(t)\) using \(O(n)\) operations.
3. Check your code against Exercise 9.1, Exercise 9.4, and Exercise 9.7, and compute \(\max_{t \in [0,1]} |e^t - p(t)|\) for the cubic Hermite interpolant of \(e^x\) on \(\{\{0,0,1,1\}\}\).

Berrut, Jean-Paul, and Lloyd N. Trefethen. 2004. “Barycentric Lagrange Interpolation.” SIAM Review 46 (3): 501–17. https://doi.org/10.1137/S0036144502417715.

Driscoll, Tobin A., and Richard J. Braun. 2017. Fundamentals of Numerical Computation. Society for Industrial; Applied Mathematics. https://doi.org/10.1137/1.9781611975086.

Epperson, James F. 1987. “On the Runge Example.” The American Mathematical Monthly 94 (4): 329–41. https://doi.org/10.2307/2323093.

Trefethen, Lloyd N. 2019. Approximation Theory and Approximation Practice, Extended Edition. Philadelphia, PA: Society for Industrial; Applied Mathematics. https://doi.org/10.1137/1.9781611975949.