Spline Toolbox

csape

Cubic spline interpolation with end conditions

Syntax

pp = csape(x,y)
pp = csape(x,y,conds,valconds)

Description

A cubic spline s (in ppform) with knot sequence x is constructed that satisfies s(x(j)) = y(:,j) for all j, as well as an additional end condition at the first and at the last data site, as specified by conds and valconds.

conds may be a string whose first character matches one of the following: 'complete' or 'clamped', 'not-a-knot', 'periodic', 'second', 'variational', with the following meanings.

'complete'	Match endslopes (as given in valconds, with default as under "default'')
'not-a-knot'	Make second and second-last sites inactive knots (ignoring valconds if given)
'periodic'	Match first and second derivatives at first site with those at last site
'second'	Match end second derivatives (as given in valconds, with default [0 0], i.e., as in 'variational')
'variational'	Set end second derivatives equal to zero (ignoring valconds if given)
default	Match endslopes to the slope of the cubic that matches the first four data at the respective end (i.e., Lagrange)

By giving conds as a 1-by-2 matrix instead, it is possible to specify different conditions at the two endpoints. Explicitly, Dⁱs is given the value valconds(j) at the left (j = 1) respectively right (j = 2) endpoint in case conds(j) = i,
i = 1,2. There are default values for conds and/or valconds.  

Available conditions are:   

clamped	Ds(e) = valconds()	conds() = 1
curved	D2s(e) = valconds()	conds() = 2
Lagrange	Ds(e) = Dp(e)	default
periodic	Drs(a) = Drs(b),r = 1,2	conds() = [0  0]
variational	D2s(e) = 0	conds() = 2 and valconds(j) = 0

Here, e = a(b) is the first (last) data site in case j = 1 (j = 2), and (in the Lagrange condition) p is the cubic polynomial that interpolates to the given data at e and the three sites nearest e.

If conds(j) is not specified or is different from 0, 1, or 2, then it is taken to be 1 and the corresponding valconds(j) is taken to be the corresponding default value.

The default value for valconds(j) is the derivative of the cubic interpolant at the nearest four sites in case conds(j) = 1, and is 0 otherwise.

It is possible (and, in the case of gridded data required) to specify valconds as part of y. Specifically, if size(y)==[d,ny] and ny==length(x)+2, then valconds is taken to be y(:,[1 end]), and y(:,i+1) is matched at x(i), i=1:length(x).

It is also possible to handle gridded data, by having x be a cell array containing m univariate meshes and, correspondingly, having y be an m-dimensional array (or an m+1-dimensional array if the function is to be vector-valued). Correspondingly, conds is a cell array with m entries, but the information normally specified by valconds is now expected to be part of y.

This command calls on a much expanded version of the Fortran routine CUBSPL in PGS.

Examples

csape(x,y) provides the cubic spline interpolant with the Lagrange end conditions, while csape(x,y,[2 2]) provides the variational, or natural cubic spline interpolant, as does csape(x,y,'v'). csape([-1 1],[-1 1],
[1 2],[3 6]) provides the cubic polynomial p for which p(-1), = -1, Dp(-1) = 3, p(1) = 1, D²p(1) = 6, i.e., p(x) = x³ . Finally, csape([-1 1],[-1 1]) provides the straight line p for which p(±1) = ±1, i.e., p(x) = x.

As a multivariate vector-valued example, here is a sphere, done as a parametric bicubic spline, 3d-valued, using prescribed slopes in one direction and periodic side conditions in the other:

x = 0:4; y=-2:2; s2 = 1/sqrt(2);
 clear v
 v(3,:,:) = [0 1 s2 0 -s2 -1 0].'*[1 1 1 1 1];
 v(2,:,:) = [1 0 s2 1 s2 0 -1].'*[0 1 0 -1 0];
 v(1,:,:) = [1 0 s2 1 s2 0 -1].'*[1 0 -1 0 1];
 sph = csape({x,y},v,{'clamped','periodic'});
 values = fnval(sph,{0:.1:4,-2:.1:2});
 surf(squeeze(values(1,:,:)),squeeze(values(2,:,:)),... 
  squeeze(values(3,:,:))); axis equal, axis off

The lines involving fnval and surf could have been replaced by the simple command: fnplt(sph). Note that v is a 3-dimensional array, with v(:,i,j) the 3-vector to be matched at (x(i),y(j)), i=1:5, j=1:5. Note further that, in accordance with conds{1} being 'clamped', size(v,2) is 7 (and not 5), with the first and last entry of v(r,:,j) specifying the end slopes to be matched.

End conditions other than the ones listed earlier can be handled along the following lines. Suppose that we want to enforce the condition

for given scalars a, b, and c, and with e := x(1). Then one could compute the cubic spline interpolant s₁ to the given data using the default end condition as well as the cubic spline interpolant s₀ to zero data and some (nontrivial) end condition at e, and then obtain the desired interpolant in the form

Here are the (not inconsiderable) details (in which the first polynomial piece of s₁ and s₀ is pulled out to avoid differentiating all of s₁ and s₀):

pp1 = csape(x,y);
dp1 = fnder(fnbrk(pp1,1)); 
pp0 = csape(x,zeros(size(y)),[1,0],[1,0]); 
dp0 = fnder(fnbrk(pp0,1));
e = x(1); 
lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e); 
lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);
pp = fncmb(pp0,(c-lam1)/lam0,pp1);

Algorithm

The relevant tridiagonal linear system is constructed and solved using MATLAB's sparse matrix capabilities.

See Also

csapi, spapi, spline

Cautionary Note

If the sequence x is not nondecreasing, both x and y will be reordered in concert to make it so. Also, if the value sequence y is vector-valued, then valconds(:,j), j=1:2, must be vectors of that same length (if explicitly given).

chbpnt

csapi