Spline Functions

Edward Neuman
Department of Mathematics
Southern Illinois University at Carbondale

Introduction

Package SplineFunctions.m  consits of several functions for computations with the univariate and multivariate splines. Topics covered include B-splines and their applications, interpolation by splines of arbitrary degree, convex interpolating splines, and quadratic histosplines. Also, a function for evaluation of  the  s-variate centered box splines on the s+1-directional meshes is included.
For the reader's convenience we give some basic definitions in the subsequent sections of this 
notebook. Below we provide definition of the polynomial splines and clarify some terms used both in this notebook and the package as well.

Definition.  Let knts = {(Subscript[t, i])}_(i = 1)^m be a sequence of real numbers that are sorted in the nondecreasing order. A spline function of degree k (k ≥ 0) with the smoothness ν (ν ≥ -1) and knots in the list knts is a piecewise polynomial function s(x) such that

                                          s(x) = (Subscript[p, i])(x)       for  (Subscript[t, i]) ≤ x ≤  (Subscript[t, i + 1])  ((Subscript[t, i]) < (Subscript[t, i + 1])),
                                    
where (Subscript[p, i]) is a polynomial of degree at most k and s(x) is ν-times continuously differentiable function on the interval ((Subscript[t, i]), (Subscript[t, i + 1])). 
  
Throughout the sequel we will refer to this representation as the pp-form. In this package the pp-form is realized as a list of lists:

                                           {{(Subscript[p, 1]), {(Subscript[t, 1]) , (Subscript[t, 2]) }}, {(Subscript[p, 2]), {(Subscript[t, 2]) , (Subscript[t, 3])}}, ... , }.
                                           
If two consecutive knots coincide, then the corresponding sublist is not included in this representation. Another form used for representing splines is called a B-form. For the spline s(x) expressed as a linear combination of the B-splines
 
                                                      s(x) = (∑_(p = i)^j) (Subscript[c, p]) (Subscript[B, p, k]) (x) 
                                                                              
its B-form is realized as a list of coefficients {(Subscript[((Subscript[c, p])}), p = i]^j) .

B-splines and Their Applications

Definition.The univariate B-spline(Subscript[B, p, k, ])   of degree k (k ≥ 0),can be realized as thelocally supported piecewise polynomial kernel of the distribution

\(\((1)\)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ f\) → k! (∫_((S^(k + 1 ))))f(  
                    (∑_(l = 0)^(k + 1)) 
                      (Subscript[λ, l]) (Subscript[t, p + l]) ) 
                  dλ,    f ϵ (Subscript[ℂ, 0]^∞)(ℝ).

Here (S^(k + 1)) stands for the k+1 - dimensional Euclidean simplex in (ℝ^(k + 1)), 𝕕λ = d(Subscript[λ, 1]) ... d(Subscript[λ, k + 1]), 
 (Subscript[λ, 0]) = 1 - (Subscript[λ, 1]) - ... - (Subscript[λ, k + 1]),  and (Subscript[t, p]) ≤ ... ≤ (Subscript[t, p + k + 1])  are the knots of the B-spline. 
 
 For the excellent exposition of these functions the reader is referred to: Carl de Boor, A Practical Guide to Splines, Springer-Verlag, 1978.
The following result (due to C. de Boor) can be used for numerically stable evaluation of B-splines:

(2)                        (Subscript[B, p, k])(x) = (((x - (Subscript[t, p]))(Subscript[B, p, k - 1])(x) + ((Subscript[t, p + k + 1]) - x)(Subscript[B, p + 1, k - 1]) (x) )/((Subscript[t, p + k + 1]) -  (Subscript[t, p]))),          k  ≥  1.

Example  1

Needs["SplineFunctions`"]

In this example we will use, on several occasions, a family of the quadratic B-splines introduced in the de Boor's book (pp.111-112):

\(k = 2; \ 
t = {0, 1, 1, 3, 4, 6, 6, 6}; 
tp = Partition[t, k + 2, 1]; 
    
le\  = \ Length[t] - k - 1; \)

\(pprep = Map[bpol[1, k, #, x]&, \ tp]; 
Map[plotpdf[#, x]&, pprep]; \)

The piecewise polynomial representations of these splines are shown below

pprep // MatrixForm

The derivatives of the above splines can be evaluated with the aid of the function diffpdf. For instance, the first order derivatives of the B-splines in the list  pprep are:

\(\(Map[diffpdf[\ pprep[\([#]\)], \ x, \ 1]&, \ Range[le]] // MatrixForm

    \)\)

The antiderivative of a B-spline is calculated using function antiderbsp. For the third B-spline of this example, we have

\(\(
ad = antiderbsp[{1, 3, 4, 6}, x]\)\)

plotpdf[ad, x];

Differentiating  we get back the B-spline function:

diffpdf[ad, x, 1] // MatrixForm

This illustrates use of the function bopbs. The square roots of the B-splines in the list pprep are calculated and next their graphs are plotted:

\(bsq = bopbs[pprep, {1/2, 1/2, 1/2, 1/2, 1/2}, x, Power]; 

    Map[plotpdf[#, x]&, bsq]; \)

It is well known that

\(TraditionalForm
    \`\((3)\)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        \(∫\_\(-∞\)\%∞\(\( B\_\(p, k\)\)(x)\) 
            dx\) = 1\/\(k + 1\)\)

For the quadratic B-splines of this example, we have

\(\(Map[integratepdf[\ pprep[\([#]\)], \ x]&, \ Range[le]]
\)\)

Moments (Subscript[μ, m]) of the B - splines are defined as follows

(4)                                                       (Subscript[μ, m])=(∫_(-∞)^∞) (x^m) 
                (Subscript[B, p, k])(x) dx,                 m = 0,1,....

Function momentbspl calculates  moments  of  the B-splines using  a method presented in:  E. Neuman, Moments and Fourier transform of B-splines, J.Comput. Appl. Math. 7(1981),  pp.51- 62.
This method does not require an explicit  knowledge of the B- spline (only its knots are needed).  For the family of B-splines used  here moments of orders one and two are calculated:

\(\(
knts = Partition[t, 4, 1]; 

    mb = Map[{momentbspl[1, #], momentbspl[2, #]}&, \ knts]\)\)

The integral inner products of the B- splines are needed when solving the least squares approximation problem from the spline spaces. We use the function innerprod to calculate
the entries of the  the Gram matrix for the family of the B-splines introduced earlier:

\(\(Table[innerprod[pprep[\([p]\)], \ pprep[\([q]\)], \ x], \ {p, 5}, 
        \ {q, 5}] // MatrixForm
\)\)

The Fourier transform of the B-splines can be found using function Fouriertbspl. For the 
B-splines used in this example, we have

Ftbspl = Map[Fouriertbspl[#, v]&, knts] // Together

The graphs of the moduli of functions in the list Ftbspl are shown below

\(\(
\(Plot[Evaluate[Abs[Ftbspl]], {v, \(-1\), 1}]; \)\)\)

Example  2

The normalized B-splines  (Subscript[N, p, k])  are defined as follows :

(5)                               (Subscript[N, p, k]) = ((Subscript[t, p + k + 1]) - (Subscript[t, p]))(Subscript[B, p, k]).

They possess a remarkable property called a partition of unity:

(∑_(p = -∞)^∞) (Subscript[N, p, k])(x) = 1     for each   x ϵ ℝ.

For a given value of x the above infinite series simplifies to a finite sum. For the family of the quadratic B-splines discussed in Example 1 we construct the corresponding normalized B-splines and next we plot the graph of their sum.

\(c = Table[t[\([p + k + 1]\)] - t[\([p]\)], \ {p, 5}]; 

    s = linearcomb[pprep, c, x]; \)

plotpdf[Rest[s],  x];

Example  3

A common application of splines is to interpolate a function or set of points in the plane. Function interpspl[ knts, xi, yi, x ] calculates a spline interpolant of degree k = n - m + 1, where m is the length of the list of knots (knts) of the interpolating spline and n is the length of the list xi housing the nodes of interpolation. If the node, say (Subscript[x, i]), is repeated ν-times, then the interpolating spline matches consecutive derivatives up to order ν-1. The corresponding numbers in the list yi are treated as the values of derivatives at (Subscript[x, i]). Let

\(knts = {0, 1}; \ xi = {0, 0, 0, 0}; \ yi = Exp[xi]; 

    interpspl[knts, xi, yi, x]\)

Resulting interpolant is a partial sum of the Maclaurin series of the exponential function (e^x).

We will now deal with the Subbotin-Marsden interpolation scheme. In this method the nodes
of interpolation are located at the endpoints of the interval of interpolation and at the 
midpoints of the intervals generated by the sequence of knots of the spline interpolant. Let

\(knts = Table[pt, {pt, 0, 1, .25}]; 

    xi = Apply[Plus, Partition[knts, 2, 1], 1]/2; 

    xi = Prepend[Append[xi, 1], 0]\)

\(spl = interpspl[knts, xi, Exp[xi], x]; 
spl // MatrixForm\)

The error of interpolation is plotted below:

\(gr1 = bopbs[spl, {\(-Exp[x]\)}, x, Plus]; 
plotpdf[gr1, x]; \)

We will interpolate now the same data using a complete cubic spline interpolant.

\(xi = {0, 0, .25, .5, .75, 1, 1}; 

    spl = interpspl[knts, xi, Exp[xi], x]; 
spl // MatrixForm\)

The error function is displayed below:

\(gr2 = bopbs[spl, {\(-Exp[x]\)}, x, Plus]; 
plotpdf[gr2, x]; \)

Example  4

In this example we will find an approximation of the Fourier transform of the function

f(x) = (e^(-(x^2))),     -1 ≤ x ≤ 1,   and  f(x) = 0, otherwise.

First, a complete cubic spline interpolating f(x) at five equally spaced points is found and next its Fourier transform, denoted by ft, is determined using function approxFt

Clear[t,  x]SetAttributes[f, Listable]\(f[x_] := Exp[\(-x^2\)]; 
knts = Table[pt, {pt, \(-1\), 1, .5}]; 

    xi = {\(-1\), \(-1\), \(-.5\), 0, .5, 1, 1}; 

    yi = {f[\(-1\)], \(f'\)[\(-1\)], f[{\(-.5\), 0, .5, 1}], \(f'\)[1]} // 
        Flatten; 
coeff = interpspl[knts, xi, yi, x, True]; 

    ft = approxFT[knts, Length[xi], coeff, t] // Together\)

The Fourier transform of f(x), denoted by fe,  is:

fe = Integrate[f[x] Exp[-I t x], {x, -1, 1}]

Plot of the modulus of the difference of ft-fe is shown below

\(\(diff = 
      If[t == 0, Abs[ft[\([2]\)] - fe /. t -> 0], 
		
        Abs[ft[\([1]\)] - fe]]; 

    Plot[diff, {t, 0 - 1, 1}, PlotRange -> All]; 
\)\)

Convex Interpolating Splines

In some problems that arise in science and engineering we often seek an approximant that preserves the shape of the data. For instance, monotonicity and/or convexity of the approximant are frequently  required in the data fitting and mapmaking. Function  ConvexSpl is used for generating a quadratic spline interpolant s(x) that preserves convexity
of the data. The interpolated data must be s-convex. The set {((Subscript[x, 1]) , (Subscript[z, 1])),  ...  , ((Subscript[x, n]),  (Subscript[z, n]))} is said to be s-convex if the associated slopes are in the ascending order. Quadratic spline generated by this function satisfies the interpolatory conditions
                                                 s((Subscript[x, i])) = (Subscript[z, i])  (i = 1, 2, ... , n)
and is continuously differentiable function on the real line ℝ. The additional knots (Subscript[ξ, i]) (i=1,2, ... , n-1) are inserted so that (Subscript[x, i]) < (Subscript[ξ, i]) < (Subscript[x, i + 1]). The algorithm used here is discussed in: J.W. Schmidt and M. Walther, Tensor Product Splines on Refined Grids in s-convex Interpolation, in: Multivariate Approximation: Recent Trends and Results, pp.189-202, Berlin, 1997.

Example  5

In this example we interpolate the exponential function (e^(-x)) by the convex quadratic spline interpolant s(x),  0 ≤ x ≤ 0.5.

\(xi = Range[0, .5, .1]; 
sp = ConvexSpl[xi, Exp[\(-xi\)], x]; 

    sp // MatrixForm\)

plotpdf[ sp, x ];

The error function s(x) - (e^(-x)) is plotted below

\(gr1 = bopbs[sp, {\(-Exp[\(-x\)]\)}, x, Plus]; 
plotpdf[gr1, x]; \)

If the data are not s-convex, then the function ConvexSpl does not return anything and the warning message is issued:

\(xi = Range[\(-1\), 1, .25]; 
ConvexSpl[xi, xi^3, x]\)

Quadratic Histosplines

Definition. A continuously differentiable quadratic spline s(x), with the knot set
                                                   x = {(Subscript[x, 1]) < (Subscript[x, 2]) < ... < (Subscript[x, n])}
 is said to be a quadratic histospline if it satisfies the following interpolatory conditions:

(1)                                   (Subscript[(Subscript[(∫_x), i]^x), i + 1])s(x) dx = (Subscript[F, i]),          1  ≤  i  ≤  n  -  1,

(2)                                    s((Subscript[x, 1])) = (Subscript[y, 1]),   s((Subscript[x, n])) = (Subscript[y, n]).

Here (Subscript[F, 1]), ... , (Subscript[F, n]), (Subscript[y, 1]), and (Subscript[y, n]) are given real numbers. For a derivation of the algorithm for computing quadratic histosplines the reader is referred to: A. Sharma and J. Tzimbalario, Quadratic splines, J. Approx. Th. 19(1977), pp. 186-193.
The histospline s(x) can be used to fit a histogram. Function QuadraticHistospline generates the interpolant s(x) in the pp-form.

Example 6

\(xi\  = \ {\(-1\), 0, 2, 3, 5, 6}; 
F = {3, \ 2, \ 4, \(-\ 1\), \ 10}; 
    
y1\  = \ 1; \ yn\  = \ 10; 
k = Length[xi]; \)

hispl   =  QuadraticHistospline[xi,  F,  y1,  yn,  x]

We will check the interpolatory conditions (1) and (2)

Map[integratepdf[{#}, x]&, hispl] - F

{y1 - hispl[[1, 1]] /. x -> xi[[1]], 
      yn - hispl[[k - 1, 1]] /. x -> xi[[k]]}

Next, we will check continuity of s(x) and s'(x) at the interior knots

\(\(cnt = 
      Table[
	Table[{
            Limit[D[hispl[\([i - 1, 1]\)], {x, j}], x -> xi[\([i]\)], 
              Direction -> 1], 
		
            Limit[D[hispl[\([i, 1]\)], {x, j}], x -> xi[\([i]\)], 
              Direction -> \(-1\)]}, 
			{i, 2, k - 1}], {j, 0, 1}]; 
    \)\)

The left and right limits of s(x) at the interior knots:

cnt[[1]]

The corresponding limits for s'(x) are:

cnt[[2]]

This produces the graphs of the histogram (g2) and its interpolating quadratic histospline (g1)

\(\(\ \ g1 = plotpdf[hispl, x, True]; 	
\ \ 
    st = Map[UnitStep[x - #[\([1]\)], #[\([2]\)] - x]&, 
		
        Partition[xi, 2, 1]]; 
\ \ 
    g2 = Plot[F.st, {x, xi[\([1]\)] - 0.5, xi[\([k]\)] + .5}, 
        PlotStyle -> {Thickness[.01], RGBColor[0, 0, 1]}, 
		
        DisplayFunction -> Identity]; \)\)

Show[g1, g2, DisplayFunction -> $DisplayFunction];

Box Splines on s+1-Directional Meshes

Definition. The box spline (Subscript[B, n])(▫ |X) can be realized as the locally supported piecewise polynomial kernel of the distribution

(1)                                 f → (∫_((I^m))) f ((∑_(j = 1)^m) (Subscript[λ, j ])(t^j)) 𝕕λ,                 f ϵ (Subscript[ℂ, 0]^∞)((ℝ^s)).

Here I = [((-1)/2), (1/2)] ,   
        (t^j) ϵ (ℝ^s) 
          ( 1 ≤  j   ≤  m) are the columns 
          of the matrix  X , and 𝕕λ=d(Subscript[λ, 1]) ... d(Subscript[λ, m]).
(See:C.de Boor and R.DeVore,Approximation by smooth multivariate splines ,
Trans.Amer.Math.Soc. 276 (1983),pp .775-788.)
		Let  (Subscript[e, j ])(1 ≤ j ≤ s) denote the 
        jth coordinate vector in (ℝ^s).Also, 
      let (Subscript[e, s + 1]) = (Subscript[e, 1]) + ... + (Subscript[e, s]).For(Subscript[n, 1]),...,(Subscript[n, s + 1]) ϵ {1, 2, ...}  let m=(Subscript[n, 1]) + ...+(Subscript[n, s + 1]).The s+1-directional mesh in  (ℝ^s) is the matrix

(2)                                                  X = ((Subscript[e, 1]), ... , (Subscript[e, s + 1])),

where each vector (Subscript[e, i]) is  repeated  (Subscript[n, i]) - 
          times. Box spline with the matrix (2) is said to be the centered box spline. Function Boxsplsdm[x, n]  evaluates the centered box splines on s + 1 - 
            directional meshes at the point x 
              ϵ(ℝ^s), 
          and the parameter n is the s+1-dimensional vector of multiplicities (Subscript[n, i]).The algorithm used here is presented in:
 Edward Neuman,A new formula for box splines on three - 
                directional meshes,Math.Comp. 62 (205),1994,pp .227-229.
 Function BoxSplsdm can also be used for computing 
                certain multiple integrals.Let 
 sinc z = (sin  z)/z denote thesincfunction. Then

\(TraditionalForm
    \`\(\((3)\)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            \(∫\_\(ℝ\^s\)\ 
              \((sinc\ t\_1)\)\^n\_1\ \((sinc\ t\_2)\)\^n\_2\)\ ...\)\ 
        \(\((sinc(\(\(t\_1 + \)\ ...\)\  + t\_s))\)\^n\_\(s + 1\)\) 
        dt\  = \ \(π\^s\) \(\(B\_n\)(0\_s | X)\), 


    where\ 𝕕t\  = \ 
      \(dt\_1\ ...\)\ dt\_s\ \ and\ \ 0\_s\ \ stands\ for\ the\ zero\ vector\ 
        in\ ℝ\^s\ \(\((Neuman, \ p  .229)\).\)\)

Example 7

Let (Subscript[n, 1]) = (Subscript[n, 2]) = (Subscript[n, 3]) = 2 and let s = 2. The corresponding box spline (on three-directional mesh) is a bivariate piecewise polynomial function of total degree four with continuous partial derivatives of order two.

tb=Table[BoxSplsdm[{x,y}, {2,2,2}], {x,-1.5,1.5,.25},{y,-1.5,1.5,.25}];

ListPlot3D[tb, ColorFunction->Hue, Boxed->False, Axes->Off];

For s = 2 and (Subscript[n, 1]) = (Subscript[n, 2]) = (Subscript[n, 3]) = 1 Mathematica cannot evaluate integral (3).

sinc[z]:=Sin[z]/z

Timing[Integrate[sinc[u] sinc[v] sinc[u+v], {v,-∞,∞},{u,-∞,∞}]]

In order to evaluate integral (3), using function BoxSplsdm,  let us denote it by J and define

J[n_List] := 
      (π^(Length[n] - 1)) BoxSplsdm[Table[0, {Length[n] - 1}], n]

Table[Timing[J[{i, i, i}]], {i, 3}] // MatrixForm

Table[Timing[J[{i, i, i, i}]], {i, 3}] // MatrixForm