## Cubic Spline Interpolation

### Cubic Spline Interpolant of Smooth Data

Suppose you want to interpolate some smooth data, e.g., to

rng(6), x = (4*pi)*[0 1 rand(1,15)]; y = sin(x);

You can use the cubic spline interpolant obtained by

cs = csapi(x,y);

and plot the spline, along with the data, with the following code:

fnplt(cs); hold on plot(x,y,'o') legend('cubic spline','data') hold off

This produces a figure like the following.

**Cubic Spline Interpolant of Smooth Data**

This is, more precisely, the cubic spline interpolant with the
not-a-knot end conditions, meaning that it is the unique piecewise
cubic polynomial with two continuous derivatives with breaks at all *interior* data
sites except for the leftmost and the rightmost one. It is the same
interpolant as produced by the MATLAB^{®} `spline`

command, `spline(x,y)`

.

### Periodic Data

The sine function is 2π-periodic. To check how well your interpolant does on that score, compute, e.g., the difference in the value of its first derivative at the two endpoints,

diff(fnval(fnder(cs),[0 4*pi])) ans = -.0100

which is not so good. If you prefer to get an interpolant whose first and second
derivatives at the two endpoints, `0`

and `4*pi`

, match, use
instead the command `csape`

which permits specification of many different kinds
of end conditions, including periodic end conditions. So, use instead

pcs = csape(x,y,'periodic');

for which you get

diff(fnval(fnder(pcs),[0 4*pi]))

Output is `ans = 0`

as the difference of end slopes. Even the difference
in end second derivatives is small:

diff(fnval(fnder(pcs,2),[0 4*pi]))

Output is `ans = -4.6074e-015`

.

### Other End Conditions

Other end conditions can be handled as well. For example,

cs = csape(x,[3,y,-4],[1 2]);

provides the cubic spline interpolant with breaks at the and with its slope at the leftmost data site equal to 3, and its second derivative at the rightmost data site equal to -4.

### General Spline Interpolation

If you want to interpolate at sites other than the breaks and/or
by splines other than cubic splines with simple knots, then you use the `spapi`

command.
In its simplest form, you would say `sp = spapi(k,x,y)`

;
in which the first argument, `k`

, specifies the *order *of
the interpolating spline; this is the number of coefficients in each
polynomial piece, i.e., 1 more than the nominal degree of its polynomial
pieces. For example, the next figure shows a linear, a quadratic,
and a quartic spline interpolant to your data, as obtained by the
statements

sp2 = spapi(2,x,y); fnplt(sp2,2), hold on sp3 = spapi(3,x,y); fnplt(sp3,2,'k--'), sp5 = spapi(5,x,y); fnplt(sp5,2,'r-.'), plot(x,y,'o') legend('linear','quadratic','quartic','data'), hold off

**Spline Interpolants of Various Orders of Smooth
Data**

Even the cubic spline interpolant obtained from `spapi`

is different from the one provided
by `csapi`

and `spline`

. To emphasize their difference,
compute and plot their second derivatives, as follows:

fnplt(fnder(spapi(4,x,y),2)), hold on, fnplt(fnder(csapi(x,y),2),2,'k--'),plot(x,zeros(size(x)),'o') legend('from spapi','from csapi','data sites'), hold off

This gives the following graph:

**Second Derivative of Two Cubic Spline Interpolants
of the Same Smooth Data**

Since the second derivative of a cubic spline is a broken line,
with vertices at the breaks of the spline, you can see clearly that `csapi`

places
breaks at the data sites, while `spapi`

does not.

### Knot Choices

It is, in fact, possible to specify explicitly just where the
spline interpolant should have its breaks, using the command ```
sp
= spapi(knots,x,y)
```

; in which the sequence `knots`

supplies,
in a certain way, the breaks to be used. For example, recalling that
you had chosen `y`

to be `sin(x)`

,
the command

ch = spapi(augknt(x,4,2),[x x],[y cos(x)]);

provides a cubic Hermite interpolant to the sine function, namely
the piecewise cubic function, with breaks at all the `x(i)`

's,
that matches the sine function in value *and* slope
at all the `x(i)`

's. This makes the interpolant continuous
with continuous first derivative but, in general, it has jumps across
the breaks in its second derivative. Just how does this command know
which part of the data value array `[y cos(x)]`

supplies
the values and which the slopes? Notice that the data site array here
is given as `[x x]`

, i.e., each data site appears
twice. Also notice that `y(i)`

is associated with
the first occurrence of `x(i)`

, and `cos(x(i))`

is
associated with the second occurrence of `x(i)`

.
The data value associated with the first appearance of a data site
is taken to be a function value; the data value associated with the
second appearance is taken to be a slope. If there were a third appearance
of that data site, the corresponding data value would be taken as
the second derivative value to be matched at that site. See Constructing and Working with B-form Splines for
a discussion of the command `augknt`

used here to
generate the appropriate "knot sequence".

### Smoothing

What if the data are noisy? For example, suppose that the given values are

noisy = y + .3*(rand(size(x))-.5);

Then you might prefer to approximate instead. For example, you might try the cubic smoothing spline, obtained by the command

scs = csaps(x,noisy);

and plotted by

fnplt(scs,2), hold on, plot(x,noisy,'o'), legend('smoothing spline','noisy data'), hold off

This produces a figure like this:

**Cubic Smoothing Spline of Noisy Data**

If you don't like the level of smoothing done by `csaps(x,y)`

,
you can change it by specifying the smoothing parameter, `p`

,
as an optional third argument. Choose this number anywhere between
0 and 1. As `p`

changes from 0 to 1, the smoothing
spline changes, correspondingly, from one extreme, the least squares
straight-line approximation to the data, to the other extreme, the
"natural" cubic spline interpolant to the data. Since `csaps`

returns
the smoothing parameter actually used as an optional second output,
you could now experiment, as follows:

[scs,p] = csaps(x,noisy); fnplt(scs,2), hold on fnplt(csaps(x,noisy,p/2),2,'k--'), fnplt(csaps(x,noisy,(1+p)/2),2,'r:'), plot(x,noisy,'o') legend('smoothing spline','more smoothed','less smoothed',... 'noisy data'), hold off

This produces the following picture.

**Noisy Data More or Less Smoothed**

At times, you might prefer simply to get the smoothest cubic
spline `sp`

that is within a specified tolerance `tol`

of
the given data in the sense that ```
norm(noisy - fnval(sp,x))^2
<= tol
```

. You create this spline with the command ```
sp
= spaps(x,noisy,tol)
```

for your defined tolerance `tol`

.

### Least Squares

If you prefer a least squares approximant, you can obtain it
by the statement `sp = spap2(knots,k,x,y)`

; in which
both the knot sequence `knots`

and the order `k`

of
the spline must be provided.

The popular choice for the order is 4, and that gives you a cubic spline. If you have no clear idea of how to choose the knots, simply specify the number of polynomial pieces you want used. For example,

sp = spap2(3,4,x,y);

gives a cubic spline consisting of three polynomial pieces.
If the resulting error is uneven, you might try for a better knot
distribution by using `newknt`

as
follows:

sp = spap2(newknt(sp),4,x,y);