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 f (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:

$$D^{m(j)}f(x(j))=y(:,j)$$

with $$m(j)=\{i<j:x(i)=x(j)\}$$ and D^mf the m-th derivative of f. In this case, the r-fold repetition of a site z in x corresponds to the prescribing of value and the first r – 1 derivatives of f at z. To match the average of all data values with the same data instead, call spapi with an additional fourth argument.

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

spapi(k,x,y) , with k a positive integer, specifies the desired spline order, k. In this case the spapi function calls the aptknt function to determine a workable, but 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 a positive integer specifying the polynomial order used in the i-th variable. The spapi function then provides a corresponding knot sequence for the i-th variable. Further, y must be an (r+m)-dimensional array, with y(:,i1,...,im) the datum to fit at the site [x{1}(i1),...,x{m}(im)], for all i1, ..., im. In contrast to the univariate case, if the spline is scalar-valued, then y can be an m-dimensional array.

spapi(...,'noderiv') with the character vector or string scalar '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

Bivariate Spline Interpolant and Osculatory Interpolation to Gridded Data

The function 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

$$\mathit{f}\left(0+\right)=2,\mathit{f}\left(1\right)=0,\mathrm{Df}\left(1\right)=1,{\mathit{D}}^2\mathit{f}\left(1\right)=2,\mathit{f}\left(2-\right)=-1.$$

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 you want to interpolate the sin(x) function at the distinct data sites by a cubic spline, and to match its slope at a subsequence x(s), then call the spapi function with these arguments:

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

The aptknt function will provide a suitable knot sequence. If you want to interpolate the same data by quintic splines, then simply change the value 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 $\mathit{g}\left(\mathit{u},\mathit{v}\right)={\mathit{u}}^3{\mathit{v}}^3$. This is helpful to see exactly where the slopes, and slopes of slopes (the cross derivatives), must be placed in the data values supplied. Since g is a bicubic polynomial, its interpolant, f, must be g itself. 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

Input Arguments

knots — Sequence of knots

vector | cell array

Knot sequence of the spline, specified as a nondecreasing vector.

k — Evaluation points

vector | matrix | cell

Points at which you want to evaluate the spline function f, specified as a vector, matrix or cell array.

x — Data sites

vector | cell array

Data sites of data values y to be fit, specified as a vector or as a cell array for multivariate data. Spline f is created with knots at each data site x such that f(x(j)) = y(:,j) for all values of j.

For multivariate, gridded data, you can specify x as a cell array that specifies the data site in each variable dimension: f(x1(i),x2(j),...xn(k)) = y(:,i,j,...,k).

y — Data values to fit

vector | matrix | array

Data values to fit during creation of the spline, specified as a vector, matrix, or array. Data values y(:,j) can be scalars, matrices, or n-dimensional arrays. Data values given at the same data site x are averaged.

Data Types: single | double

Output Arguments

spline — Spline structure

spline structure

Spline, returned as a structure with these fields.

Form — Form of spline

pp | B-

Form of the spline, returned as pp or B-. pp indicates that the spline is given in piecewise polynomial form, B- indicates it is given in B-form.

Knots — Knot locations of spline

vector | cell array

Knot positions of the spline, returned as a vector or as a cell array of vectors for multivariate data. Vectors contain strictly increasing elements that represent the start and end of each of the intervals over which the polynomial pieces are defined.

Coefs — Coefficients of polynomials

matrix | array

Coefficients of polynomials for each piece, returned as a matrix or as an array for multivariate data.

Number — Number of polynomial pieces

scalar | vector

Number of polynomial pieces describing the spline, returned as a scalar or as a vector of numbers of pieces in each variable for multivariate data.

Order — Order of polynomials

scalar | vector

Order of the polynomial function describing each polynomial piece of the spline, returned as a scalar or as a vector containing the order in each variable for multivariate data.

Dim — Dimensionality

scalar

Dimensionality of the target function, returned as a scalar.

Limitations

The given (univariate) knots and sites must satisfy the Schoenberg-Whitney conditions for the interpolant to be defined. If the site sequence x is nondecreasing, then

$$\text{knots}(j)<x(j)<\text{knots}(j+k),\text{ all }j$$

with equality possible at knots(1) and knots(end)). In the multivariate case, these conditions must hold in each variable separately.

Algorithms

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

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

Version History

Introduced in R2006b