#%config InlineBackend.figure_format = 'svg'
from pylab import *
from scipy.interpolate import barycentric_interpolate
Consider a grid for the interval [ a , b ] [a,b] [ a , b ]
a = x 0 < x 1 < … < x N = b a = x_0 < x_1 < \ldots < x_N = b a = x 0 < x 1 < … < x N = b We have seem piecewise polynomial approximations which lead to globally continuous approximations. If we are given { f ( x i ) , f ′ ( x i ) } \{f(x_i), f'(x_i)\} { f ( x i ) , f ′ ( x i )} we can perform piecewise Hermite interpolation leading to globally C 1 \cts^1 C 1 approximations. But this requires knowledge of f ′ f' f ′ . Alternately, spline approximations can lead to globally smooth approximations without knowing derivative information.
The derivative of a spline of order m m m is a spline of order m − 1 m-1 m − 1 , and similarly for anti- derivatives. Note that
m = 2 m=2 m = 2 corresponds to piecewise linear function which is continuous, andm = 3 m=3 m = 3 is a piecewise quadratic function which is continuous and has continuous derivatives Kincaid & Cheney, 2002 .6.10.1 Cubic spline interpolation ¶ Consider the case of m = 4 m=4 m = 4 where s ( x ) s(x) s ( x ) is a cubic polynomial in every
sub-interval [ x i − 1 , x i ] [x_{i-1},x_i] [ x i − 1 , x i ] and we are given the values of some
function
y i = f ( x i ) , i = 0 , 1 , … , N y_i = f(x_i), \qquad i=0,1,\ldots,N y i = f ( x i ) , i = 0 , 1 , … , N Let us count the number
of unknowns and number of equations available. We can write, for
i = 1 , 2 , … , N i=1,2,\ldots,N i = 1 , 2 , … , N
s ( x ) = a i + b i x + c i x 2 + d i x 3 , x ∈ [ x i − 1 , x i ] s(x) = a_i + b_i x + c_i x^2 + d_i x^3, \qquad x \in [x_{i-1},x_i] s ( x ) = a i + b i x + c i x 2 + d i x 3 , x ∈ [ x i − 1 , x i ] Hence we have 4 N 4N 4 N unknown coefficients { a i , b i , c i , d i } \{a_i,b_i,c_i,d_i\} { 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 , i = 0 , 1 , … , N s(x_i) = y_i, \qquad i=0,1,\ldots,N s ( x i ) = y i , i = 0 , 1 , … , N and its derivatives must be
continuous at the interior points
s ( r ) ( x i − ) = s ( r ) ( x i + ) , i = 1 , 2 , … , N − 1 , r = 0 , 1 , 2 s^{(r)}(x_i^-) = s^{(r)}(x_i^+), \qquad i=1,2,\ldots,N-1, \qquad r=0,1,2 s ( r ) ( x i − ) = s ( r ) ( x i + ) , i = 1 , 2 , … , N − 1 , r = 0 , 1 , 2 This gives us ( N + 1 ) + 3 ( N − 1 ) = 4 N − 2 (N+1)+3(N-1) = 4N-2 ( N + 1 ) + 3 ( N − 1 ) = 4 N − 2 equations, and we are short of two
equations since we have 4 N 4N 4 N unknown coefficients. There are different
ways in which the two missing equations can be supplied.
6.10.2 Construction of cubic spline ¶ Let
M i = s ′ ′ ( x i ) , i = 0 , 1 , … , N M_i = s''(x_i), \qquad i=0,1,\ldots,N M i = s ′′ ( x i ) , i = 0 , 1 , … , N Since s ( x ) s(x) s ( x ) is cubic on [ x i , x i + 1 ] [x_i,x_{i+1}] [ x i , x i + 1 ] , then s ′ ′ ( x ) s''(x) s ′′ ( x ) is linear and thus
s ′ ′ ( x ) = ( x i + 1 − x ) M i + ( x − x i ) M i + 1 h i , i = 0 , 1 , … , N − 1 s''(x) = \frac{(x_{i+1}-x)M_i + (x-x_i)M_{i+1}}{h_i}, \qquad i=0,1,\ldots,N-1 s ′′ ( x ) = h i ( x i + 1 − x ) M i + ( x − x i ) M i + 1 , i = 0 , 1 , … , N − 1 where
h i = x i + 1 − x i h_i = x_{i+1}-x_i h i = x i + 1 − x i By the above construction s ′ ′ ( x ) s''(x) s ′′ ( x ) is continuous on [ a , b ] [a,b] [ a , b ] . Integrating twice, we get
s ( x ) = 1 6 h 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 ) , x ∈ [ x i , x i + 1 ] \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*} s ( x ) = 6 h i 1 [( 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 ) , x ∈ [ x i , x i + 1 ] with C i , D i , M i C_i,D_i,M_i C i , D i , M i to be determined. Using the interpolation conditions s ( x i ) = y i s(x_i)=y_i s ( x i ) = y i and s ( x i + 1 ) = y i + 1 s(x_{i+1}) = y_{i+1} s ( x i + 1 ) = y i + 1 gives
C i = y i h i − h i M i 6 , D i = y i + 1 h i − h i M i + 1 6 , i = 0 , 1 , … , N − 1 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 C i = h i y i − 6 h i M i , D i = h i y i + 1 − 6 h i M i + 1 , i = 0 , 1 , … , N − 1 At this stage, the only unknowns are the M i M_i M i values. Note that by construction, s ( x ) s(x) s ( x ) and s ′ ′ ( x ) s''(x) s ′′ ( x ) are continuous on [ a , b ] [a,b] [ a , b ] . We will now enforce continuity of s ′ ( x ) s'(x) s ′ ( x ) at the interior nodes. On [ x i , x i + 1 ] [x_i,x_{i+1}] [ x i , x i + 1 ]
s ′ ( x ) = 1 2 h i [ − ( x i + 1 − x ) 2 M i + ( x − x i ) 2 M i + 1 ] + y i + 1 − y i h i − ( M i + 1 − M i ) h i 6 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} s ′ ( x ) = 2 h i 1 [ − ( x i + 1 − x ) 2 M i + ( x − x i ) 2 M i + 1 ] + h i y i + 1 − y i − 6 ( M i + 1 − M i ) h i and on [ x i − 1 , x i ] [x_{i-1},x_i] [ x i − 1 , x i ]
s ′ ( x ) = 1 2 h i − 1 [ − ( x i − x ) 2 M i − 1 + ( x − x i − 1 ) 2 M i ] + y i − y i − 1 h i − 1 − ( M i − M i − 1 ) h i − 1 6 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} s ′ ( x ) = 2 h i − 1 1 [ − ( x i − x ) 2 M i − 1 + ( x − x i − 1 ) 2 M i ] + h i − 1 y i − y i − 1 − 6 ( M i − M i − 1 ) h i − 1 so that continuity of derivative at x = x i x=x_i x = x i , s ′ ( x i − ) = s ′ ( x i + ) s'(x_i^-) = s'(x_i^+) s ′ ( x i − ) = s ′ ( x i + ) , yields
h i − 1 6 M i − 1 + h i − 1 + h i 3 M i + h i 6 M i + 1 = y i + 1 − y i h i − y i − y i − 1 h i − 1 , 1 ≤ i ≤ N − 1 \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 6 h i − 1 M i − 1 + 3 h i − 1 + h i M i + 6 h i M i + 1 = h i y i + 1 − y i − h i − 1 y i − y i − 1 , 1 ≤ i ≤ N − 1 We have N − 1 N-1 N − 1 equations for the N + 1 N+1 N + 1 unknowns M 0 , M 1 , … , M N M_0, M_1, \ldots,M_N M 0 , M 1 , … , M N .
6.10.3 End-point derivative conditions: complete cubic spline ¶ To obtain the two extra equations, we impose
s ′ ( x 0 ) = y 0 ′ , s ′ ( x N ) = y N ′ s'(x_0) = y_0', \qquad s'(x_N) = y_N' s ′ ( x 0 ) = y 0 ′ , s ′ ( x N ) = y N ′ where y 0 ′ = f ′ ( x 0 ) y_0'=f'(x_0) y 0 ′ = f ′ ( x 0 ) , y N ′ = f ′ ( x N ) y_N'=f'(x_N) y N ′ = f ′ ( x N ) are assumed to be known. This yields the equations
h 0 3 M 0 + h 0 6 M 1 = y 1 − y 0 h 0 − y 0 ′ h N − 1 6 M N − 1 + h N − 1 3 M N = y N ′ − y N − y N − 1 h N − 1 \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} 3 h 0 M 0 + 6 h 0 M 1 6 h N − 1 M N − 1 + 3 h N − 1 M N = h 0 y 1 − y 0 − y 0 ′ = y N ′ − h N − 1 y N − y N − 1 We now have N + 1 N+1 N + 1 equations for the N + 1 N+1 N + 1 unknowns { M 0 , M 1 , … , M N } \{M_0, M_1, \ldots, M_N\} { M 0 , M 1 , … , M N } . They can be written in matrix form
where
b = [ y 1 − y 0 h 0 − y 0 ′ { y i + 1 − y i h i − y i − y i − 1 h i − 1 } i = 1 N − 1 y N ′ − y N − y N − 1 h N − 1 ] , M = [ M 0 M 1 ⋮ M N ] 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} b = ⎣ ⎡ h 0 y 1 − y 0 − y 0 ′ { h i y i + 1 − y i − h i − 1 y i − y i − 1 } i = 1 N − 1 y N ′ − h N − 1 y N − y N − 1 ⎦ ⎤ , M = ⎣ ⎡ M 0 M 1 ⋮ M N ⎦ ⎤ A = [ h 0 3 h 0 6 0 ⋯ diag { h i − 1 6 , h i − 1 + h i 3 h i 6 } i = 1 N − 1 ⋯ 0 h N − 1 6 h N − 1 3 ] 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] A = ⎣ ⎡ 3 h 0 6 h 0 0 ⋯ diag { 6 h i − 1 , 3 h i − 1 + h i 6 h i } i = 1 N − 1 ⋯ 0 6 h N − 1 3 h N − 1 ⎦ ⎤ The matrix A A 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 8 N 8N 8 N arithmetic operations (Thomas Algorithm). The resulting cubic spline is called the complete cubic spline interpolant, and we denote it by s c ( x ) s_c(x) s c ( x ) .
Theorem 1 (Error of complete cubic spline)
Let f ∈ C 4 [ a , b ] f \in \cts^4[a,b] f ∈ C 4 [ a , b ] and given a partition
τ N : a = x 0 < x 1 < … < x N = b \tau_N: \quad a=x_0 < x_1 < \ldots < x_N = b τ N : a = x 0 < x 1 < … < x N = b define
∣ τ N ∣ : = max 1 ≤ i ≤ N ( x i − x i − 1 ) |\tau_N| := \max_{1 \le i \le N} (x_i - x_{i-1}) ∣ τ N ∣ := 1 ≤ i ≤ N max ( x i − x i − 1 ) Let s c ( x ) s_c(x) s c ( x ) be the
complete cubic spline interpolant of f ( x ) f(x) f ( x ) on the partition τ N \tau_N τ N
s c ( x i ) = f ( x i ) , i = 0 , 1 , … , N s c ′ ( a ) = f ′ ( a ) s c ′ ( b ) = f ′ ( b ) \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} s c ( x i ) s c ′ ( a ) s c ′ ( b ) = f ( x i ) , i = 0 , 1 , … , N = f ′ ( a ) = f ′ ( b ) Then for some constants c j c_j c j
max a ≤ x ≤ b ∣ f ( j ) ( x ) − s c ( j ) ( x ) ∣ ≤ c j ∣ τ N ∣ 4 − j max a ≤ x ≤ b ∣ f ( 4 ) ( x ) ∣ , j = 0 , 1 , 2 \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 a ≤ x ≤ b max ∣ f ( j ) ( x ) − s c ( j ) ( x ) ∣ ≤ c j ∣ τ N ∣ 4 − j a ≤ x ≤ b max ∣ f ( 4 ) ( x ) ∣ , j = 0 , 1 , 2 With the additional assumption
∣ τ N ∣ min 1 ≤ i ≤ N ( x i − x i − 1 ) ≤ γ < ∞ , N → ∞ \frac{ |\tau_N|}{\min_{1 \le i \le N}(x_i - x_{i-1})} \le \gamma < \infty, \qquad N \to \infty min 1 ≤ i ≤ N ( x i − x i − 1 ) ∣ τ N ∣ ≤ γ < ∞ , N → ∞ the above bound holds for j = 3 j=3 j = 3 also. Acceptable constants are
c 0 = 5 384 , c 1 = 1 24 , c 2 = 3 8 c_0 = \frac{5}{384}, \qquad c_1 = \frac{1}{24}, \qquad c_2 = \frac{3}{8} c 0 = 384 5 , c 1 = 24 1 , c 2 = 8 3 For proof, see Boor, 1978 .
Taking j = 0 j=0 j = 0 we see that for a uniform grid of size h = ∣ τ N ∣ h = |\tau_N| h = ∣ τ N ∣
max a ≤ x ≤ b ∣ f ( x ) − s c ( x ) ∣ ≤ 5 384 h 4 max a ≤ x ≤ b ∣ f ( 4 ) ( x ) ∣ \max_{a \le x \le b}|f(x) - s_c(x)| \le \frac{5}{384} h^4 \max_{a \le x \le b} |f^{(4)} (x)| a ≤ x ≤ b max ∣ f ( x ) − s c ( x ) ∣ ≤ 384 5 h 4 a ≤ x ≤ b max ∣ f ( 4 ) ( x ) ∣ The rate of convergence is 4; if h h h is reduced to h / 2 h/2 h /2 , the errors reduces by a factor of 1 / 16 1/16 1/16 .
For uniform refinement, γ = 1 \gamma = 1 γ = 1 and third derivative s c ′ ′ ′ s_c''' s c ′′′ converges at O ( h ) \order{h} O ( h ) . Note that every time we differentiate s c s_c s c , we lose one order of convergence rate.
We will approximate the function
f ( x ) = ( 1 − x 2 ) 2 sin ( 4 π x ) exp ( sin ( 2 π x ) ) , x ∈ [ − 1 , + 1 ] f(x) = (1-x^2)^2 \sin(4\pi x) \exp(\sin(2\pi x)), \qquad x \in [-1,+1] f ( x ) = ( 1 − x 2 ) 2 sin ( 4 π x ) exp ( sin ( 2 π x )) , x ∈ [ − 1 , + 1 ] for which
f ′ ( − 1 ) = f ′ ( + 1 ) = 0 f'(-1) = f'(+1) = 0 f ′ ( − 1 ) = f ′ ( + 1 ) = 0 which will be used in constructing the cubic spline below. The function looks like this.
Sourcexmin, 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.4 Optimality 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.
Let g ∈ C 2 [ a , b ] g \in \cts^2[a,b] g ∈ C 2 [ a , b ] and satisfies the interpolation conditions
g ( x i ) = f ( x i ) , i = 0 , 1 , … , N g ′ ( x 0 ) = f ′ ( x 0 ) g ′ ( x N ) = f ′ ( x N ) \begin{aligned}
g(x_i) &= f(x_i), \quad i=0,1,\ldots,N \\
g'(x_0) &= f'(x_0) \\
g'(x_N) &= f'(x_N)
\end{aligned} g ( x i ) g ′ ( x 0 ) g ′ ( x N ) = f ( x i ) , i = 0 , 1 , … , N = f ′ ( x 0 ) = f ′ ( x N ) Then
∫ a b ∣ s c ′ ′ ( x ) ∣ 2 d x ≤ ∫ a b ∣ g ′ ′ ( x ) ∣ 2 d x \int_a^b |s_c''(x)|^2 \ud x \le \int_a^b |g''(x)|^2 \ud x ∫ a b ∣ s c ′′ ( x ) ∣ 2 d x ≤ ∫ a b ∣ g ′′ ( x ) ∣ 2 d x with
equality iff g ( x ) ≡ s c ( x ) g(x) \equiv s_c(x) g ( x ) ≡ s c ( x ) . Thus s c s_c s c oscillates least among all
smooth functions satisfying the interpolation conditions.
Let
k ( x ) = s c ( x ) − g ( x ) k(x) = s_c(x) - g(x) k ( x ) = s c ( x ) − g ( x ) Lets compute the oscillations in g g g .
∫ a b ∣ g ′ ′ ( x ) ∣ 2 d x = ∫ a b [ s c ′ ′ ( x ) − k ′ ′ ( x ) ] 2 d x = ∫ a b [ s c ′ ′ ( x ) ] 2 d x − 2 ∫ a b s c ′ ′ ( x ) k ′ ′ ( x ) d x + ∫ a b [ k ′ ′ ( x ) ] 2 d x \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} ∫ a b ∣ g ′′ ( x ) ∣ 2 d x = ∫ a b [ s c ′′ ( x ) − k ′′ ( x ) ] 2 d x = ∫ a b [ s c ′′ ( x ) ] 2 d x − 2 ∫ a b s c ′′ ( x ) k ′′ ( x ) d x + ∫ a b [ k ′′ ( x ) ] 2 d x The second term is
∫ a b s c ′ ′ ( x ) k ′ ′ ( x ) d x = ∑ i = 0 N − 1 ∫ x i x i + 1 s c ′ ′ ( x ) k ′ ′ ( x ) d x \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 ∫ a b s c ′′ ( x ) k ′′ ( x ) d x = i = 0 ∑ N − 1 ∫ x i x i + 1 s c ′′ ( x ) k ′′ ( x ) d x Perform integration by parts twice to remove derivative from k k k
∫ x i x i + 1 s c ′ ′ ( x ) k ′ ′ ( x ) d x = [ s c ′ ′ k ′ ] x i x i + 1 − ∫ x i x i + 1 s c ′ ′ ′ ( x ) k ′ ( x ) d x = [ s c ′ ′ k ′ ] x i x i + 1 − [ s c ′ ′ ′ k ] x i x i + 1 + ∫ x i x i + 1 s c ′ ′ ′ ′ ( x ) k ( x ) d x \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} ∫ x i x i + 1 s c ′′ ( x ) k ′′ ( x ) d x = [ s c ′′ k ′ ] x i x i + 1 − ∫ x i x i + 1 s c ′′′ ( x ) k ′ ( x ) d x = [ s c ′′ k ′ ] x i x i + 1 − [ s c ′′′ k ] x i x i + 1 + ∫ x i x i + 1 s c ′′′′ ( x ) k ( x ) d x But s c ′ ′ ′ ′ = 0 s_c''''=0 s c ′′′′ = 0 since s c s_c s c is a cubic polynomial and
k ( x i ) = s c ( x i ) − g ( x i ) = f ( x i ) − f ( x i ) = 0 , 0 ≤ i ≤ N k(x_i) = s_c(x_i) - g(x_i) = f(x_i) - f(x_i) = 0, \qquad 0 \le i \le N k ( x i ) = s c ( x i ) − g ( x i ) = f ( x i ) − f ( x i ) = 0 , 0 ≤ i ≤ N Hence
∫ a b s c ′ ′ ( x ) k ′ ′ ( x ) d x = ∑ i = 0 N − 1 [ s c ′ ′ k ′ ] x i x i + 1 \int_a^b s_c''(x) k''(x) \ud x = \sum_{i=0}^{N-1} [s_c'' k']_{x_i}^{x_{i+1}} ∫ a b s c ′′ ( x ) k ′′ ( x ) d x = i = 0 ∑ N − 1 [ s c ′′ k ′ ] x i x i + 1 Both s c ′ ′ s_c'' s c ′′ and k ′ k' k ′ are continuous at x i x_i x i , 1 ≤ i ≤ N − 1 1 \le i \le N-1 1 ≤ i ≤ N − 1 . Hence
we get
∫ a b s c ′ ′ ( x ) k ′ ′ ( x ) d x = − s c ′ ′ ( x 0 ) k ′ ( x 0 ) + s c ′ ′ ( x N ) k ′ ( x N ) \int_a^b s_c''(x) k''(x) \ud x = -s_c''(x_0) k'(x_0) + s_c''(x_N) k'(x_N) ∫ a b s c ′′ ( x ) k ′′ ( x ) d x = − s c ′′ ( x 0 ) k ′ ( x 0 ) + s c ′′ ( x N ) k ′ ( x N ) But
k ′ ( x 0 ) = s c ′ ( x 0 ) − g ′ ( x 0 ) = f ′ ( x 0 ) − f ′ ( x 0 ) = 0 k'(x_0) = s_c'(x_0) - g'(x_0) = f'(x_0) - f'(x_0) = 0 k ′ ( x 0 ) = s c ′ ( x 0 ) − g ′ ( x 0 ) = f ′ ( x 0 ) − f ′ ( x 0 ) = 0 and
similarly k ′ ( x N ) = 0 k'(x_N) = 0 k ′ ( x N ) = 0 , so that
∫ a b s c ′ ′ ( x ) k ′ ′ ( x ) d x = 0 \int_a^b s_c''(x) k''(x) \ud x = 0 ∫ a b s c ′′ ( x ) k ′′ ( x ) d x = 0 This leads to
∫ a b ∣ g ′ ′ ( x ) ∣ 2 d x = ∫ a b [ s c ′ ′ ( x ) ] 2 d x + ∫ a b [ s c ′ ′ ( x ) − g ′ ′ ( x ) ] 2 d x ≤ ∫ a b [ s c ′ ′ ( x ) ] 2 d x \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 ∫ a b ∣ g ′′ ( x ) ∣ 2 d x = ∫ a b [ s c ′′ ( x ) ] 2 d x + ∫ a b [ s c ′′ ( x ) − g ′′ ( x ) ] 2 d x ≤ ∫ a b [ s c ′′ ( x ) ] 2 d x We get equality iff
s c ′ ′ ( x ) − g ′ ′ ( x ) = 0 , x ∈ [ a , b ] s_c''(x) - g''(x) = 0, \qquad x \in [a,b] s c ′′ ( x ) − g ′′ ( x ) = 0 , x ∈ [ a , b ] i.e., if s c ( x ) − g ( x ) s_c(x) - g(x) s c ( x ) − g ( x ) is affine wrt x x x . The interpolation conditions then imply that s c ( x ) ≡ g ( x ) s_c(x) \equiv g(x) s c ( x ) ≡ g ( x ) .
6.10.5 Not-a-knot condition ¶ If the end derivative conditions f ′ ( a ) f'(a) f ′ ( a ) , f ′ ( b ) f'(b) f ′ ( b ) are not known, then we
need other conditions to complete the system of equations. We can demand
that s ′ ′ ′ ( x ) s'''(x) s ′′′ ( x ) is continuous at x 1 x_1 x 1 and x N − 1 x_{N-1} x N − 1 . This is equivalent
to requiring that s ( x ) s(x) s ( x ) be a cubic spline function with knots
{ x 0 , x 2 , x 3 , … , x N − 3 , x N − 2 , x N } \{ x_0, x_2, x_3, \ldots, x_{N-3}, x_{N-2}, x_N\} { x 0 , x 2 , x 3 , … , x N − 3 , x N − 2 , x N } while still interpolating at all node points
{ x 0 , x 1 , x 2 , … , x N − 1 , x N } \{x_0,x_1,x_2,\ldots,x_{N-1},x_N\} { x 0 , x 1 , x 2 , … , x N − 1 , x N } This again leads to a tri-diagonal linear system of equations.
6.10.6 Natural cubic spline ¶ In the natural cubic spline approximation, the two end-point conditions
s ′ ′ ( x 0 ) = s ′ ′ ( x N ) = 0 s''(x_0) = s''(x_N) = 0 s ′′ ( x 0 ) = s ′′ ( x N ) = 0 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.
Kincaid, D., & Cheney, W. (2002). Numerical Analysis: Mathematics of Scientific Computing (3rd ed.). American Mathematics Society. de Boor, C. (1978). A Practical Guide to Splines . Springer New York, NY.