## 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); ```