Spline functions - Numerical Analysis

#%config InlineBackend.figure_format = 'svg'
from pylab import *
from scipy.interpolate import barycentric_interpolate

Consider a grid for the interval $[a,b]$

a = x_0 < x_1 < \ldots < x_N = b

(1)

We have seem piecewise polynomial approximations which lead to globally continuous approximations. If we are given $\{f(x_i), f'(x_i)\}$ we can perform piecewise Hermite interpolation leading to globally $\cts^1$ approximations. But this requires knowledge of $f'$ . Alternately, spline approximations can lead to globally smooth approximations without knowing derivative information.

The derivative of a spline of order $m$ is a spline of order $m-1$ , and similarly for anti- derivatives. Note that

$m=2$ corresponds to piecewise linear function which is continuous, and
$m=3$ is a piecewise quadratic function which is continuous and has continuous derivatives Kincaid & Cheney, 2002.

6.10.1Cubic spline interpolation¶

Consider the case of $m=4$ where $s(x)$ is a cubic polynomial in every sub-interval $[x_{i-1},x_i]$ and we are given the values of some function

y_i = f(x_i), \qquad i=0,1,\ldots,N

(2)

Let us count the number of unknowns and number of equations available. We can write, for $i=1,2,\ldots,N$

s(x) = a_i + b_i x + c_i x^2 + d_i x^3, \qquad x \in [x_{i-1},x_i]

(3)

Hence we have $4N$ unknown coefficients $\{a_i,b_i,c_i,d_i\}$ . The spline function must satisfy some constraints; firstly, it must interpolate the given function values

s(x_i) = y_i, \qquad i=0,1,\ldots,N

(4)

and its derivatives must be continuous at the interior points

s^{(r)}(x_i^-) = s^{(r)}(x_i^+), \qquad i=1,2,\ldots,N-1, \qquad r=0,1,2

(5)

This gives us $(N+1)+3(N-1) = 4N-2$ equations, and we are short of two equations since we have $4N$ unknown coefficients. There are different ways in which the two missing equations can be supplied.

6.10.2Construction of cubic spline¶

Let

M_i = s''(x_i), \qquad i=0,1,\ldots,N

(6)

Since $s(x)$ is cubic on $[x_i,x_{i+1}]$ , then $s''(x)$ is linear and thus

s''(x) = \frac{(x_{i+1}-x)M_i + (x-x_i)M_{i+1}}{h_i}, \qquad i=0,1,\ldots,N-1

(7)

where

h_i = x_{i+1}-x_i

(8)

By the above construction $s''(x)$ is continuous on $[a,b]$ . Integrating twice, we get

\begin{align*} s(x) =& \frac{1}{6h_i}[ (x_{i+1}-x)^3 M_i + (x-x_i)^3 M_{i+1}] \\ & + C_i (x_{i+1}-x) + D_i (x-x_i), \quad\quad x \in [x_i,x_{i+1}] \end{align*}

(9)

with $C_i,D_i,M_i$ to be determined. Using the interpolation conditions $s(x_i)=y_i$ and $s(x_{i+1}) = y_{i+1}$ gives

C_i = \frac{y_i}{h_i} - \frac{h_iM_i}{6}, \qquad D_i = \frac{y_{i+1}}{h_i} - \frac{h_i M_{i+1}}{6}, \qquad i=0,1,\ldots,N-1

(10)

At this stage, the only unknowns are the $M_i$ values. Note that by construction, $s(x)$ and $s''(x)$ are continuous on $[a,b]$ . We will now enforce continuity of $s'(x)$ at the interior nodes. On $[x_i,x_{i+1}]$

s'(x) = \frac{1}{2 h_i}[-(x_{i+1}-x)^2 M_i + (x-x_i)^2 M_{i+1}] + \frac{y_{i+1}-y_i} {h_i} - \frac{(M_{i+1}-M_i)h_i}{6}

(11)

and on $[x_{i-1},x_i]$

s'(x) = \frac{1}{2 h_{i-1}}[-(x_{i}-x)^2 M_{i-1} + (x-x_{i-1})^2 M_{i}] + \frac{y_{i}- y_{i-1}}{h_{i-1}} - \frac{(M_{i}-M_{i-1})h_{i-1}}{6}

(12)

so that continuity of derivative at $x=x_i$ , $s'(x_i^-) = s'(x_i^+)$ , yields

\frac{h_{i-1}}{6}M_{i-1} + \frac{h_{i-1}+h_i}{3}M_i + \frac{h_i}{6}M_{i+1} = \frac{y_{i+1}-y_i}{h_i} - \frac{y_i - y_{i-1}}{h_{i-1}}, \quad 1 \le i \le N-1

(13)

We have $N-1$ equations for the $N+1$ unknowns $M_0, M_1, \ldots,M_N$ .

6.10.3End-point derivative conditions: complete cubic spline¶

To obtain the two extra equations, we impose

s'(x_0) = y_0', \qquad s'(x_N) = y_N'

(14)

where $y_0'=f'(x_0)$ , $y_N'=f'(x_N)$ are assumed to be known. This yields the equations

\begin{aligned} \frac{h_0}{3} M_0 + \frac{h_0}{6} M_1 &= \frac{y_1-y_0}{h_0} - y_0' \\ \frac{h_{N-1}}{6} M_{N-1} + \frac{h_{N-1}}{3}M_N &= y_N' - \frac{y_N - y_{N-1}} {h_{N-1}} \end{aligned}

(15)

We now have $N+1$ equations for the $N+1$ unknowns $\{M_0, M_1, \ldots, M_N\}$ . They can be written in matrix form

AM = b

(16)

where

b = \begin{bmatrix} \frac{y_1-y_0}{h_0} - y_0' \\ \left\{ \frac{y_{i+1}-y_i}{h_i} - \frac{y_i - y_{i-1}}{h_{i-1}} \right\}_{i=1}^{N-1} \\ y_N' - \frac{y_N - y_{N-1}}{h_{N-1}} \end{bmatrix}, \qquad M = \begin{bmatrix} M_0 \\ M_1 \\ \vdots \\ M_N \end{bmatrix}

(17)

A = \left[ \begin{aligned} & \frac{h_0}{3} \quad \frac{h_0}{6} \quad 0 \quad \cdots \\ & \textrm{diag}\left\{ \frac{h_{i-1}}{6}, \quad \frac{h_{i-1}+h_i}{3} \quad \frac{h_i} {6} \right\}_{i=1}^{N-1} \\ & \qquad\quad\quad \cdots \quad 0 \quad \frac{h_{N-1}}{6} \quad \frac{h_{N-1}}{3} \end{aligned} \right]

(18)

The matrix $A$ is symmetric, positive definite and diagonally dominant. Hence the problem is uniquely solvable. The matrix is tridiagonal and can be solved rapidly using about $8N$ arithmetic operations (Thomas Algorithm). The resulting cubic spline is called the complete cubic spline interpolant, and we denote it by $s_c(x)$ .

Theorem 1 (Error of complete cubic spline)

Let $f \in \cts^4[a,b]$ and given a partition

\tau_N: \quad a=x_0 < x_1 < \ldots < x_N = b

(19)

define

|\tau_N| := \max_{1 \le i \le N} (x_i - x_{i-1})

(20)

Let $s_c(x)$ be the complete cubic spline interpolant of $f(x)$ on the partition $\tau_N$

\begin{aligned} s_c(x_i) &= f(x_i), \qquad i=0,1,\ldots,N \\ s_c'(a) &= f'(a) \\ s_c'(b) &= f'(b) \end{aligned}

(21)

Then for some constants $c_j$

\max_{a \le x \le b}|f^{(j)}(x) - s_c^{(j)}(x)| \le c_j |\tau_N|^{4-j} \max_{a \le x \le b} |f^{(4)}(x)|, \qquad j=0,1,2

(22)

With the additional assumption

\frac{ |\tau_N|}{\min_{1 \le i \le N}(x_i - x_{i-1})} \le \gamma < \infty, \qquad N \to \infty

(23)

the above bound holds for $j=3$ also. Acceptable constants are

c_0 = \frac{5}{384}, \qquad c_1 = \frac{1}{24}, \qquad c_2 = \frac{3}{8}

(24)

For proof, see Boor, 1978.

Example 1

We will approximate the function

f(x) = (1-x^2)^2 \sin(4\pi x) \exp(\sin(2\pi x)), \qquad x \in [-1,+1]

(26)

for which

f'(-1) = f'(+1) = 0

(27)

which will be used in constructing the cubic spline below. The function looks like this.

Source

xmin, xmax = -1.0, +1.0
f = lambda x: (1-x**2)**2*sin(4*pi*x)*exp(sin(2*pi*x))

xe = linspace(xmin,xmax,200)
ye = f(xe)
plot(xe,ye,label='f(x)')
legend(), xlabel('x');

Let us first try polynomial interpolation using Chebyshev points.

N = 20
x = cos(linspace(0,pi,N+1))
y = f(x)
yp = barycentric_interpolate(x,y,xe)
plot(xe,ye,x,y,'o',xe,yp,'r-')
title("Chebyshev interpolation with "+str(N+1)+" points")
legend(('Exact','Data','Polynomial'));

We see that polynomial interpolation is rather oscillatory. Of course this will become better as we increase the number of points/degree. Using uniformly spaced points for polynomial interpolation would have given us even bad approximation. We will next try cubic spline interpolation. First we implement a function which constructs the matrix and right hand side.

def spline_matrix(x,y):
    N = len(x) - 1
    h = x[1:] - x[0:-1]
    A = zeros((N+1,N+1))
    A[0,0] = h[0]/3; A[0,1] = h[0]/6
    for i in range(1,N):
        A[i,i-1] = h[i-1]/6; A[i,i] = (h[i-1]+h[i])/3; A[i,i+1] = h[i]/6
    A[N,N-1] = h[N-1]/6; A[N,N] = h[N-1]/3

    r = zeros(N+1)
    r[0] = (y[1] - y[0])/h[0]
    for i in range(1,N):
        r[i] = (y[i+1]-y[i])/h[i] - (y[i]-y[i-1])/h[i-1]
    r[N] = -(y[N]-y[N-1])/h[N-1]
    return A,r

We will use uniformly spaced points for spline interpolation. We first construct the matrix and right hand side and solve the matrix problem. Finally, we evaluate the spline in each interval to make a plot.

x = linspace(xmin,xmax,N+1)
y = f(x)

A,r = spline_matrix(x,y)
m = solve(A,r)

plot(xe,ye,x,y,'o')
for i in range(0,N):
    h  = x[i+1] - x[i]
    xs = linspace(x[i],x[i+1],10)
    ys = ((x[i+1]-xs)**3*m[i] + (xs-x[i])**3*m[i+1])/(6*h) \
         + ((x[i+1]-xs)*y[i] + (xs-x[i])*y[i+1])/h \
         - ((x[i+1]-xs)*m[i] + (xs-x[i])*m[i+1])*h/6
    plot(xs,ys,'r-')

title("Cubic spline interpolation with "+str(N+1)+" points")
legend(('Exact','Data','Spline'));

We see that spline gives a better approximation and does not have too much oscillation.

6.10.4Optimality property¶

The complete cubic spline gives approximations which have less oscillations compared to other types of approximation. This is shown in the Python example. We first prove the optimality result.

Proof 1

Let

k(x) = s_c(x) - g(x)

(30)

Lets compute the oscillations in $g$ .

\begin{aligned} \int_a^b |g''(x)|^2 \ud x &= \int_a^b [s_c''(x) - k''(x)]^2 \ud x \\ &= \int_a^b [s_c''(x)]^2 \ud x - 2 \int_a^b s_c''(x) k''(x) \ud x + \int_a^b [k''(x)]^2 \ud x \end{aligned}

(31)

The second term is

\int_a^b s_c''(x) k''(x) \ud x = \sum_{i=0}^{N-1} \int_{x_i}^{x_{i+1}} s_c''(x) k''(x) \ud x

(32)

Perform integration by parts twice to remove derivative from $k$

\begin{aligned} \int_{x_i}^{x_{i+1}} s_c''(x) k''(x) \ud x &= [s_c'' k']_{x_i}^{x_{i+1}} - \int_{x_i} ^{x_{i+1}} s_c'''(x) k'(x) \ud x \\ &= [s_c'' k']_{x_i}^{x_{i+1}} - [s_c''' k]_{x_i}^{x_{i+1}} + \int_{x_i}^{x_{i+1}} s_c''''(x) k(x) \ud x \end{aligned}

(33)

But $s_c''''=0$ since $s_c$ is a cubic polynomial and

k(x_i) = s_c(x_i) - g(x_i) = f(x_i) - f(x_i) = 0, \qquad 0 \le i \le N

(34)

Hence

\int_a^b s_c''(x) k''(x) \ud x = \sum_{i=0}^{N-1} [s_c'' k']_{x_i}^{x_{i+1}}

(35)

Both $s_c''$ and $k'$ are continuous at $x_i$ , $1 \le i \le N-1$ . Hence we get

\int_a^b s_c''(x) k''(x) \ud x = -s_c''(x_0) k'(x_0) + s_c''(x_N) k'(x_N)

(36)

But

k'(x_0) = s_c'(x_0) - g'(x_0) = f'(x_0) - f'(x_0) = 0

(37)

and similarly $k'(x_N) = 0$ , so that

\int_a^b s_c''(x) k''(x) \ud x = 0

(38)

This leads to

\int_a^b |g''(x)|^2 \ud x = \int_a^b [s_c''(x)]^2 \ud x + \int_a^b [s_c''(x) - g''(x)]^2 \ud x \le \int_a^b [s_c''(x)]^2 \ud x

(39)

We get equality iff

s_c''(x) - g''(x) = 0, \qquad x \in [a,b]

(40)

i.e., if $s_c(x) - g(x)$ is affine wrt $x$ . The interpolation conditions then imply that $s_c(x) \equiv g(x)$ .

6.10.5Not-a-knot condition¶

If the end derivative conditions $f'(a)$ , $f'(b)$ are not known, then we need other conditions to complete the system of equations. We can demand that $s'''(x)$ is continuous at $x_1$ and $x_{N-1}$ . This is equivalent to requiring that $s(x)$ be a cubic spline function with knots

\{ x_0, x_2, x_3, \ldots, x_{N-3}, x_{N-2}, x_N\}

(42)

while still interpolating at all node points

\{x_0,x_1,x_2,\ldots,x_{N-1},x_N\}

(43)

This again leads to a tri-diagonal linear system of equations.

6.10.6Natural cubic spline¶

In the natural cubic spline approximation, the two end-point conditions

s''(x_0) = s''(x_N) = 0

(44)

are used. This leads to a tridiagonal system of equations. The cubic spline approximations converge slowly near the end-points. The natural cubic spline also has an optimality property, for proof, see the problems collection.