Spline interpolation

`spline = spapi(knots,x,y) `

spapi(k,x,y)

spapi({knork1,...,knorkm},{x1,...,xm},y)

spapi(...,'noderiv')

`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 this is taken in the osculatory
sense, i.e., in the sense that *D ^{m}*

`m(j)`

:
= #{ i < j : `x(i)`

= `x(j)`

},
and `x`

corresponds
to the prescribing of value and the first `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 `i`

th variable, thus leaving it to `spapi`

to
provide 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 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 character vector `'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.

`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

*f*(0+) = 2, *f*(1)
= 0, *Df*(1) = 1, *D*^{2}*f*(1)
= 2, *f*(2–) = –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 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 *g*(*u*,*v*)
= *u*^{3}*v*^{3},
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 *g* is a bicubic polynomial,
its interpolant, *f*, must be *g* 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

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

$$\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.

`spcol`

is called on to provide the almost-block-diagonal collocation matrix
(*B _{j}*,

`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.