In [1]:
%matplotlib inline
import numpy
import matplotlib.pyplot as plt
This gets complicated quickly!
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)$.
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)$.
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$.
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()
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()
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()
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:
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()
Defined by a recurrence relation
$$T_k(x) = 2 x T_{k-1}(x) - T_{k-2}(x)$$
Leading coefficient of $x^N$ in $T_N(x)$ is $2^{N-1}$ for $N \geq 1$
Extreme values:
$$|T_N(x)| \leq 1 ~~~~\text{for}~~~~~ -1 \leq x \leq 1$$
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}$$
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()
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.
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()
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()
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 - Use the data to define $d_k$ and preserve monotonicity if possible (avoid overshoots)
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()
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()