In [1]:
%matplotlib inline
import numpy
import matplotlib.pyplot as plt

Interpolation

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

Applications

  • Data filling
  • Function approximation
  • Fundamental component of other algorithms
    • Root finding (secant method)
    • Optimization, minima/maxima (successive parabolic interpolation)
    • Numerical integration and differentiation

Polynomial Interpolation (1-D)

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

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

Example 1: 2 Points

Given points are $(x_0, y_0)$ and $(x_1, y_1)$ which will lead to a line:

Define $P_1(x) = p_1 x + p_0$ and use the two points to find $p_0$ and $p_1$:

$$y_0 = p_1 x_0 + p_0 ~~~~ \Rightarrow ~~~~ p_0 = y_0 - p_1 x_0$$
$$y_1 = p_1 x_1 + p_0 ~~~~ \Rightarrow ~~~~ y_1 = p_1 x_1 + y_0 - p_1 x_0 ~~~~ \Rightarrow ~~~~ p_1 = \frac{y_1 - y_0}{x_1 - x_0} ~~~~ \Rightarrow ~~~~ p_0 = y_0 - \frac{y_1 - y_0}{x_1 - x_0} x_0$$
$$P_1(x) = \frac{y_1 - y_0}{x_1 - x_0} x + y_0 - \frac{y_1 - y_0}{x_1 - x_0} x_0 = \frac{y_1 - y_0}{x_1 - x_0} (x - x_0) + y_0$$

Example 2: 3 Points

Given points are $(x_0, y_0)$, $(x_1, y_1)$, and $(x_2, y_2)$ which will lead to quadratic polynomial:

Define $P_2(x) = p_0 x^2 + p_1 x + p_2$ leading to the equations

$$y_0 = p_2 x_0^2 + p_1 x_0 + p_0$$$$y_1 = p_2 x_1^2 + p_1 x_1 + p_0$$$$y_2 = p_2 x_2^2 + p_1 x_2 + p_0$$

This gets complicated quickly!

Proof

Let

$$P_N(x) = \sum^N_{n=0} p_n x^n $$

or $$P_N(x) = p_0 + p_1 x + \cdots + p_{N - 1} x^{N - 1} + p_{N} x^N$$

and require $P_N(x_i) = y_i$ for $i=0,1,\ldots,N$ and $x_i \neq x_j ~~~ \forall i,j$.

Assume there exists another polynomial

$$Q_N(x) = \sum^N_{n=0} q_n x^n$$

that passes through the same set of points such that $Q_N(x_i) = y_i$. Now compute $T_N(x) = P_N(x) - Q_N(x)$:

We know that by assumption that $T_N(x_i) = 0$ but what about for all $x$?

$$T_N(x) = P_N(x) - Q_N(x) = \sum^N_{n=0} p_n x^n - q_n x^n = \sum^N_{n=0} (p_n - q_n) x^n$$

But if $T_N(x_i) = 0$ implies that $p_n - q_n = 0$ individually and therefore $P_N(x) = Q_N(x)$.

Monomial Basis

Let $P_N(x) = \sum^N_{n=0} p_n x^n$, $P_N(x)$ can be represented by a linear combination of the monomials

$$1, x, x^2, x^3, \ldots, x^{N-1}, x^N$$

with weights

$$p_0, p_1, p_2, p_3, \ldots, p_{N-1},~~ \text{and}~~ p_N$$

respectively.

Example 3: Monomial Basis

Consider $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:

$$P_3(x_0) = p_0 + p_1 x_0 + p_2 x_0^2 + p_3 x_0^3 = y_0$$$$P_3(x_1) = p_0 + p_1 x_1 + p_2 x_1^2 + p_3 x_1^3 = y_1$$$$P_3(x_2) = p_0 + p_1 x_2 + p_2 x_2^2 + p_3 x_2^3 = y_2$$$$P_3(x_3) = p_0 + p_1 x_3 + p_2 x_3^2 + p_3 x_3^3 = y_3$$

Lets rewrite these as a matrix equation:

$$\vec{x} = \begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \end{bmatrix} ~~~~ \vec{y} = \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \end{bmatrix} ~~~~ \vec{p} = \begin{bmatrix} p_0 \\ p_1 \\ p_2 \\ p_3 \end{bmatrix}$$

Define the Vandermonde matrix as

$$ V = \begin{bmatrix} 1 & x_0 & x_0^2 & x_0^3 \\ 1 & x_1 & x_1^2 & x_1^3 \\ 1 & x_2 & x_2^2 & x_2^3 \\ 1 & x_3 & x_3^2 & x_3^3 \end{bmatrix} $$

which allows us to write the system of linear equations as $V \vec{p} = \vec{y}$:

$$\begin{bmatrix} 1 & x_0 & x_0^2 & x_0^3 \\ 1 & x_1 & x_1^2 & x_1^3 \\ 1 & x_2 & x_2^2 & x_2^3 \\ 1 & x_3 & x_3^2 & x_3^3 \end{bmatrix} \begin{bmatrix} p_0 \\ p_1 \\ p_2 \\ p_3 \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \end{bmatrix}$$

Vandermonde matrices in general are defined as

$$V = \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^N \\ 1 & x_1 & x_1^2 & \cdots & x_1^N \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_m & x_m^2 & \cdots & x_m^N \\ \end{bmatrix} $$

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 $P_N(x)$.

Finding $p_i$

Finding the coefficients of $P_N(x)$ can be done by solving the system outlined above but there are other approaches:

  • Use numpy.polyfit(x, y, x.shape[0] - 1)
  • Use numpy.vander(x, N=None) to construct the matrix and use a linear solver routine.
  • Use a Lagrangian basis

Lagrangian 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 $P_N(x)$ can be written as

$$P_N(x) = \sum^{N}_{i=0} y_i \ell_i(x)$$

where

$$\ell_i(x) = \prod^{N}_{j=0, j \neq i} \frac{x - x_j}{x_i - x_j} = \frac{x - x_0}{x_i - x_0} \frac{x - x_1}{x_i - x_1} \cdots \frac{x - x_{i-1}}{x_i - x_{i-1}}\frac{x - x_{i+1}}{x_i - x_{i+1}} \cdots \frac{x - x_{N}}{x_i - x_{N}}$$

Note that $\ell_i(x_i) = 1$ and $\forall j\neq i, ~~ \ell_i(x_j) = 0$.

Example 4: $N = 1$ Lagrange Polynomial

Given 2 points $(x_0, y_0)$ and $(x_1, y_1)$ the Lagrange form of $P_N(x)$ is given by

$$\ell_0(x) = \frac{x - x_1}{x_0 - x_1}$$

and

$$\ell_1(x) = \frac{x - x_0}{x_1 - x_0}$$

so that

$$P_1(x) = \ell_0(x) \cdot y_0 + \ell_1(x) \cdot y_1 = \frac{x - x_1}{x_0 - x_1} \cdot y_0 + \frac{x - x_0}{x_1 - x_0} \cdot y_1$$

One important aspect of Lagrange polynomials to note is that the $\ell_i(x)$ functions are exactly 1 when $x = x_i$ and that every other $\ell_j(x)$ where $j \neq i$ is 0.


In [2]:
data = numpy.array([[-1.5, -0.5], [0.0, 0.0]])
N = data.shape[0] - 1
M = data.shape[0]
x = numpy.linspace(-2.0, 2.0, 100)

# ====================================================
# Compute the Lagrange basis (\ell_i(x))
lagrange_basis = numpy.ones((N + 1, x.shape[0]))
for i in xrange(N + 1):
    for j in xrange(N + 1):
        if i != j:
            print i, data[i, 0], data[j, 0]
            lagrange_basis[i, :] *= (x - data[j, 0]) / (data[i, 0] - data[j, 0])

# Calculate full polynomial
P = numpy.zeros(x.shape[0])
for n in xrange(N + 1):
    P += lagrange_basis[n, :] * data[n, 1]

# ====================================================
# Plot individual basis functions
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)
for i in xrange(N + 1):
    axes.plot(x, lagrange_basis[i, :], label="$\ell_{%s}(x)$" % i)
#     axes.plot(data[i, 0], data[i, 1], 'ko')

axes.set_title("Lagrange Basis $\ell_i(x)$")
axes.set_xlabel("x")
axes.set_ylabel("$\ell_i(x)$")
axes.legend(loc=8)

# Plot full polynomial P_N(x)
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)
axes.plot(x, P, label="$P_{%s}(x)$" % N)
for point in data:
    axes.plot(point[0], point[1], 'ko')
axes.set_title("$P_N(x)$")
axes.set_xlabel("x")
axes.set_ylabel("$P_N(x)$")
    
plt.show()


0 -1.5 0.0
1 0.0 -1.5

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

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


In [3]:
data = numpy.empty((6,2))
data[:, 0] = numpy.linspace(-1, 1, 6)
data[:, 1] = numpy.sin(2.0 * numpy.pi * data[:, 0])
N = data.shape[0] - 1
M = data.shape[0]
x = numpy.linspace(-1.0, 1.0, 100)

# ====================================================
# Compute the Lagrange basis (\ell_i(x))
lagrange_basis = numpy.ones((N + 1, x.shape[0]))
for i in xrange(N + 1):
    for j in xrange(N + 1):
        if i != j:
            lagrange_basis[i, :] *= (x - data[j, 0]) / (data[i, 0] - data[j, 0])

# Calculate full polynomial
P = numpy.zeros(x.shape[0])
for n in xrange(N + 1):
    P += lagrange_basis[n, :] * data[n, 1]

# ====================================================
# Plot individual basis functions
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)
for i in xrange(N + 1):
    axes.plot(x, lagrange_basis[i, :], label="$\ell_{%s}(x)$" % i)

axes.set_title("Lagrange Basis $\ell_i(x)$")
axes.set_xlabel("x")
axes.set_ylabel("$\ell_i(x)$")
# axes.legend(loc=8)

# Plot full polynomial P_N(x)
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)
axes.plot(x, P, label="$P_{%s}(x)$" % N)
axes.plot(x, numpy.sin(2.0 * numpy.pi * x), 'r--', label="True $f(x)$")
for point in data:
    axes.plot(point[0], point[1], 'ko')
axes.set_title("$P_N(x)$")
axes.set_xlabel("x")
axes.set_ylabel("$P_N(x)$")
    
plt.show()


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.


In [4]:
def f(x):
    return 1.0 / (1.0 + 25.0 * x**2)

x = numpy.linspace(-1, 1, 100)

data = numpy.empty((6, 2))
data[:, 0] = numpy.linspace(-1, 1, 6)
data[:, 1] = f(data[:, 0])
N = data.shape[0] - 1

# Calculate interpolant
lagrange_basis = numpy.ones((N + 1, x.shape[0]))
for i in xrange(N + 1):
    for j in xrange(N + 1):
        if i != j:
            lagrange_basis[i, :] *= (x - data[j, 0]) / (data[i, 0] - data[j, 0])

# Calculate full polynomial
P = numpy.zeros(x.shape[0])
for n in xrange(N + 1):
    P += lagrange_basis[n, :] * data[n, 1]
    
# Plot the results
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.plot(x, P, 'b', label="$P_6(x)$")
axes.plot(x, f(x), 'k', label="True $f(x)$")
axes.plot(data[:, 0], data[:, 1], 'ro', label="data")
          
axes.set_title("Interpolation of Runge's function")
axes.set_xlabel("x")
axes.set_ylabel("y")
axes.legend(loc=1)

plt.show()


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

$$f(x) = P_N(x) + R_N(x)$$

where $P_N(x)$ is the interpolating polynomial and

$$R_N(x) = Q(x) \frac{f^{(N+1)}(c)}{(N+1)!} ~~~~ \text{with}~~~~ c \in [-1,1]$$

and

$$Q(x) = (x-x_0)(x-x_1)\cdots(x-x_N).$$

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

Chebyshev Polynomials

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

First 5 polynomials $$T_0(x) = 1$$

$$T_1(x) = x$$$$T_2(x) = 2 x^2 - 1$$$$T_3(x) = 4 x^3 - 3 x$$$$T_4(x) = 8x^4 - 8x^2 + 1$$

In [5]:
x = numpy.linspace(-1, 1, 100)

fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.plot(x, numpy.ones(x.shape), label="$T_0$")
axes.plot(x, x, label="$T_1$")
axes.plot(x, 2.0 * x**2 - 1.0, label="$T_2$")
axes.plot(x, 4.0 * x**3 - 3.0 * x, label="$T_3$")
axes.plot(x, 8.0 * x**4 - 8.0 * x**2 + 1, label="$T_4$")

axes.set_ylim((-1.1, 1.1))
axes.set_title("Chebyshev Polynomials")
axes.set_xlabel("x")
axes.set_ylabel("$T_N(x)$")
axes.legend(loc=1)

plt.show()


Properties:

  1. Defined by a recurrence relation

    $$T_k(x) = 2 x T_{k-1}(x) - T_{k-2}(x)$$

  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 ~~~~\text{for}~~~~~ -1 \leq x \leq 1$$

  4. Minimax principle: The polynomial

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

    is a monic polynomial with the property that

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

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

Error Analysis Redux

Recall that in order to minimize the error with Lagrange's Theorem we needed to minimize $|Q(x)|$ for $x \in [-1, 1]$. This implies that in order to minimize $R_N(x)$ with Chebyshev polynomials set $Q(x) = T(x)$. We only know $Q(x)$ through its roots $x_0, x_1, \ldots, x_N$ so we will require that the points used (nodes) in the interpolation by the zeros of $T_{N+1}(x)$.

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 ) ~~~~~ \text{for}~~~~ k=0, 1, \ldots, N-1$$

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


In [6]:
x = numpy.linspace(0, numpy.pi, 100)
N = 15

fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.plot(numpy.cos(x), numpy.sin(x), 'r--')
axes.plot(numpy.linspace(-1.1, 1.1, 100), numpy.zeros(x.shape), 'r')
for k in xrange(1, N + 1):
    location = [numpy.cos((2.0 * k - 1.0) * numpy.pi / (2.0 * N)),
                numpy.sin((2.0 * k - 1.0) * numpy.pi / (2.0 * N))]
    axes.plot(location[0], location[1], 'ko')
    axes.plot(location[0], 0.0, 'ko')
    axes.plot([location[0], location[0]], [0.0, location[1]], 'k--')
    
axes.set_xlim((-1.1, 1.1))
axes.set_ylim((-0.1, 1.1))
plt.show()



In [7]:
def f(x):
    return 1.0 / (1.0 + 25.0 * x**2)

# Parameters
x = numpy.linspace(-1, 1, 100)
num_points = 6

# ============================================================
#  Equidistant nodes
equidistant_data = numpy.empty((num_points, 2))
equidistant_data[:, 0] = numpy.linspace(-1, 1, num_points)
equidistant_data[:, 1] = f(equidistant_data[:, 0])
N = equidistant_data.shape[0] - 1

lagrange_basis = numpy.ones((N + 1, x.shape[0]))
for i in xrange(N + 1):
    for j in xrange(N + 1):
        if i != j:
            lagrange_basis[i, :] *= (x - equidistant_data[j, 0]) / (equidistant_data[i, 0] - equidistant_data[j, 0])

# Calculate full polynomial
P_lagrange = numpy.zeros(x.shape[0])
for n in xrange(N + 1):
    P_lagrange += lagrange_basis[n, :] * equidistant_data[n, 1]

# ============================================================
#  Chebyshev nodes
chebyshev_data = numpy.empty((num_points, 2))
chebyshev_data[:, 0] = numpy.cos((2.0 * numpy.arange(N + 1) - 1.0) * numpy.pi / (2.0 * N))
chebyshev_data[:, 1] = f(chebyshev_data[:, 0])

# Project data onto Chebyshev polynomials
coeff = numpy.polynomial.chebyshev.chebfit(chebyshev_data[:, 0], chebyshev_data[:, 1], N)

# Find needed Chebyshev polynomials
T = numpy.empty((N, x.shape[0]))
T[0, :] = numpy.ones(x.shape)
T[1, :] = x
for k in xrange(2, N):
    T[k, :] = 2.0 * x * T[k-1, :]- T[k-2, :]
    
# Could plot them here but actually just want to evaluate the polynomial at x
P_cheby = numpy.zeros(x.shape)
for k in xrange(N):
    P_cheby += T[k, :] * coeff[k]

# ============================================================
# Plot the results
fig = plt.figure()
fig.set_figwidth(fig.get_figwidth() * 2.0)

axes = fig.add_subplot(1, 2, 1)
axes.plot(x, P_lagrange, 'b', label="$P_%s(x)$" % N)
axes.plot(x, f(x), 'k', label="True $f(x)$")
axes.plot(equidistant_data[:, 0], equidistant_data[:, 1], 'ro', label="data")
axes.set_title("Interpolation of Runge's function at Equispaced Points")
axes.set_xlabel("x")
axes.set_ylabel("y")
axes.legend(loc=1)

axes = fig.add_subplot(1, 2, 2)
axes.plot(x, P_cheby, 'b', label="$P_%s(x)$" % N)
axes.plot(x, f(x), 'k', label="True $f(x)$")
axes.plot(chebyshev_data[:, 0], chebyshev_data[:, 1], 'ro', label="data")
          
axes.set_title("Interpolation of Runge's function at Chebyshev Points")
axes.set_xlabel("x")
axes.set_ylabel("y")
axes.legend(loc=1)


plt.show()


/Users/mandli/src/envs/claw/lib/python2.7/site-packages/numpy/polynomial/chebyshev.py:1761: RankWarning: The fit may be poorly conditioned
  warnings.warn(msg, pu.RankWarning)

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)$a nd $(x_{k+1}, y_{k+1})$ define the segment as

$$\mathcal{P}_k(x) = \frac{y_{k+1} - y_k}{x_{k+1} - x_k} (x - x_k) + y_k$$

The final interpolant $P(x)$ is then defined on $[x_k, x_{k+1}]$ using this function.


In [8]:
data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 0.0], [6.0, 0.5], [9.0, -2.0], [9.5, -3.0]])
x = numpy.linspace(0.0, 10, 100)

# Lagrange Basis
N = data.shape[0] - 1

lagrange_basis = numpy.ones((N + 1, x.shape[0]))
for i in xrange(N + 1):
    for j in xrange(N + 1):
        if i != j:
            lagrange_basis[i, :] *= (x - data[j, 0]) / (data[i, 0] - data[j, 0])

# Calculate full polynomial
P_lagrange = numpy.zeros(x.shape[0])
for n in xrange(N + 1):
    P_lagrange += lagrange_basis[n, :] * data[n, 1]

# C^0 Piece-wise linear
# P_pw_linear = numpy.interp(x, data[:, 0], data[:, 1])
P_linear = numpy.zeros(x.shape)
for n in xrange(1, N + 1):
    P_linear += ((data[n, 1] - data[n - 1, 1]) / (data[n, 0] - data[n - 1, 0]) * (x - data[n - 1, 0]) 
                         + data[n - 1, 1]) * (x > data[n - 1, 0]) * (x <= data[n, 0])
# Add end points for continuity
P_linear += numpy.ones(x.shape) * data[0, 1] * (x < data[0, 0])
P_linear += numpy.ones(x.shape) * data[-1, 1] * (x >= data[-1, 0])

# Plot
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.plot(data[:,0], data[:,1], 'ko')
axes.plot(x, P_lagrange, 'b--')
axes.plot(x, P_linear, 'r')

axes.set_title("Interpolated Data - $C^0$ Linear")
axes.set_xlabel("x")
axes.set_ylabel("$P_1(x)$")
axes.set_xlim([0.0, 10.0])
axes.set_ylim([-3.0, 15.0])

plt.show()


Piece-Wise Overlapping Polynomials

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


In [9]:
data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 0.0], [6.0, 0.5], [9.0, -2.0], [9.5, -3.0]])
x = numpy.linspace(0.0, 10, 100)

# C^0 Piece-wise quadratic
P_quadratic = numpy.zeros(x.shape)
for k in xrange(1, N + 1, 2):
    p = numpy.polyfit(data[k - 1:k + 2, 0], data[k - 1:k + 2, 1], 2)
    P_quadratic += numpy.polyval(p, x) * (x > data[k - 1, 0]) * (x <= data[k + 1, 0])
# Add end points for continuity
P_quadratic += numpy.ones(x.shape) * data[0, 1] * (x < data[0, 0])
P_quadratic += numpy.ones(x.shape) * data[-1, 1] * (x >= data[-1, 0])
    
# Plot
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.plot(data[:,0], data[:,1], 'ko')
axes.plot(x, P_quadratic, 'r')

axes.set_title("Interpolated Data - $C^0$ Quadratic")
axes.set_xlabel("x")
axes.set_ylabel("$P_3(x)$")
axes.set_xlim([0.0, 10.0])
axes.set_ylim([-3.0, 15.0])

plt.show()


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.

$$\mathcal{P}_k(x) = p_0 + p_1 x + p_2 x^2 + p_3 x^3$$$$\mathcal{P}_k(x_k) = y_k, ~~~ \mathcal{P}_k(x_{k+1}) = y_{k+1}$$

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.

$$\frac{\text{d}}{\text{dx}} \mathcal{P}_k(x_k) = d_k, ~~~ \frac{\text{d}}{\text{dx}} \mathcal{P}_k(x_{k+1}) = d_{k+1}$$

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

$$p_0 + p_1 x_k + p_2 x_k^2 + p_3 x_k^3 = y_k$$$$p_0 + p_1 x_{k+1} + p_2 x_{k+1}^2 + p_3 x_{k+1}^3 = y_{k+1}$$$$p_1 + p_2 x_k + p_3 x_k^2 = d_k$$$$p_1 + p_2 x_{k+1} + p_3 x_{k+1}^2 = d_{k+1}$$

Can also rewrite as a system similar to before

$$\begin{bmatrix} 1 & x_k & x_k^2 & x_k^3 \\ 1 & x_{k+1} & x_{k+1}^2 & x_{k+1}^3 \\ 0 & 1 & x_k & x_k^2 \\ 0 & 1 & x_{k+1} & x_{k+1}^2 \end{bmatrix} \begin{bmatrix} p_0 \\ p_1 \\ p_2 \\ p_3 \end{bmatrix} = \begin{bmatrix} y_k \\ y_{k+1} \\ d_k \\ d_{k+1} \end{bmatrix}$$

We can also generalize this by letting $s \in [0, 1]$ and think of the problem with $(0, y_k)$ and $(1, y_{k+1})$. This simplifies the above system to

$$\begin{bmatrix} 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & 1 \end{bmatrix} \begin{bmatrix} p_0 \\ p_1 \\ p_2 \\ p_3 \end{bmatrix} = \begin{bmatrix} y_k \\ y_{k+1} \\ d_k \\ d_{k+1} \end{bmatrix}$$

which can be solved to find

$$\mathcal{P}(s) = (1-s)^2 (1 + 2s) y_k + s^2 (3 - 2 s) y_{k+1} + s (1 - s)^2 d_k - s^2 (1 - s)d_{k+1}$$$$\mathcal{P}'(s) = 6s(s-1) y_k + 6s(1-s) y_{k+1} + (s-1)(3s-1) d_k - s(3s-2) d_{k+1}$$$$\mathcal{P}''(s) = 6 (1-2s)(y_{k+1} - y_k) + (6s - 4) d_k + (6s-2) d_{k+1}$$

Now, how to choose $d_k$?

PCHIP

PCHIP - Use the data to define $d_k$ and preserve monotonicity if possible (avoid overshoots)

Rules
  1. Definte $\delta_k$ as slope of $k$th segment
  2. If $\text{sign}(\delta_k) \neq \text{sign}(\delta_{k-1})$ then $d_k = 0$
  3. Else if $\delta_k == 0$ or $\delta_{k-1} == 0$ then $d_k = 0$
  4. Else use weighted harmonic mean:

    Let $w_1 = 2 h_k + h_{k+1}$, $w_2 = h_k + 2 h_{k+1}$, and $h_k = x_k - x_{k-1}$ then

    $$\frac{1}{d_k} = \frac{1}{w_1 + w_2} \left(\frac{w_1}{\delta_{k-1}} + \frac{w_2}{\delta_k} \right )$$


In [10]:
import scipy.interpolate as interpolate

data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 0.0], [6.0, 0.5], [9.0, -2.0], [9.5, -3.0]])
x = numpy.linspace(0.0, 10, 100)

# C^1 Piece-wise PCHIP
P_pchip = interpolate.pchip_interpolate(data[:, 0], data[:, 1], x)

# Plot
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.plot(data[:,0], data[:,1], 'ro')
axes.plot(x, P_pchip, 'r')

axes.set_title("Interpolated Data - $C^1$ Cubic PCHIP")
axes.set_xlabel("x")
axes.set_ylabel("$P_3(x)$")
axes.set_xlim([0.0, 10.0])
axes.set_ylim([-3.0, 15.0])

plt.show()


Cubic Splines

Now enfore continuity on second derivatives as well:

$$\mathcal{P}''_{k-1}(x_k) = \mathcal{P}''_{k-1}(x_k)$$

From our generalization before we know

$$\mathcal{P}''(s) = 6 (1-2s)(y_{k+1} - y_k) + (6s - 4) d_k + (6s-2) d_{k+1}$$

and our constraint now becomes

$$\mathcal{P}''_{k-1}(1) = \mathcal{P}''_{k-1}(0)$$$$\mathcal{P}''_{k-1}(1) = 6 (1-2 \cdot 1)(y_{k} - y_{k-1}) + (6\cdot 1 - 4) d_{k-1} + (6\cdot 1-2) d_{k}$$$$\mathcal{P}''_{k-1}(0) = 6 (1-2 \cdot 0)(y_{k+1} - y_k) + (6\cdot 0 - 4) d_k + (6\cdot 0-2) d_{k+1}$$$$-6(y_{k} - y_{k-1}) + 2 d_{k-1} + 4 d_{k} = 6 (y_{k+1} - y_k) - 4 d_k -2 d_{k+1}$$

We now have constraints on choosing the $d_k$ values. Note that we still need to prescribe them at the boundaries of the full interval.


In [11]:
import scipy.interpolate as interpolate

data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 0.0], [6.0, 0.5], [9.0, -2.0], [9.5, -3.0]])
x = numpy.linspace(0.0, 10, 100)

# C^2 Piece-wise Splines
P_spline = interpolate.UnivariateSpline(data[:, 0], data[:, 1], s=0)

# Plot
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.plot(data[:,0], data[:,1], 'ro')
axes.plot(x, P_spline(x), 'r')

axes.set_title("Interpolated Data - $C^2$ Cubic Splines")
axes.set_xlabel("x")
axes.set_ylabel("$P_3(x)$")
axes.set_xlim([0.0, 10.0])
axes.set_ylim([-3.0, 15.0])

plt.show()



In [12]:
import scipy.interpolate as interpolate

data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 0.0], [6.0, 0.5], [9.0, -2.0], [9.5, -3.0]])
x = numpy.linspace(0.0, 10, 100)

# Lagrange Basis
N = data.shape[0] - 1

lagrange_basis = numpy.ones((N + 1, x.shape[0]))
for i in xrange(N + 1):
    for j in xrange(N + 1):
        if i != j:
            lagrange_basis[i, :] *= (x - data[j, 0]) / (data[i, 0] - data[j, 0])

# Calculate full polynomial
P_lagrange = numpy.zeros(x.shape[0])
for n in xrange(N + 1):
    P_lagrange += lagrange_basis[n, :] * data[n, 1]

# C^0 Piece-wise linear
# P_pw_linear = numpy.interp(x, data[:, 0], data[:, 1])
P_linear = numpy.zeros(x.shape)
for n in xrange(1, N + 1):
    P_linear += ((data[n, 1] - data[n - 1, 1]) / (data[n, 0] - data[n - 1, 0]) * (x - data[n - 1, 0]) 
                         + data[n - 1, 1]) * (x > data[n - 1, 0]) * (x <= data[n, 0])
    
# Add end points for continuity
P_linear += numpy.ones(x.shape) * data[0, 1] * (x < data[0, 0])
P_linear += numpy.ones(x.shape) * data[-1, 1] * (x >= data[-1, 0])
# C^0 Piece-wise quadratic
P_quadratic = numpy.zeros(x.shape)
for k in xrange(1, N + 1, 2):
    p = numpy.polyfit(data[k - 1:k + 2, 0], data[k - 1:k + 2, 1], 2)
    P_quadratic += numpy.polyval(p, x) * (x > data[k - 1, 0]) * (x <= data[k + 1, 0])
# Add end points for continuity
P_quadratic += numpy.ones(x.shape) * data[0, 1] * (x < data[0, 0])
P_quadratic += numpy.ones(x.shape) * data[-1, 1] * (x >= data[-1, 0])
    
# C^1 Piece-wise PCHIP
P_pchip = interpolate.pchip_interpolate(data[:, 0], data[:, 1], x)

# C^2 Piece-wise Splines
P_spline = interpolate.UnivariateSpline(data[:, 0], data[:, 1], s=0)
    
# Plot
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.plot(data[:,0], data[:,1], 'ko', label="Data")
axes.plot(x, P_lagrange, 'y', label="Lagrange")
axes.plot(x, P_linear, 'g', label="PW Linear") 
axes.plot(x, P_quadratic, 'r', label="PW Quadratic")
axes.plot(x, P_pchip, 'c', label="PW Cubic - PCHIP")
axes.plot(x, P_spline(x), 'b', label="PW Cubic - Spline")

axes.set_title("Interpolated Data - Method Comparisons")
axes.set_xlabel("x")
axes.set_ylabel("$P(x)$")
axes.legend()
axes.set_xlim([0.0, 10.0])
axes.set_ylim([-3.0, 15.0])

plt.show()