Main Content

B-spline collocation matrix

`colmat = spcol(knots,k,tau) `

colmat = spcol(knots,k,tau,arg1,arg2,...)

`colmat = spcol(knots,k,tau) `

returns the
matrix, with `length(tau)`

rows and
`length(knots)-k`

columns, whose
(*i*,*j*)th entry is

$${D}^{m}{}^{(i)}{B}_{j}(\text{tau}(i))$$

This is the value at tau(*i*) of the
*m*(*i*)th derivative of the
*j*th B-spline of order `k`

for the knot sequence `knots`

. Here, `tau`

is
a sequence of sites, assumed to be *nondecreasing*, and
*m* = knt2mlt(tau), i.e.,
*m*(*i*) is
#{*j* < *i*:tau(*j*) = tau(*i*)},
all *i*.

`colmat = spcol(knots,k,tau,arg1,arg2,...) `

also returns that matrix, but gives you the opportunity to specify some
aspects.

If one of the `argi`

is a character vector or string scalar with
the same first two letters as in `'slvblk'`

, the matrix is returned
in the almost block-diagonal format (specialized for splines) required by `slvblk`

(and understood by
`bkbrk`

).

If one of the `argi`

is a character vector or string scalar with
the same first two letters as in `'sparse'`

, then the matrix is
returned in the `sparse`

format of MATLAB^{®}.

If one of the `argi`

is a character vector or string scalar with
the same first two letters as in `'noderiv'`

, multiplicities are
ignored, i.e., *m*(*i*) is taken to be 1 for all
*i*.

To solve approximately the non-standard second-order ODE

$${D}^{2}y(t)=5\cdot (y(t)-\mathrm{sin}(2t))$$

on the interval [0..π], using cubic splines with 10 polynomial pieces, you can use
`spcol`

in the following way:

tau = linspace(0,pi,101); k = 4; knots = augknt(linspace(0,pi,11),k); colmat = spcol(knots,k,brk2knt(tau,3)); coefs = (colmat(3:3:end,:)/5-colmat(1:3:end,:))\(-sin(2*tau).'); sp = spmak(knots,coefs.');

You can check how well this spline satisfies the ODE by computing and plotting the
residual,
*D*^{2}*y*(*t*)
– 5· (*y*(*t*) – sin(2*t*)),
on a fine mesh:

t = linspace(0,pi,501); yt = fnval(sp,t); D2yt = fnval(fnder(sp,2),t); plot(t,D2yt - 5*(yt-sin(2*t))) title(['residual error; to be compared to max(abs(D^2y)) = ',... num2str(max(abs(D2yt)))])

The statement `spcol([1:6],3,.1+[2:4])`

provides the
matrix

ans = 0.5900 0.0050 0 0.4050 0.5900 0.0050 0 0.4050 0.5900

in which the typical row records the values at 2.1, or 3.1, or 4.1, of all
B-splines of order 3 for the knot sequence `1:6`

. There are three
such B-splines. The first one has knots 1,2,3,4, and its values are recorded in the
first column. In particular, the last entry in the first column is zero since it
gives the value of that B-spline at 4.1, a site to the right of its last
knot.

If you add the character vector or string scalar `'sl'`

as an
additional input to `spcol`

, then you can ask `bkbrk`

to extract detailed
information about the block structure of the matrix encoded in the resulting output
from `spcol`

. Thus, the statement
`bkbrk(spcol(1:6,3,.1+2:4,'sl'))`

gives:

block 1 has 2 row(s) 0.5900 0.0050 0 0.4050 0.5900 0.0050 next block is shifted over 1 column(s) block 2 has 1 row(s) 0.4050 0.5900 0.0050 next block is shifted over 2 column(s)

The sequence `tau`

is assumed to be nondecreasing.

This is the most complex command in this toolbox since it has to deal with various
ordering and blocking issues. The recurrence relations are used to generate, simultaneously, the values of
all B-splines of order `k`

having anyone of the
`tau(i)`

in their support.

A separate calculation is carried out for the (presumably few) sites at which
derivative values are required. These are the sites `tau(i)`

with
*m*(*i*) > 0. For these, and for every
order *k* – *j*, *j* =
*j*_{0},
*j*_{0} – 1,...,0, with
*j*_{0} equal to max(*m*),
values of all B-splines of that order are generated by recurrence and used to
compute the *j*th derivative at those sites of all B-splines of
order `k`

.

The resulting rows of B-spline values (each row corresponding to a particular
`tau(i)`

) are then assembled into the overall (usually rather
sparse) matrix.

When the optional argument `'sl'`

is present, these rows are
instead assembled into a convenient almost block-diagonal form that takes advantage
of the fact that, at any site `tau(i)`

, at most
`k`

B-splines of order `k`

are nonzero. This
fact (together with the natural ordering of the B-splines) implies that the
collocation matrix is almost block-diagonal, i.e., has a staircase shape, with the individual
blocks or steps of varying height but of uniform width `k`

.

The command `slvblk`

is designed to take advantage of this
storage-saving form available when used, in `spap2`

,
`spapi`

, or `spaps`

, to help determine the
B-form of a piecewise-polynomial function from interpolation or other approximation
conditions.