spcol

B-spline collocation matrix

Syntax

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

Description

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} (tau (i))$

This is the value at tau(i) of the m(i)th derivative of the jth 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.

Examples

To solve approximately the non-standard second-order ODE

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

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²y(t) – 5· (y(t) – sin(2t)), 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)

Limitations

The sequence tau is assumed to be nondecreasing.

Algorithms

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₀, j₀ – 1,...,0, with j₀ equal to max(m), values of all B-splines of that order are generated by recurrence and used to compute the jth 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.