Cubic spline interpolation with end conditions
pp = csape(x,y)
pp = csape(x,y,conds)
pp = csape(x,y)
is the
ppform of a cubic spline s with knot sequence x
that
satisfies s(x(j)) = y(:,j)
for
all j
, as well as an additional end condition at
the ends (meaning the leftmost and at the rightmost data site), namely
the default condition listed below. The data values y(:,j)
may
be scalars, vectors, matrices, even NDarrays. Data values at the
same data site are averaged.
pp = csape(x,y,conds)
lets you choose the end conditions to be used, from a rather large
and varied catalog, by proper choice of conds
.
If needed, you supply the corresponding end condition values as additional
data values, with the first (last) data value taken as the end condition
value at the left (right) end. In other words, in that case, s(x(j))
matches y(:,j+1)
for
all j
, and the variable endcondvals used in the
detailed description below is set to y(:,[1 end])
.
For some choices of conds
, these end condition
values need not be present and/or are ignored when present.
conds
may be a character vector whose
first character matches one of the following: 'complete'
or '
clamped'
, 'notaknot'
, 'periodic'
, 'second'
, 'variational'
,
with the following meanings.
 Match endslopes (as given, with default as under “default”). 
Make second and secondlast sites inactive knots (ignoring end condition values if given).  
 Match first and second derivatives at left end with those at right end. 
 Match end second derivatives (as given, with default
[0 0], i.e., as in 
 Set end second derivatives equal to zero (ignoring end condition values 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 1by2 matrix instead,
it is possible to specify different conditions
at the two ends. Explicitly, the ith derivative, D^{i}s,
is given the value endcondvals(:,j)
at the left (j is 1) respectively right (j is 2) end in case conds(j)
is i,i = 1:2. There are default
values for conds
and/or endcondvals.
clamped  Ds(e) = endcondvals(:, 

curved  D^{2}s(e)
= endcondvals(:, 

Lagrange  Ds(e) = Dp(e)  default 
periodic  D^{r}s(a) = D^{r}s(b), r = 1,2 

variational  D^{2}s(e) = 0 

Here, e is a (e is b),
i.e., the left (right) end, in case j
is 1 (j
is
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 endcondvals(:,j)
is
taken to be the corresponding default value.
The default value for endcondvals(:,j
) is
the derivative of the cubic interpolant at the nearest four sites
in case conds(j)
is 1, and is 0 otherwise.
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 mdimensional
array (or an m+rdimensional
array if the function is to be rvalued). Correspondingly, conds
is
a cell array with m entries, and end condition
values may be correspondingly supplied in each of the m variables.
This, as the last example below, of bicubic spline interpolation,
makes clear, may require you to supply end conditions for end conditions.
This command calls on a much expanded version of the Fortran
routine CUBSPL
in PGS.
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],[3 1 1 6],[1 2])
provides
the cubic polynomial p for which Dp(–1)
= 3, p(–1) = –1, p(1)
= 1, D^{2}p(1)
= 6, i.e., p(x) = x^{3}.
Finally, csape([1 1],[1 1])
provides the straight
line p for which p(±1) = ±1, i.e., p(x)
= x.
End conditions other than the ones listed earlier can be handled along the following lines. Suppose that you want to enforce the condition
$$\lambda (s):=aDs(e)+b{D}^{2}s(e)=c$$
for given scalars a, b,
and c, and with e equal to x
(1).
Then one could compute the cubic spline interpolant s_{1} to
the given data using the default end condition as well as the cubic
spline interpolant s_{0} to
zero data and some (nontrivial) end condition at e,
and then obtain the desired interpolant in the form
$$s={s}_{1}+\left((c\lambda )({s}_{1})\right)/\lambda ({s}_{0}){s}_{0}$$
Here are the (not inconsiderable) details (in which the first polynomial piece of s_{1} and s_{0} is pulled out to avoid differentiating all of s_{1} and s_{0}):
% Data: x and y [x, y] = titanium(); % Scalars a, b, and c a = 2; b = 1; c = 0; % End condition at left e = x(1); % The cubic spline interpolant s1 to the % given data using the default end % condition s1 = csape(x,y); % The cubic spline interpolant s0 to % zero data and some (nontrivial) end % condition at e s0 = csape(x,[1,zeros(1,length(y)),0],[1,0]); % Compute the derivatives of the first % polynomial piece of s1 and s0 ds1 = fnder(fnbrk(s1,1)); ds0 = fnder(fnbrk(s0,1)); % Compute interpolant with desired end conditions lam1 = a*fnval(ds1,e) + b*fnval(fnder(ds1),e); lam0 = a*fnval(ds0,e) + b*fnval(fnder(ds0),e); pp = fncmb(s0,(clam1)/lam0,s1);
Plot to see the results:
fnplt( pp, [594, 632] ) hold on fnplt( s1, 'b', [594, 632] ) plot( x, y, 'ro', 'MarkerFaceColor', 'r' ) hold off axis( [594, 632, 0.62, 0.655] ) legend 'Desired endconditions' ... 'Default endconditions' 'Data' ... Location SouthEast
As a multivariate vectorvalued example, here is a sphere, done as a parametric bicubic spline, 3Dvalued, using prescribed slopes in one direction and periodic end conditions in the other:
x = 0:4; y=2:2; s2 = 1/sqrt(2); v = zeros( 3, 7, 5 ); v(1,:,:) = [1 0 s2 1 s2 0 1].'*[1 0 1 0 1]; v(2,:,:) = [1 0 s2 1 s2 0 1].'*[0 1 0 1 0]; v(3,:,:) = [0 1 s2 0 s2 1 0].'*[1 1 1 1 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 3dimensional array, with v(:,i+1,j)
the
3vector 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.
Here is a bivariate example that shows the need for supplying end conditions of end conditions when supplying end conditions in both variables. You reproduce the bicubic polynomial g(x,y) = x^3y^3 by complete bicubic interpolation. You then derive the needed data, including end condition values, directly from g in order to make it easier for you to see just how the end condition values must be placed. Finally, you check the result.
sites = {[0 1],[0 2]}; coefs = zeros(4, 4); coefs(1,1) = 1; g = ppmak(sites,coefs); Dxg = fnval(fnder(g,[1 0]),sites); Dyg = fnval(fnder(g,[0 1]),sites); Dxyg = fnval(fnder(g,[1 1]),sites); f = csape(sites,[Dxyg(1,1), Dxg(1,:), Dxyg(1,2); ... Dyg(:,1), fnval(g,sites), Dyg(:,2) ; ... Dxyg(2,1), Dxg(2,:), Dxyg(2,2)], ... {'complete','complete'}); if any(squeeze(fnbrk(f,'c'))coefs) disp( 'this is wrong' ) end
csape
recognizes that you supplied explicit
end condition values by the fact that you supplied exactly two more
data values than data sites. In particular, even when using different
end conditions at the two ends, if you wish to supply an end condition
value at one end, you must also supply one for the other end.
The relevant tridiagonal linear system is constructed and solved using the sparse matrix capabilities of MATLAB^{®}.