Interpolation¶

Interpolation is a fundamental numerical problem that is central to many other numerical algorithms. Simply put, given a finite number of points where values are known, find an interpolant, a (usually) continuous function that returns values everywhere and is guaranteed to pass through the known data.

Objectives:¶

• Define an interpolant

• Understand Polynomial Interpolation and Approximation

  • Define the interpolating polynomial $P_N(x)$ that exactly interpolates $N+1$ points

  • Calculation (in the monomial basis)

  • Uniqueness

  • Other bases (Lagrange, Newton?)

  • Error Analysis

  • Chebyshev Polynomials and optimal interpolation

• Understand Other interpolants

  • Piecewise Polynomial interpolation

  • Overlapping Polynomial interpolation

  • Cubic Splines and other smooth interpolants

  • Higher dimensional Interpolation schemes and scipy.interpolate

Interpolation (vs Fitting)¶

Definition: Given a discrete set of values $y_i$ at locations $x_i$, an interpolant is a (piece-wise) continuous function $f(x)$ that passes exactly through the data (i.e. $f(x_i) = y_i$).

A visual example for 3 random points

Comment: In general a polynomial of degree $N$ can be used to interpolate $N+1$ data points. There are many choices of functions to use to interpolate values, but here we focus on polynomials.

Applications¶

• Data filling

• Function approximation

• Fundamental component of other algorithms

  • Root finding (secant method)

  • Optimization, minima/maxima (successive parabolic interpolation)

  • Numerical integration and differentiation

  • The Finite Element Method

Polynomial Interpolation¶

Theorem: There is a unique polynomial of degree $N$, $P_N(x)$, that passes exactly through $N + 1$ values $y_0, y_1, \ldots, y_N$ at distinct points $x_0, x_1, \ldots, x_N$.

Consequence of the number of unknowns in $P_N(x)$.

Example 1: 2 Points¶

Given two points $(x_0, y_0)$ and $(x_1, y_1)$, There is a unique line

that connects them. We simply need to use the data to find $p_0$ and $p_1$:

We first note that we have two equations and two unknowns. The two equations can be found by assuming the function $P_1(x)$ interpolates the two data points

We can also (and should) write this problem as a small $2\times2$ linear algebra problem $A\mathbf{x}=\mathbf{b}$

Question: What are the unknowns, and where does the data sit in $A$ and $\mathbf{b}$?

With a bit of algebra you should be able to show that the solution of this problem is

and

which is just the equation of the straight line with slope $p_1$, that connects the two points.

Example 2: 3 Points¶

Given three points $(x_0, y_0)$, $(x_1, y_1)$, and $(x_2, y_2)$ leads to a quadratic polynomial:

Again, setting $P_2(x_i) = y_i$ yields 3 equations in 3 unknowns

Which reduces to the $3\times3$ linear system $A(\mathbf{x})\mathbf{p} = \mathbf{y}$

A more general approach to solving the system will be explored later, but first it is important to determine whether or not the system even has a solution.

Proof - Uniqueness of Polynomial Interpolants¶

Let

or

and require $\mathcal{P}_N(x_i) = y_i$ for $i=0,1,\ldots,N$ and $x_i \neq x_j ~~~ \forall i,j$.

Preliminaries: Monomial Basis¶

We can think of $\mathcal{P}_N(x) = \sum^N_{n=0} p_n x^n$ as a polynomial, or more fundamentally as a linear combination of a set of simpler functions, the monomials

with weights

respectively.

Linear independence of the Monomials¶

The monomials, form a linearly independent set of functions such that no monomial $x^n$ can be written as a linear combination of any other monomial. We can see this graphically, for the first few monomials

But more fundamentally. A set of functions is linearly independent if the only linear combination that add to form the zero function, e.g.

is if all the coefficients $p_i = 0$, $\forall i=0,\ldots N$

Theorem: The monomials $x^0,\ldots, x^n$ are linear independent.

Proof: consider $P_N(x) = 0$ for all $x$. Since the polynomials (and monomials) are differentiable at least $n$ times, differentiate $n$ times to yield

which implies $p_n=0$.

Using this result and differentiating $n-1$ times shows $p_{n-1}=0$, which by induction gives all $p_i = 0$.

Put another way, the only $n$th degree polynomial that is zero everywhere is if all coefficients are zero.

The Fundamental theorem of algebra¶

Every $n$th degree polynomial has exactly $n$ complex roots, i.e.

for $a_i\in \mathbb{C}$. Therefore, a non-trivial $n$th order polynomial can only be zero at at most $n$ points.

Proof - Uniqueness of Polynomial Interpolants¶

Let

interpolate the $N+1$ points $y_i$ at $x_i$.

i.e.

and $x_i \neq x_j ~~~ \forall i,j$.

Assume there exists another polynomial

that passes through the same set of points such that $Q_N(x_i) = y_i$. Now compute $T_N(x) = \mathcal{P}_N(x) - Q_N(x)$:

Now, by construction, $T_N(x_i) = 0$ which implies that it is equal to zero at $n+1$ points. However,

is a $n$th order polynomial which has at most $n$ real roots. The only way to reconcile this is if $T_n(x) = 0$, for all $x$, and therefore $p_n - q_n = 0$ individually and therefore $\mathcal{P}_N(x) = Q_N(x)$.

Example 3: Monomial Basis¶

Consider $\mathcal{P}_3(x) = p_0 + p_1 x + p_2 x^2 + p_3 x^3$ with the four data points $(x_i, y_i), ~~ i = 0,1,2,3$. We have four equations and four unknowns as expected:

Lets rewrite these as a matrix equation: First define the following vectors

When we write the system in matrix/vector form the matrix that arises is called a Vandermonde matrix:

We can now write the system of linear equations

as $V(\mathbf{x}) \mathbf{p} = \mathbf{y}$:

Note: the columns of $V$ are simply the monomial functions sampled at the discrete points $x_i$. i.e $V_j = \mathbf{x}^j$.

Because the monomials are linearly independent, so are the columns of $V$

• What happens if we have redundant data? Either $(x_i, y_i)$ is repeated or for one $i$ we have two values of $y$.

• What if we have more points then the order of polynomial we want?

• How does this relate to solving the above linear system of equations?

Vandermonde matrices in general are defined as

where $V$ is a $m \times n$ matrix with points $(x_i, y_i)$ for $i = 0, 1, 2, 3, \ldots m$ and for an order $N$ polynomial $\mathcal{P}_N(x)$.

Finding $p_i$¶

Finding the coefficients of $\mathcal{P}_N(x)$ can be done by solving the system outlined above. There are functions in numpy that can do this for us such as:

• numpy.polyfit(x, y, x.shape[0] - 1)

• numpy.vander(x, N=None) to construct the matrix and use a linear solver routine.

We can also use a different basis that might be easier to use.

Basis¶

Def: A basis for a $N$ dimensional vector space is a set of linearly independent vectors that span the space.

The monomials, $1,x,\ldots, x^n$, form the usual basis for the vector space of $n$th degree polynomials $P_N(x)$.

Example $P_2(x)$ is the space of all quadratic functions. i.e. $P_2(x) = \mathrm{span}< 1,x,x^2>$

i.e for every vector $\mathbf{p}\in\mathbb{R}^3$, there is a unique quadratic function in $P_2(x)$. (we say $P_2$ is isomorphic to $\mathbb{R}^3$ and is a three dimensional function space).

However, the monomials are not the only basis for $P_N$

Lagrange Basis¶

Given $N+1$ points $(x_0,y_0), (x_1,y_1), \ldots, (x_{N},y_{N})$ again assuming the $x_i$ are all unique, the interpolating polynomial $\mathcal{P}_N(x)$ can also be written as

where

are the Lagrange Polynomials

Lagrange Polynomials¶

A Key property of the Lagrange polynomials is that

which is why the weights in $P_N(x)$ are simply the $y$ values of the interpolant

Visualizing the Lagrange Polynomials¶

Solving for the coefficients of $P_N(x)$¶

In general, if

where $\phi_j(x)$ is any basis function for $P_N$ (i.e. monomial, Lagrange, and there are many more). Then finding the unique set of weights for the interpolating polynomial through $N+1$ distinct data points $(x_i, y_i)$, just reduces to solving $N+1$ linear equations $P_N(x_i) = y_i$.

For the monomial basis this reduces to the linear system

What is the matrix $V$, for the Lagrange Basis? What are the weights $\mathbf{w}$?

A little pen-and-paper work here can be instructive$\ldots$

Linear Independence of the Lagrange Polynomials¶

Because the weights of each basis function in the Lagrange basis is just the $y$ value at the interpolation points, it is straightforward to show that the Lagrange polynomials are linearly independent. I.e. the statement

is equivalent to interpolating the zero function, where all the $w_j =0$

Example : $N = 1$ 1st order general Lagrange Polynomial as linear interpolant¶

Given 2 points $(x_0, y_0)$ and $(x_1, y_1)$ the Lagrange form of $\mathcal{P}_N(x)$ is given by

and

and the 1st order interpolating Polynomial is simply

The behavior of $\mathcal{P}_1(x)$ becomes clearer if we note that on the interval $x\in[x_0,x_1]$ that

is simply the fractional distance across the interval.

We should also note that

(show this)

Then the interpolating polynomial is simply $$

$$

for $s\in[0,1]$ which is just the linear line segment that connects points $y_0$ and $y_1$

As a specific example we will plot $\ell_0(x)$, $\ell_1(x)$ and $\mathcal{P}_1$ for the interval $x=[1,3]$ and $y_0=2$, $y_1=3$

Example: Interpolate six points from $\sin(2\pi x)$¶

Use six evenly spaced points to approximate $\sin$ on the interval $x \in [-1, 1]$. What is the behavior as $N \rightarrow \infty$? Also plot the absolute error between $f(x)$ and the interpolant $P_N(x)$.

Example 6: Runge’s Function¶

Interpolate $f(x) = \frac{1}{1 + 25 x^2}$ using 6 points of your choosing on $x \in [-1, 1]$.

Try it with 11 points.

Keep increasing the number of points and see what happens.

Example 7: Weierstrass “Monster” Function¶

Defined as

such that

This function is continuous everywhere but not differentiable anywhere.

Rules of Thumb¶

• Avoid high-order interpolants when possible! Keep increasing the number of points and see what happens.

• Avoid extrapolation - Increase the range of $x$ in the above example and check how good the approximation is beyond our sampling interval

Error Analysis¶

Theorem: Lagrange Remainder Theorem - Let $f(x) \in C^{N+1}[-1, 1]$, then

where $\mathcal{P}_N(x)$ is the interpolating polynomial and

where

is a monic polynomial of order $N+1$ with exactly $N+1$ roots at the nodes $x_i$

A few things to note:

• For Taylor’s theorem note that $Q(x) = (x - x_0)^{N+1}$ and the error only vanishes at $x_0$.

• For Lagrange’s theorem the error vanishes at all $x_i$.

• To minimize $R_N(x)$ requires minimizing $|Q(x)|$ for $x \in [-1, 1]$.

Minimizing $R_N(x)$¶

Minimizing the error $R_N(x)$ in Lagrange’s theorem is equivalent to minimizing $|Q(x)|$ for $x \in [-1, 1]$.

Minimizing error $\Leftrightarrow$ picking roots of $Q(x)$ or picking the points where the interpolant data is located. How do we this?

Chebyshev Polynomials¶

Chebyshev polynomials $T_N(x)$ are another basis that can be used for interpolation.

First 5 polynomials

In general, the Chebyshev polynomials are generated by a recurrence relation

Chebyshev nodes¶

Chebyshev polynomials have many special properties and locations including the location of their roots and extrema known as Chebyshev nodes

• Chebyshev nodes of the 1st kind (roots)

• Chebyshev nodes of the 2nd kind (extrema)

Properties of Chebyshev Polynomials¶

1. Defined by a recurrence relation

   $$T_0(x) = 1,\quad T_1(x) = x$$

   (61)
   $$T_k(x) = 2 x T_{k-1}(x) - T_{k-2}(x),\quad k\geq2$$

   (62)

2. Leading coefficient of $x^N$ in $T_N(x)$ is $2^{N-1}$ for $N \geq 1$

3. Extreme values:

   $$|T_N(x)| \leq 1 \quad \text{for} \quad -1 \leq x \leq 1$$

   (63)

Properties of Chebyshev Polynomials¶

4. Minimax principle: The polynomial

   $$T(x) = \frac{T_{N+1}(x)}{2^N}$$

   (64)

   is a monic polynomial, a univariate function with the leading coefficient equal to 1, with the property that

   $$\max |T(x)| \leq \max |Q(X)| \quad \text{for} \quad x \in [-1, 1], \quad \text{and}$$

   (65)
   $$\max |T(x)| = \frac{1}{2^N}$$

   (66)

Recall that the remainder term in the Lagrange Remainder Theorem was

with

Error Analysis Redux¶

Given that the Chebyshev polynomials are a minimum on the interval $[-1, 1]$ we would like $Q(x) = T(x)$.

Since we only know the roots of $Q(x)$ (the points where the interpolant data is located) we require these points to be the roots of the Chebyshev polynomial $T_{N+1}(x)$ therefore enforcing $T(x) = Q(x)$.

The zeros of $T_N(x)$ in the interval $[-1, 1]$ can be shown to satisfy

These nodal points (sampling the function at these points) can be shown to minimize interpolation error.

Summary¶

1. Minimizing the error in Lagrange’s theorem is equivalent to minimizing

   $$|Q(x)| \quad \text{for} \quad x \in [-1, 1].$$

   (70)

2. We know Chebyshev polynomials are a minimum on the interval $[-1, 1]$ so we would like to have $Q(x)=T(x)$.

3. Since we only know the roots of $Q(x)$ (the points where the interpolant data is located) we require these points to be the roots of the Chebyshev polynomial $T_{N+1}(x)$ therefore enforcing $T(x) = Q(x)$.

4. The zeros of $T_N(x)$ in the interval $[-1, 1]$ can be shown to satisfy

   $$x_k = \cos\left( \frac{(2k - 1) \pi}{2 N} \right ) \quad \text{for} \quad k=1, \ldots, N$$

   (71)

   These nodal points (sampling the function at these points) can be shown to minimize interpolation error.

Notes¶

• The Chebyshev nodes minimize interpolation error for any polynomial basis (due to uniqueness of the interpolating polynomial, any polynomial that interpolates these points are identical regardless of the basis).

• Chebyshev nodes uniquely define the Chebyshev polynomials.

• The boundedness properties of Chebyshev polynomials are what lead us to the roots as a minimization but there are other uses for these orthogonal polynomials.

• There are two kinds of Chebyshev nodes and therefore two definitions.

Example: Chebyshev Interpolation of Runge’s function¶

Comparison between interpolation at Chebyshev Nodes vs equally spaced points

Affine Transformation of intervals¶

While the chebyshev polynomials (and their roots) are defined on the interval

it turns out we can use them to approximate functions on any arbitrary interval $x\in[a,b]$ as long as we can map the interval $[-1,1]$ onto $[a,b]$. This is easy to do given the affine transformation

defined by a scaling $\beta$ and a shift $\alpha$

To solve for $\alpha$ and $\beta$, we simply transform the two endpoints such that $$

$$

or

or solving

Affine Transformation of intervals¶

Thus if we know the Chebyshev nodes for $s\in[-1, 1]$

then the correct interpolation points on $[a,b]$ are simply $x_k = x(s_k)$

For example: if $s_k=0$ (in the center of the interval), then

Note¶

Affine Transformation of the interval is useful for may problems and applications.

Piece-Wise Polynomial Interpolation¶

Given $N$ points, use lower order polynomial interpolation to fit the function in pieces. We can choose the order of the polynomials and the continuity.

• $C^0$: Interpolant is continuous

  • Linear interpolation

  • Quadratic interpolation

• $C^1$: Interpolation and 1st derivative are continuous

  • Cubic Hermite polynomials (PCHiP)

• $C^2$: Interpolation, 1st and 2nd derivatives are continuous

  • Cubic splines

Piece-Wise Linear¶

Given a segment between point $(x_k, y_k)$and $(x_{k+1}, y_{k+1})$ we can do linear interpolation for all $x\in[x_k, x_{k+1}]$ as

where

is the fractional distance of $x$ within the segment (i.e. $s=\ell_1(x)$ for $x\in[x_k, x_{k+1}]$)

Piece-Wise Quadratic $C^0$ Polynomials¶

Use every three points $(x_{k+1}, y_{k+1})$, $(x_{k}, y_{k})$, and $(x_{k-1}, y_{k-1})$, to find quadratic interpolant and define final interpolant $P(x)$ using the quadratic interpolant $\mathcal{P}_k(x)$ on $[x_{k-1}, x_{k+1}]$.

Transformation of intervals¶

The previous algorithms are quite direct but can look a bit messy because every interval (or set of points) appears different...An important trick that is used in finite-element analysis is to transform each interval (or element)

to the unit interval

by the affine transformation

and do all the interpolation in the transformed frame.

Example: Linear Interpolation¶

for Linear interpolation over an arbitrary interval $[x_k, x_{k+1}]$ we can use a Lagrange interpolant

where

and do all the interpolation in the transformed frame.

Substituting

into the definitions of $\ell_0(x)$ and $\ell_1(x)$ show that within this element

and

Piece-Wise $C^1$ Cubic Interpolation¶

For the previous two cases we had discontinous 1st derivatives! We can make this better by constraining the polynomials to be continuous at the boundaries of the piece-wise intervals.

Given a segment between points $(x_k, y_k)$ and $(x_{k+1}, y_{k+1})$ we want to fit a cubic function between the two points.

Now we have 4 unknowns but only two data points! Constraining the derivative at each interval end will lead to two new equations and therefore we can solve for the interpolant.

where we need to prescribe the $d_k$’s. Since we know the polynomial we can write these 4 equations as

Rewriting this as a system we get