Interpolation

Numerical Analysis

Urbain Vaes

NYU Paris

Existence and uniqueness

Theorem (Existence and uniqueness)

Given data \(x_0 < x_1< \dotsc < x_n\) and \(y_0, \dotsc, y_n\), there exists a unique \(p \in \mathcal P_n\) such that \[ \forall i \in \{0, \dotsc, n\}, \qquad p(x_i) = y_i. \]

Proof of existence: Define Lagrange polynomials \[ \forall i \in \{0, \dotsc, n\}, \qquad L_i(x) := \prod_{\substack{j = 0 \\ j \neq i}}^{n} \frac {x - x_j} {x_i - x_j}. \] Then set \(p(x) = \sum_{i=0}^{n} y_i L_i(x)\) and verify that \(p\) is an interpolating polynomial, which proves existence.

Proof of uniqueness: assume \(q \in \mathcal P_n\) is another interpolating polynomial. Then \(p - q \in \mathcal P_n\) and \[ \forall i \in \{0, \dotsc, n\}, \qquad p(x) - q(x) = 0. \] Thus \(p - q\) is a polynomial of degree at most \(n\) with \(n + 1\) roots, so \(p - q = 0\).

Illustration: Lagrange polynomials

Key property of Lagrange polynomials: \[ L_i(x_j) = \begin{cases} 1 & \text{if $i = j$} \\ 0 & \text{if $i \neq j$} \end{cases} \]

using Plots
using LaTeXStrings

# Interpolation nodes
x = [1, 2, 3]

function L(i, z)
    result = 1.0
    for j in 1:length(x)
        if i != j
            result *= (z - x[j]) / (x[i] - x[j])
        end
    end
    return result
end

scatter(x, ones(3), title="Lagrange polynomials")
plot!(legend=:outertopright, size=(900, 900))
plot!(z -> L(1, z), label=L"L_1")
plot!(z -> L(2, z), label=L"L_2")
plot!(z -> L(3, z), label=L"L_3")

Interpolation error

Theorem

Assumptions

  • \(u\colon [a, b] \to \mathbf R\) is a \(C^{n+1}([a, b])\) function

  • \(a=x_0< x_1< \dotsc < x_n=b\) are \(n+1\) distinct interpolation nodes

  • \(\widehat u\) interpolates \(u\) at \(x_0, x_1, \dotsc, x_n\), i.e. \(\widehat u(x_i) = u(x_i)\) for all \(i ∈ \{0,\dotsc,n\}\)

Then, it holds that \(\quad ∀\, x ∈ [a,b],\quad ∃\, ξ=ξ(x) ∈ [a,b]\) \[ e_n(x) := u(x) - \widehat u(x) = \frac{u^{(n+1)}(\xi)}{(n+1)!} (x-x_0) \dotsc (x - x_n) \]

Corollary

Theorem

Assumptions

  • \(u\colon [a, b] \to \mathbf R\) is a smooth function

  • \(a=x_0< x_1< \dotsc < x_n=b\) are \(n+1\) distinct interpolation nodes

  • \(\widehat u\) interpolates \(u\) at \(x_0, x_1, \dotsc, x_n\).

Then it holds that \[ E_n := \sup_{x \in [a, b]} \bigl\lvert e_n(x) \bigr\rvert \leqslant\frac{C_{n+1}}{4(n+1)} h^{n+1}, \qquad C_{n+1} := \sup_{x \in [a, b]} \left\lvert u^{(n+1)}(x) \right\rvert. \] where \(h\) is the maximum spacing between two successive interpolation nodes.

Chebyshev roots

Question? What nodes \((x_0, \dotsc x_n)\) to pick to minimize \((x-x_0) \dotsc (x - x_n)\)? Recall \[ u(x) - \widehat u(x) = \frac{u^{(n+1)}(\xi)}{(n+1)!} (x-x_0) \dotsc (x - x_n) \]

Theorem

Assume that \(p\) is a monic polynomial of degree \(n \geqslant 1\): \[ p(x) = \alpha_0 + \alpha_1 x + \dotsb + \alpha_{n-1} x^{n-1} + x^n. \] Then it holds that \[ \label{eq:chebychev_lower_bound} \sup_{x \in [-1, 1]} \bigl\lvert p(x) \bigr\rvert \geqslant\frac{1}{2^{n-1}} =: E. \] In addition, the lower bound is achieved for \(p_*(x) = 2^{-(n-1)} T_n(x)\), where \(T_n\) is the Chebyshev polynomial of degree \(n\): \[ \label{eq:chebyshev_polynomial} T_n(x) = \cos(n\arccos x) \qquad (-1 \leqslant x \leqslant 1). \]

\(\leadsto\) Idea: use roots of \(T_n\) as interpolation nodes

Sine function with equidistant nodes

\[ u(x)=\sin{x} \] \[ x_k=a + (b-a) \frac{k}{n} \quad (0 ≤ k ≤ n) \]

Runge function with equidistant nodes

\[ u(x)=\frac{1}{1+25x^2} \] \[ x_k=a + (b-a) \frac{k}{n} \quad (0 ≤ k ≤ n) \]

Optimization of the nodes: Chebyshev nodes

Idea: Take \(x_0 < x_1 < \dotsc < x_n\) to be the roots of \(T_n\):

\[ x_k=-\cos \Bigl(\pi \frac{k + \frac{1}{2}}{n} \Bigl), \qquad (0 ≤ k < n) \]

Piecewise interpolation

Piecewise continuous interpolation

Let \(e_n\) denote the interpolation error and recall that \[\begin{align} &\max_{x ∈ [a, b]} \bigl\lvert e_n(x) \bigr\rvert \leqslant\frac{C_{n+1}}{4(n+1)} h^{n+1}, \qquad C_{n+1} = \max_{x ∈ [a, b]} \left\lvert u^{(n+1)}(x) \right\rvert \end{align}\] with \(h = \max_{i ∈ \{0, \dotsc, n-1\}} \lvert x_{i+1} - x_i\rvert\)

Idea to control the error :

  • Divide interval \([a,b]\) into \(n\) subinteravl of size \(h=\frac{b-a}{n}\)

  • Interpolate with a polynomial of degree \(m\) over \([x_{i},x_{i+1}]\)

    \[ \max_{x ∈ [x_{i},x_{i+1}]} \bigl\lvert e_n(x) \bigr\rvert \leqslant\frac{C_{m+1}}{4(m+1)} \left(\frac{h}{m}\right)^{m+1} \quad \Rightarrow \quad \max_{x ∈ [a, b]} \bigl\lvert e_n(x) \bigr\rvert \leqslant\frac{C_{m+1}}{4(m+1)} \left(\frac{h}{m}\right)^{m+1} \]

  • \(m\) is fixed and small but \(h\) can vary. The error scales as \(C h^{m+1}\).

Example (1/3)

Example (2/3)

Example (3/3)

Least squares approximation

Problem statement

  • \(n+1\) distinct nodes \(a=x_0< x_1< \dotsc < x_n=b\)

  • \(n+1\) values \(u_0, u_1, \dotsc, u_n\)

  • a vector space \({\rm Span}(φ_0,φ_1,\dotsc,φ_m)\) of continuous functions on \([a,b]\) with \(\color{red}{m<n}\)

We wish to find \(\widehat u(x) = α_0 φ_0(x) + \dotsb + α_m φ_m(x)\) such that \[ \forall i \in \{0, \dotsc, n\}, \qquad \widehat u(x_i) \, {\color{red}{\approx}} \, u_i \] Matrix form of these equations: \[ \mathsf A \mathbf{\boldsymbol{\alpha }}:= \begin{pmatrix} \varphi_0(x_0) & \varphi_1(x_0) & \dots & \varphi_m(x_0) \\ \varphi_0(x_1) & \varphi_1(x_1) & \dots & \varphi_m(x_1) \\ \varphi_0(x_2) & \varphi_1(x_2) & \dots & \varphi_m(x_2) \\ \vdots & \vdots & & \vdots \\ \varphi_0(x_{n-2}) & \varphi_1(x_{n-2}) & \dots & \varphi_m(x_{n-2}) \\ \varphi_0(x_{n-1}) & \varphi_1(x_{n-1}) & \dots & \varphi_m(x_{n-1}) \\ \varphi_0(x_n) & \varphi_1(x_n) & \dots & \varphi_m(x_n) \end{pmatrix} \begin{pmatrix} \alpha_0 \\ \alpha_1 \\ \vdots \\ \alpha_m \end{pmatrix} {\color{red}{\approx}} \begin{pmatrix} u_0 \\ u_1 \\ u_2 \\ \vdots \\ u_{n-2} \\ u_{n-1} \\ u_n \end{pmatrix} =: \mathbf{\boldsymbol{b}} \]

Solution via minimization problem

Key idea: Minimize \[ J(\mathbf{\boldsymbol{\alpha}})= \frac{1}{2}\,\sum_{i=0}^{n} {\lvert {\widehat u(x_i) - u_i} \rvert}^2 = \frac{1}{2}\,\sum_{i=0}^{n} \left( \sum_{j=0}^{m} \alpha_j \varphi_j(x_i) - u_i\right)^2 =\frac{1}{2}\,{\lVert {\mathsf A \mathbf{\boldsymbol{\alpha}}-\mathbf{\boldsymbol{b}}} \rVert}^2 \]

To this end, we seek \(\alpha\) such that \[ \nabla J(\mathbf{\boldsymbol{\alpha}})=\frac{1}{2}\,\nabla \Bigl( (\mathsf A \mathbf{\boldsymbol{\alpha }}- \mathbf{\boldsymbol{b}})^T(\mathsf A \mathbf{\boldsymbol{\alpha }}- \mathbf{\boldsymbol{b}}) \Bigr)=\mathbf{\boldsymbol{0}}, \] which leads to the equation \[ \nabla J(\mathbf{\boldsymbol{\alpha}})=\mathsf A^T(\mathsf A \mathbf{\boldsymbol{\alpha }}- \mathbf{\boldsymbol{b}}) =\mathbf{\boldsymbol{0}} \] Rearranging, we obtain the following linear system \[ \mathsf A^T\mathsf A \mathbf{\boldsymbol{\alpha }}= \mathsf A^T\mathbf{\boldsymbol{b}}. \]

Remark

In Julia, the result of the caluclation \(\alpha=(\mathsf A^T\mathsf A)^{-1} \mathsf A^T\mathbf{\boldsymbol{b}}\) can be obtained via

α = A\b

Polynomials.jl library

Example use of fit

using Plots, Polynomials
xs = range(0, 10, length=10)
ys = @. exp(-xs)
f1 = fit(xs, ys) # degree = length(xs) - 1
f2 = fit(xs, ys, 2) # degree = 2

scatter(xs, ys, label="Data")
plot!(f1, extrema(xs)..., label="Fit")
plot!(f2, extrema(xs)..., label="Quadratic Fit")
width, height = 750, 500 ; plot!(size = (width, height), legend = :topright)