MATH661

Lecture 16: Piecewise Interpolation

1.Splines

Instead of adopting basis functions defined over the entire sampling interval $[x_{0}, x_{n}]$ as exemplified by the monomial or Lagrange bases, approximations of $f : ℝ \to ℝ$ can be constructed with different branches over each subinterval, by introducing $S_{i} : [x_{i - 1}, x_{i}] \to ℝ$ , and the approximation

p (t) = {\begin{cases} S_{1} (t) & x_{0} ⩽ t < x_{1} \\ S_{2} (t) & x_{1} ⩽ t < x_{2} \\ ⋮ & ⋮ \\ S_{n} (t) & x_{n - 1} ⩽ t < x_{n} \\ S_{n + 1} (t) & t = x_{n} \end{cases} . .

The interpolation conditions $p (x_{i}) = y_{i}$ lead to constraints

S_{i} (x_{i - 1}) = y_{i - 1} .

The form of $S (t)$ can be freely chosen, and though most often $S (t)$ is a low-degree polynomial, the spline functions may have any convenient form, e.g., trigonometric or arcs of circle. The accuracy of the $p (t)$ approximant is determined by the choice of form of $S (t)$ , and by the sample points. It is useful to introduce a quantitative measure of the sampling through the following definitions.

Definition. ${x_{0}, x_{1}, \dots, x_{n}}$ is a partition of the interval $[a, b] \subset ℝ$ if $x_{i} \in ℝ$ , $i = 0, 1, \dots, n$ , satisfy

$a = x_{0} < x_{1} < \dots < x_{n - 1} < x_{n} = b .$

Definition. The norm of partition $X = {x_{0}, x_{1}, \dots, x_{n}}$ of the interval $[a, b] \subset ℝ$ is

$|| X || = {max}_{1 ⩽ i ⩽ n} | x_{i} - x_{i - 1} | .$

Constant splines (degree 0).

A simple example is given by the constant functions

S_{i} (t) = y_{i - 1}

. Arbitrary accuracy of the approximation can be achieved in the limit of

n \to \infty

|| X || \to 0

. Over each subinterval the polynomial error formula gives

f (t) - S_{i} (t) = f^{'} (ξ_{t}) (t - x_{i - 1}),

so overall

| f (t) - p (t) | ⩽ {|| f^{'} ||}_{\infty} || X ||,

which becomes

| f (t) - p (t) | ⩽ {|| f^{'} ||}_{\infty} h,

for equidistant partitions $x_{i} = x_{0} + i h$ , $h = (x_{n} - x_{0}) / n$ . The interpolant $p (t)$ converges to $f (t)$ linearly (order of convergence is 1)

Linear splines (degree 1).

A piecewise linear interpolant is obtained by

S_{i} (t) = \frac{t - x_{i - 1}}{x_{i} - x_{i - 1}} (y_{i} - y_{i - 1}) + y_{i - 1} .

The interpolation error is bounded by

| f (t) - p (t) | ⩽ \frac{1}{2} {|| f^{'} ||}_{\infty} h^{2},

for an equidistant partition, exhibiting quadratic convergence.

Quadratic splines (degree 2).

A piecewise quadratic interpolant is formulated as

S_{i} (t) = b_{i} {(t - x_{i - 1})}^{2} + c_{i} (t - x_{i - 1}) + y_{i - 1} .

The interpolation conditions are met since $S_{i} (x_{i - 1}) = y_{i - 1}$ . The additional parameters of this higher order spline interpolant can be determined by enforcing additional conditions, typically continuity of function and derivative at the boundary between two subintervals

\begin{array}{ll} S_{i} (x_{i}) = b_{i} h_{i}^{2} + c_{i} h_{i} = y_{i}, & i = 1, 2, \dots, n \\ S_{i}^{'} (x_{i}) = 2 b_{i} h_{i} + c_{i} = 2 b_{i + 1} h_{i + 1} + c_{i + 1} = S_{i + 1}^{'} (x_{i}) & i = 1, 2, \dots, n - 1 \end{array} .

An additional condition is required to close the system, for example $S_{n}^{'} (x_{i}) = y_{n}^{'}$ (known end slope), or $S_{n}^{'} (x_{i}) = 0$ (zero end slope), or $S_{n}^{'} (x_{i}) = S_{n}^{'} (x_{i - 1})$ (constant end-slope). The coefficients $b_{i}, c_{i}$ are conveniently determined by observing that $S_{i}^{'} (t)$ is linear over interval $[x_{i - 1}, x_{i}]$ of length $h_{i} = x_{i} - x_{i - 1}$ , and is given by

S_{i}^{'} (t) = \frac{t - x_{i - 1}}{h_{i}} (s_{i} - s_{i - 1}) + s_{i - 1} = \frac{s_{i - 1}}{h_{i}} (x_{i} - t) + \frac{s_{i}}{h_{i}} (t - x_{i - 1}),

with $s_{i} = y_{i}^{'}$ , the slope of the interpolant at $x_{i} .$ The continuity of first derivative conditions $S_{i}^{'} (x_{i}) = S_{i + 1}^{'} (x_{i})$ are satisfied, and integration gives

S_{i} (t) = \frac{s_{i}}{2 h_{i}} {(t - x_{i - 1})}^{2} - \frac{s_{i - 1}}{2 h_{i}} {(x_{i} - t)}^{2} + A_{i} .

The interpolation condition $S_{i} (x_{i - 1}) = y_{i - 1}$ , determines the constant of integration $A_{i}$

A_{i} - \frac{s_{i - 1} h_{i}}{2} = y_{i - 1} \Rightarrow A_{i} = y_{i - 1} + \frac{s_{i - 1} h_{i}}{2},

Imposing the continuity of function condition $S_{i} (x_{i}) = S_{i + 1} (x_{i})$ gives

\frac{s_{i} h_{i}}{2} + y_{i - 1} + \frac{s_{i - 1} h_{i}}{2} = - \frac{s_{i} h_{i + 1}}{2} + y_{i} + \frac{s_{i} h_{i + 1}}{2},

s_{i - 1} + s_{i} = \frac{2}{h_{i}} (y_{i} - y_{i - 1}), i = 1, 2, \dots, n,

a bidiagonal system for the slopes that is solved by backward substituion in $𝒪 (2 n)$ operations. For $i = 1$ , the $s_{0}$ value arising in the system has to be given by an end condition, and the overall system $𝑩 𝒔 = 𝒅$ is defined by

𝑩 = [\begin{array}{lllll} 1 \\ 1 & 1 \\ 1 & 1 \\ ⋱ & ⋱ \\ 1 & 1 \end{array}], 𝒅 = [\begin{array}{l} \frac{2}{h_{1}} (y_{1} - y_{0}) - s_{0} \\ \frac{2}{h_{2}} (y_{2} - y_{1}) \\ ⋮ \\ \frac{2}{h_{n}} (y_{n} - y_{n - 1}) \end{array}], 𝒔 \in ℝ^{n}, 𝑩 \in ℝ^{n \times n} .

The interpolation error is bounded by

| f (t) - p (t) | ⩽ \frac{1}{2} {|| f^{'} ||}_{\infty} h^{2},

for an equidistant partition, exhibiting quadratic convergence.

Cubic splines (degree 3).

The approach outlined above can be extended to cubic splines, of special interest since continuity of curvature is achieved at the nodes, a desirable feature in many applications. The second derivative is linear

S_{i}^{''} (t) = \frac{z_{i - 1}}{h_{i}} (x_{i} - t) + \frac{z_{i}}{h_{i}} (t - x_{i - 1}),

with $z_{i - 1} = S_{i}^{''} (x_{i - 1})$ , $z_{i} = S_{i}^{''} (x_{i})$ the curvature at the endpoints of the $[x_{i - 1}, x_{i}]$ subinterval. Double integration gives

S_{i} (t) = \frac{z_{i - 1}}{6 h_{i}} {(x_{i} - t)}^{3} + \frac{z_{i}}{6 h_{i}} {(t - x_{i - 1})}^{3} + A_{i} (t - x_{i - 1}) + B_{i} (x_{i} - t) .

The interpolation conditions $S_{i} (x_{i - 1}) = y_{i - 1}$ , $S_{i - 1} (x_{i}) = y_{i}$ , gives the integration constants

A_{i} = \frac{y_{i}}{h_{i}} - \frac{z_{i} h_{i}}{6}, B_{i} = \frac{y_{i - 1}}{h_{i}} - \frac{z_{i - 1} h_{i}}{6}

and continuity of first derivative, $S_{i}^{'} (x_{i}) = S_{i + 1}^{'} (x_{i})$ , subsequently leads to a tridiagonal system for the curvatures

h_{i} z_{i - 1} + 2 (h_{i} + h_{i - 1}) z_{i} + h_{i + 1} z_{i + 1} = \frac{6 (y_{i + 1} - y_{i})}{h_{i + 1}} - \frac{6 (y_{i} - y_{i - 1})}{h_{i}}, i = 1, 2, \dots, n - 1 .

End conditions are required to close the system. Common choices include:

Zero end-curvature, also known as the natural end conditions: $z_{0} = z_{n} = 0$ .
Curvature extrapolation: $z_{0} = z_{1}$ , $z_{n} = z_{n - 1}$
Analytical end conditions given by the function curvature: $z_{0} = f^{''} (x_{0})$ , $z_{n} = f^{''} (x_{n})$ .

2. $B$ -splines

The above analytical approach becomes increasingly unwieldy for higher degree piecewise polynomials. An alternative approach is to systematically generate basis sets of desired polynomial degree over each subinterval. The starting point in this basis-spline ( $B - spline$ ) approach is the piecewise constant functions

B_{j, 0} (t) = {\begin{cases} 1 & x_{j} ⩽ t < x_{j + 1} \\ 0 & otherwise \end{cases} .,

leading to the interpolant

f (t) ≅ p (t) = \sum_{j = 0}^{n} y_{j} B_{j, 0} (t),

(1)

of $f : ℝ \to ℝ$ , as sampled by data set $𝒟 = {(x_{i}, y_{i} = f (x_{i})), i = 0, 1, \dots, n}$ , $a = x_{0} < x_{1} < \dots < x_{n} = b$ . The set

ℬ_{0} (t; 𝒙) = {B_{0, 0} (t), B_{1, 0} (t), \dots, B_{n, 0} (t)}

constitutes a basis for all piecewise constant approximants of real functions on the interval $[x_{0}, x_{n}]$ . Higher degree basis sets $ℬ_{k} (t; 𝒙)$ , $k > 0$ , are defined recursively through

B_{j, k} (t) = w_{j, k} (t) B_{j, k - 1} (t) + (1 - w_{j + 1, k} (t)) B_{j + 1, k - 1} (t),

with the weight function

w_{j, k} (t) = \frac{t - x_{j}}{x_{j + k} - x_{j}} .

$\circ$

Figure 1. $B$ -spline sets $ℬ_{0}, ℬ_{1}, ℬ_{2}, ℬ_{3}, ℬ_{4}$ with $𝒙 = [\begin{array}{llllll} 0 & 1 & 2 & 3 & 4 & 5 \end{array}]$

In[3]:=

B[j_, 0, t_] := If[j <= t < j + 1, 1, 0];
B[j_, k_, t_] := B[j, k, t] = 
 (t-j)/k B[j,k-1,t] + (1-(t-j-1)/k) B[j+1,k-1,t];
bases=Table[B[j,k,t],{j,0,5},{k,0,4}]; bases[[3,1]]

$If [2 \leq t < 3, 1, 0]$

In[5]:=

plots=Plot[bases,{t,0,10},PlotRange->All];
Export[HomeDirectory[] <> "/courses/MATH661/images/Bsplines.eps", plots]

$/home/mitran/courses/MATH661/images/Bsplines.eps$

In[6]:=

As the degree $k$ increases, the support of $B_{j, k} (t)$ increases to the interval $[x_{j}, x_{j + k + 1}]$ . This is the $B$ -spline analog of the additional end conditions in traditional spline formulations, and leads to the set

ℬ_{k} (t; 𝒙) = {B_{0, k} (t), B_{1, k} (t), \dots, B_{n, k} (t)}

defining a basis for splines of degree $k$ only on a subinterval within $[x_{0}, x_{n}]$ . Consider the piecewise linear case $k = 1$ , (Fig. 2). The set $ℬ_{1}$ forms a basis for piecewise linear functions if over each subinterval $[x_{j}, x_{j + 1}]$ an arbitrary linear function $S_{1} (t)$ can be expressed as a linear combination

S_{1} (t) = a + b t = \sum_{i = 0}^{n} c_{i} B_{i, 1} (t) .

Over $[x_{j}, x_{j + 1}]$ only $B_{j - 1, 1} (t), B_{j} (t)$ are not identically zero, hence

S_{1} (t) = c_{j - 1} B_{j - 1, 1} (t) + c_{j} B_{j, 1} (t) .

For the end interval $[x_{0}, x_{1}]$ , a definition of $B_{- 1, 1} (t)$ would be required,

S_{1} (t) = c_{- 1} B_{- 1, 1} (t) + c_{0} B_{0, 1} (t),

not available within the chosen $𝒙$ data set. At the other end interval $[x_{n - 1}, x_{n}]$ ,

S_{1} (t) = c_{n - 1} B_{n - 1, 1} (t) + c_{n} B_{n, 1} (t),

invokes $B_{n, 1}$ which requires $B_{n + 1, 0} (t)$ , again not available within the chosen data set. One can either include samples outside the $[a, b]$ interval or restrict the spline domain of definition. Again, this is analogous with the treatment of end conditions in traditional splines:

Sampling outside of the $[a, b]$ range seeks additional information on the function being interpolated $f$ , as for instance imposed by the condition $S^{'} (a) = f (a)$ in traditional splines;
Restricting the definition domain corresponds to inferring information on the behavior of $f$ in the end intervals as in the condition $S^{'} (x_{0}) = S^{'} (x_{1})$ in traditional splines.

Denote by $𝒮_{k} (t; 𝒙)$ the set of splines $S : [x_{0}, x_{n}] \to ℝ$ , that are piecewise polynomials of degree $k$ on the partition $𝒙$ of $[x_{0}, x_{n}]$ . The $k = 0$ , piecewise constant interpolant (1) is specified by $n + 1$ coefficients, the components of $𝒚 \in ℝ^{n + 1}$ , hence

\dim 𝒮_{0} (t; 𝒙) = n + 1,

i.e., the dimension of the space of piecewise-constant splines is equal to the number of sample points. As the degree $k$ increases, additional end conditions are required to specify a spline interpolation and

\dim 𝒮_{k} (t; 𝒙) = n + 1 + k,

requiring a basis set

ℬ_{.} (t; 𝒙) = {B_{- k, k} (t), \dots, B_{0, k} (t), B_{1, k} (t), \dots, B_{n, k} (t)} .

$\circ$

Figure 2. $B$ -spline set $ℬ_{1}$ for $𝒙 = [\begin{array}{llllll} 0 & 1 & 2 & 3 & 4 & 5 \end{array}]$

In[8]:=

plots=Plot[bases[[;;,2]],{t,0,10},PlotRange->All,GridLines->Automatic];
Export[HomeDirectory[] <> "/courses/MATH661/images/B1splines.eps", plots]

$/home/mitran/courses/MATH661/images/B1splines.eps$

In[9]:=

$\circ$

Algorithm B-spline evaluation (inefficient, does not account for known zero values of $B$ )

Input: $K \in ℕ, 𝒕 \in ℝ^{m}$ , $𝒙 \in ℝ^{k + n + 1}$

$𝑩 = 𝟎 \in ℝ^{m \times (k + n + 1)}$

for $i = 1 : m$

for $j = 1 : k + n$

if $x_{j} ⩽ t_{i} < x_{j + 1}$ then $B [i, j] = 1$ end

end

if $t_{i} \approx x_{k + n + 1}$ then $B [i, k + n + 1] = 1$ end

end

for $k = 1 : K$

for $j = 1 : k + n$

$w = (t - x_{j}) / (x_{j + k} - x_{j})$

$B [:, j] = w B [:, j] + (1 - w) B [:, j + 1]$

end

return $𝑩$

∴

function Bspline(t,x,K)
  np1=length(x); n=np1-1; m=length(t);B=zeros(m,np1)
  for i=1:m
    for j=1:n
      if (x[j]<=t[i]<x[j+1]) B[i,j]=1.0; end
    end
    if (t[i]≈x[np1]) B[i,np1]=1.0; end    
  end
  for k=1:K
    for j=1:n-k      
      w0 = (t .- x[j])/(x[j+k]-x[j])
      w1 = 1 .- (t .- x[j+1])/(x[j+1+k]-x[j+1])
      B[:,j] = w0 .* B[:,j] + w1 .* B[:,j+1]
    end
  end
  return B
end;

∴	n=10; k=0; x=collect(0:n); h=0.1; t=collect(0:h:n); B=Bspline(t,x,k); clf(); for j=1:n-k plot(t,B[:,j],"o"); end; grid("on")

∴	n=10; k=1; x=collect(-k:n+k); t=collect(0:h:n); B=Bspline(t,x,k); clf(); for j=1:n+k plot(t,B[:,j]); end; grid("on")

∴	n=10; k=2; x=collect(-k:n+k); t=collect(0:h:n); B=Bspline(t,x,k); clf(); for j=1:n+k plot(t,B[:,j]); end; grid("on")

∴	n=10; k=3; x=collect(-k:n+k); t=collect(0:h:n); B=Bspline(t,x,k); clf(); for j=1:n+k plot(t,B[:,j]); end; grid("on")

∴	n=10; k=3; x=collect(-k:n+k); t=collect(-k:h:n+k); B=Bspline(t,x,k); clf(); for j=1:n+k plot(t,B[:,j]); end; grid("on")

∴	n=10; k=2; x=collect(-k:n+k); t=collect(-k:h:n+k); B=Bspline(t,x,k); clf(); for j=1:n+k plot(t,B[:,j]); end; grid("on")

∴	n=10; k=4; x=collect(-k:n+k); t=collect(0:h:n); B=Bspline(t,x,k); clf(); for j=1:n+k plot(t,B[:,j]); end; grid("on")

∴	n=10; k=5; x=collect(-k:n+k); t=collect(0:h:n); B=Bspline(t,x,k); clf(); for j=1:n+k plot(t,B[:,j]); end; grid("on")

∴	n=10; k=6; x=collect(-k:n+k); t=collect(0:h:n); B=Bspline(t,x,k); clf(); for j=1:n+k plot(t,B[:,j]); end; grid("on")

∴	n=10; k=7; x=collect(-k:n+k); t=collect(0:h:n); B=Bspline(t,x,k); clf(); for j=1:n+k plot(t,B[:,j]); end; grid("on")

∴

A $B$ -spline interpolant of degree $k$ is given by a linear combination of the basis set $ℬ_{k} (t; 𝒙)$

f (t) ≅ p_{k} (t) = \sum_{j = - k}^{n} c_{i} B_{j, k} (t) .

$\circ$ The interpolation conditions $y_{i} = p (x_{i})$ lead to an underdetermined linear system for $k > 0$

𝑩 𝒄 = 𝒚, 𝑩 = [\begin{array}{lllll} B_{- k, k} (𝒙) & \dots & B_{0, k} (𝒙) & \dots & B_{n, k} (𝒙) \end{array}] \in ℝ^{(n + 1) \times (k + n + 1)},

analogous to the $k$ degrees of freedom in specification of end conditions for $𝒮_{k} (𝒙)$ .

∴	n=5; k=0; x=-k:n; t=0:n; B=Bspline(t,x,k)

$[\begin{array}{cccccc} 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 1.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 1.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 \end{array}]$ (2)

∴	n=5; k=1; x=-k:n; t=0:n; B=Bspline(t,x,k)

$[\begin{array}{ccccccc} 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 1.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 1.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 \end{array}]$ (3)

∴	n=5; k=2; x=-k:n; t=0:n; B=Bspline(t,x,k)

$[\begin{array}{cccccccc} 0.5 & 0.5 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.5 & 0.5 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.5 & 0.5 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.5 & 0.5 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.5 & 1.0 & 1.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 \end{array}]$ (4)

∴	n=5; k=3; x=-k:n; t=0:n; B=round.(10000*Bspline(t,x,k))/10000

$[\begin{array}{ccccccccc} 0.1667 & 0.6667 & 0.1667 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.1667 & 0.6667 & 0.1667 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.1667 & 0.6667 & 0.1667 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.1667 & 0.6667 & 0.5 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.1667 & 0.5 & 1.0 & 1.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 \end{array}]$ (5)

∴

rank(B)

$6$

∴

Lecture 16: Piecewise Interpolation

1.Splines

Constant splines (degree 0).

Linear splines (degree 1).

Quadratic splines (degree 2).

Cubic splines (degree 3).

2.B-splines

2. $B$ -splines