Learn more about Spline Toolbox

spapi - Spline interpolation

Syntax

spline = spapi(knots,x,y) spapi(k,x,y) spapi({knork1,...,knorkm},{x1,...,xm},y) spapi(...,'noderiv')

Description

spline = spapi(knots,x,y) returns the spline (if any) of order

k = length(knots) - length(x)

with knot sequence knots for which

(*)   f(x(j)) = y(:,j), all j.

If some of the entries of x are the same, then this is taken in the osculatory sense, i.e., in the sense that , with m(j) := #{ i<j : x(i) = x(j) }, and the -th derivative of . Thus -fold repetition of a site in x corresponds to the prescribing of value and the first derivatives of at . If you don't want this, call spapi with an additional (fourth) argument, in which case, at each data site, the average of all data values with the same data site is matched.

The data values, y(:,j), may be scalars, vectors, matrices, or even ND-arrays.

spapi(k,x,y) , with k a positive integer, merely specifies the desired spline order, k, in which case aptknt is used to determine a workable (though not necessarily optimal) knot sequence for the given sites x. In other words, the command spapi(k,x,y) has the same effect as the more explicit command spapi(aptknt(x,k),x,y).

spapi({knork1,...,knorkm},{x1,...,xm},y) returns the B-form of a tensor-product spline interpolant to gridded data. Here, each knorki is either a knot sequence, or else is a positive integer specifying the polynomial order to be used in the ith variable, thus leaving it to spapi to provide a corresponding knot sequence for the ith variable. Further, y must be an (r+m)-dimensional array, with y(:,i1,...,im) the datum to be fitted at the site [x{1}(i1),...,x{m}(im)], all i1, ..., im , unless the spline is to be scalar-valued, in which case, in contrast to the univariate case, y is permitted to be an m-dimensional array.

spapi(...,'noderiv') with the string 'noderiv' as a fourth argument, has the same effect as spapi(...) except that data values sharing the same site are interpreted differently. With the fourth argument present, the average of the data values with the same data site is interpolated at such a site. Without it, data values with the same data site are interpreted as values of successive derivatives to be matched at such a site, as described above, in the first paragraph of this Description.

Examples

spapi([0 0 0 0 1 2 2 2 2],[0 1 1 1 2],[2 0 1 2 -1])produces the unique cubic spline f on the interval [0..2] with exactly one interior knot, at 1, that satisfies the five conditions

, , , ,

These include 3-fold matching at 1, i.e., matching there to prescribed values of the function and its first two derivatives.

Here is an example of osculatory interpolation, to values y and slopes s at the sites x by a quintic spline:

sp = spapi(augknt(x,6,2),[x,x,min(x),max(x)],[y,s,ddy0,ddy1]);

with ddy0 and ddy1 values for the second derivative at the endpoints.

As a related example, if the function sin(x) is to be interpolated at the distinct data sites x by a cubic spline, and its slope is also to be matched at a subsequence x(s), then this can be accomplished by the command

sp = spapi(4,[x x(s)], [sin(x) cos(x(s))]);

in which a suitable knot sequence is supplied with the aid of aptknt. In fact, if you wanted to interpolate the same data by quintic splines, simply change the 4 to 6.

As a bivariate example, here is a bivariate interpolant.

x = -2:.5.2; y = -1:.25:1; [xx, yy] = ndgrid(x,y);
z = exp(-(xx.^2+yy.^2));
sp = spapi({3,4},{x,y},z); fnplt(sp)

As an illustration of osculatory interpolation to gridded data, here is complete bicubic interpolation, with the data explicitly derived from the bicubic polynomial , to make it easy for you to see exactly where the slopes and slopes of slopes (i.e., cross derivatives) must be placed in the data values supplied. Since our is a bicubic polynomial, its interpolant, , must be itself. We test this.

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 = spapi({4,4}, {sites{1}([1,2,1,2]),sites{2}([1,2,1,2])}, ...
          [fnval(g,sites), Dyg ; ...
           Dxg.'         , Dxyg]);
if any( squeeze( fnbrk(fn2fm(f,'pp'), 'c') ) - coefs )
'something went wrong', end

Algorithm

spcol is called on to provide the almost-block-diagonal collocation matrix (with repeats in x denoting derivatives, as described above), and slvblk solves the linear system (*), using a block QR factorization.

Gridded data are fitted, in tensor-product fashion, one variable at a time, taking advantage of the fact that a univariate spline fit depends linearly on the values being fitted.

Limitations

The given (univariate) knots and sites must satisfy the Schoenberg-Whitney conditions for the interpolant to be defined. Assuming the site sequence x to be nondecreasing, this means that we must have

, all