1.6.2 Cubic splines | SigmaField

Corners

Linear interpolation leaves a corner at every inner point, where the slope jumps. In the chain of lesson 1.6.1 it jumps from 2 to $-1$ at $x = 1$ . Many quantities change smoothly, so a corner is usually a fault of the method, not a feature of the data.

The slope of a curve at a point is how steeply it climbs there. Measure it with a very short chord, the straight line joining two nearby points of the curve: its rise over run, as in lesson 1.2.1. Lesson 2.2.1 makes this exact.

A corner is a point where the slope just to the left and the slope just to the right disagree.

One polynomial through every point

A tempting fix is one polynomial through all the points, which has no corners anywhere. Through $n$ points with different $x$ values there is exactly one polynomial of degree at most $n - 1$ .

For 2 points it is a line, and for 3 a parabola. np.polyfit(xs, ys, n - 1) finds it.

It usually wiggles. Take the S-shaped curve $y = \frac{4}{1 + e^{-1.5(x - 6)}}$ , which climbs from 0 to 4 around $x = 6$ . Sample it at the 13 whole numbers from 0 to 12, and fit the polynomial of degree 12 through them.

Every data point lies between 0 and 4, yet the polynomial dips to $-0.514$ near $x = 0.3$ and rises to 4.514 near $x = 11.7$ . To pass through every point it must turn many times, and the turns grow near the ends. More points make it worse, not better:

points	7	9	13	17
polynomial, worst error	0.297	0.329	0.515	0.918
cubic spline, worst error	0.21	0.09	0.0181	0.0044

The worst error is the largest gap between the estimate and the true curve. The cubic spline, built in the next section, does the opposite: each extra point brings it closer.

Three ways through the same 13 points. The chain of lines has corners, the cubic spline is smooth, and the one polynomial swings outside the range of the data near both ends.

Try it The curve is sampled with 7, 9, 13 and 17 points between 0 and 12. Each line prints the number of points, then the worst error of the polynomial and of the spline, checked at 12,001 inputs. The last line is the 13-point polynomial's lowest and highest values, and where they occur.

import numpy as np from scipy.interpolate import CubicSpline def curve(x): return 4 / (1 + np.exp(-1.5 * (x - 6))) x = np.linspace(0, 12, 12001) for n in [7, 9, 13, 17]: xs = np.linspace(0, 12, n) poly = np.polyval(np.polyfit(xs, curve(xs), n - 1), x) spline = CubicSpline(xs, curve(xs), bc_type="natural")(x) print(n, round(np.abs(poly - curve(x)).max(), 3), round(np.abs(spline - curve(x)).max(), 4)) xs = np.arange(13.0) poly = np.polyval(np.polyfit(xs, curve(xs), 12), x) print(round(poly.min(), 3), round(poly.max(), 3), round(x[poly.argmin()], 1), round(x[poly.argmax()], 1))

⌘/Ctrl + Enter runs the cell

Cubic splines

A cubic spline keeps the low degree of linear interpolation but drops the corners. It joins each pair of neighbouring points with its own cubic, a polynomial of degree 3 (lesson 1.3.1). The cubics are chosen so that they meet smoothly.

Rule: Cubic spline

\begin{gathered} \text{one cubic for each gap,} \\ \text{through both of its points,} \\ \text{with no jump in slope} \\ \text{or in bend at a join} \end{gathered}

With $n$ points there are $n - 1$ gaps, so $n - 1$ cubics. SciPy's CubicSpline(xs, ys) finds them, and the result is called like a function.

A worked example: the tent. The points $(0, 0)$ , $(1, 1)$ and $(2, 0)$ make a tent. Linear interpolation gives two straight pieces with a sharp corner at the top. The natural spline (defined below) uses two cubics: on the first gap it is $s(x) = 1.5x - 0.5x^3$ , and the second is its mirror image.

By algebra, $s(0) = 0$ and $s(1) = 1.5 - 0.5 = 1$ , so the piece passes through both of its points. At the middle of the gap, $s(0.5) = 0.75 - 0.0625 = 0.6875$ . The table compares the spline with the straight piece.

$x$	0	0.25	0.5	0.75	1
straight piece	0	0.25	0.5	0.75	1
spline $s(x)$	0	0.3672	0.6875	0.9141	1

The spline bulges above the straight piece and rounds off the top of the tent.

-0.5x^3 + 0x^2 + 1.5x + 0

import numpy as np from scipy.interpolate import CubicSpline xs, ys = [0, 1, 2], [0, 1, 0] spline = CubicSpline(xs, ys, bc_type="natural") print(spline.c[:, 0]) for x in [0, 0.25, 0.5, 0.75, 1]: print(x, np.interp(x, xs, ys), round(float(spline(x)), 4))

⌘/Ctrl + Enter runs the cell

What "smooth" means at a join

At $x = 1$ the two straight pieces have slopes 1 and $-1$ , so the slope jumps by $-2$ : that is the corner. For the spline, take short chords on each side of 1. With a run of 0.01 their slopes are 0.01495 on the left and $-0.01495$ on the right.

With a run of 0.001 they are 0.0015 and $-0.0015$ . Both sides settle on the same slope, 0, so the spline crosses the top of the tent level, with no corner.

The tent two ways. The straight pieces meet at a corner, where the slope jumps from 1 to $-1$ . The spline is level at the top from both sides.

Bend is how fast the slope itself changes. At the top of the tent the spline bends downwards, like a cap; a curve that bends upwards looks like a cup. A spline also matches the bend at each join, so the curve never switches suddenly from turning gently to turning sharply.

So "smooth" at a join means no jump in slope and no jump in bend. SciPy reports both: spline(x, 1) is the slope at $x$ and spline(x, 2) the bend. Lessons 2.2.2 and 2.4.2 name them the derivative and the second derivative.

a

import numpy as np from scipy.interpolate import CubicSpline xs, ys = [0, 1, 2], [0, 1, 0] spline = CubicSpline(xs, ys, bc_type="natural") def chords(f, a, h): left = (f(a) - f(a - h)) / h right = (f(a + h) - f(a)) / h return round(float(left), 5), round(float(right), 5) print(chords(lambda x: np.interp(x, xs, ys), 1, 0.01)) print(chords(spline, 1, 0.01), chords(spline, 1, 0.001)) print(float(spline(1, 1)), round(float(spline(0.999999, 2)), 4), round(float(spline(1.000001, 2)), 4)) print(float(spline(0, 2)), float(spline(2, 2)))

⌘/Ctrl + Enter runs the cell

End conditions

Count what a spline must satisfy: with $n$ points there are $n - 1$ cubics of 4 coefficients each. Passing through its two points fixes 2 of each cubic's 4. Matching slope and bend fixes 2 more at each of the $n - 2$ inner joins.

points $n$	3	4	5	10
coefficients, $4(n - 1)$	8	12	16	36
through the points, $2(n - 1)$	4	6	8	18
smooth joins, $2(n - 2)$	2	4	6	16
still free	2	2	2	2

Two conditions are always missing, one for each end, where a piece has no neighbour to match. This choice is the end condition. It changes the curve, most of all near the ends.

Natural: the bend is 0 at both ends, so the curve runs straight as it leaves the data. Write bc_type="natural".
Not-a-knot, SciPy's default: the first two pieces are one and the same cubic, and so are the last two.
Clamped: the slope at each end is fixed; bc_type="clamped" sets it to 0.

Rule: Natural spline

\text{bend} = 0 \text{ at both ends}

The two missing conditions are spent on making the curve straight where the data stop. It suits a curve that really does level off or straighten out at its ends.

For the tent, the natural spline gives 0.6875 at 0.5 and the default gives 0.75. With only 3 points the default's two pieces are one parabola, $2x - x^2$ . Both pass through all three points.

The right choice depends on the curve. Sample $y = \frac{10}{1 + x}$ at 0, 1, …, 5. It bends sharply near 0, so forcing a straight end there costs accuracy: the natural spline's worst error is 0.5851, the default's 0.2785.

\frac{10}{1 + x}

import numpy as np from scipy.interpolate import CubicSpline natural = CubicSpline([0, 1, 2], [0, 1, 0], bc_type="natural") default = CubicSpline([0, 1, 2], [0, 1, 0]) print(float(natural(0.5)), float(default(0.5))) def f(x): return 10 / (1 + x) xs = np.arange(6.0) x = np.linspace(0, 5, 501) for bc in ["natural", "not-a-knot"]: spline = CubicSpline(xs, f(xs), bc_type=bc) print(bc, round(np.abs(spline(x) - f(x)).max(), 4))

⌘/Ctrl + Enter runs the cell

Splines against straight lines

On smooth data a spline is usually far more accurate. For $\frac{10}{1 + x}$ at 0 to 5, linear interpolation's worst error is 0.8578 and the default spline's 0.2785, about three times smaller.

A spline has three costs that straight lines do not.

It depends on every point. Raising the point at 4 by 0.1 moves the spline by $-0.01$ at 2.5, a gap away. The straight lines there do not move; case study 1.6.4 meets this with a bad quote.
It runs away outside the data. Past the last point CubicSpline keeps using the last cubic. It gives 0.6556 at 7, where the true value is 1.25, and $-0.4444$ at 8, although $\frac{10}{1 + x}$ is never negative. Pass extrapolate=False to get nan instead.
It can overshoot, as the next section shows.

Try it The first line is the worst error of np.interp, then of the spline. Next come the spline's values at 7 and 8 beside the true ones, and the change at 2.5 after raising one point, spline then lines. The last line refuses to extrapolate.

import numpy as np from scipy.interpolate import CubicSpline def f(x): return 10 / (1 + x) xs = np.arange(6.0) ys = f(xs) spline = CubicSpline(xs, ys) x = np.linspace(0, 5, 501) print(round(np.abs(np.interp(x, xs, ys) - f(x)).max(), 4), round(np.abs(spline(x) - f(x)).max(), 4)) print(np.round(spline([7, 8]), 4), np.round(f(np.array([7, 8])), 4)) raised = ys.copy() raised[4] += 0.1 print(round(float(CubicSpline(xs, raised)(2.5) - spline(2.5)), 4), np.interp(2.5, xs, raised) - np.interp(2.5, xs, ys)) print(CubicSpline(xs, ys, extrapolate=False)(7))

⌘/Ctrl + Enter runs the cell

Overshoot

To avoid corners, a spline may swing above or below the data. Take the points $(0, 0)$ , $(1, 0)$ , $(2, 1)$ and $(3, 1)$ , a step up. The natural spline through them dips to about $-0.128$ and rises to about 1.128, outside the range of the data.

That matters when values must stay within limits, such as a quantity that cannot be negative. Linear interpolation is then safer: between two points it stays between their values, as lesson 1.6.1 showed. SciPy also has splines built to keep the shape of the data, such as PchipInterpolator.

Try it The plot compares both methods on the step. The spline is smooth but overshoots at both ends of the rise; the straight pieces never leave the range of the data. The printed numbers are the spline's lowest and highest values over 301 inputs.

import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import CubicSpline xs = np.array([0, 1, 2, 3.0]) ys = np.array([0, 0, 1, 1.0]) spline = CubicSpline(xs, ys, bc_type="natural") x = np.linspace(0, 3, 301) print(round(spline(x).min(), 3), round(spline(x).max(), 3)) plt.plot(x, np.interp(x, xs, ys), label="linear") plt.plot(x, spline(x), label="cubic spline") plt.plot(xs, ys, "o", color="black") plt.legend() plt.show()

⌘/Ctrl + Enter runs the cell

Check yourself

Interpolated values must never fall below the smallest value in the data. Which method guarantees it?

A natural cubic splineLinear interpolationBoth

Common mistakes

A higher degree is not a better fit. One polynomial through many points wiggles, and more points make it worse.
No gaps is not the same as smooth. The chain of lines has no gaps, but its slope jumps at every corner. Smooth means no jump in slope or in bend.
The end condition changes the answer, most of all near the ends: 0.6875 against 0.75 for the tent.
A spline can overshoot. Use linear interpolation when values must stay inside the data's range.
Never trust a spline outside the data. Its last cubic runs away: $-0.44$ at 8 for a curve that is never negative.
One point moves the whole spline. A bad data point spreads into neighbouring gaps; with straight lines it stays in its own two.

Cubic splines

0 of 11 correct

1.
A cubic spline through 8 points has how many cubic pieces?
2.
What is the highest degree of the one polynomial through 10 points with different $x$ values?
3.
How many inner joins does a spline through 6 points have?
4.
A spline through 6 points uses 5 cubics. How many coefficients is that in all?
5.
A spline through 5 points has 16 coefficients. Passing through the points and smooth joins fix how many of them?
6.
The first piece of the tent spline is $1.5x - 0.5x^3$ . Evaluate it at $x = 0.6$
7.
Find the slope of the chord of $1.5x - 0.5x^3$ from 0.9 to 1
8.
The chain of lines through $(0, 1)$ , $(1, 2)$ and $(2, 0)$ has slopes 1 and $-2$ . By how much does the slope change at $x = 1$ ?
9.
What is the bend of a natural spline at its two end points?
10.
With 3 points, SciPy's default spline is one parabola. Through $(0, 0)$ , $(1, 1)$ and $(2, 0)$ it is $2x - x^2$ . What does it give at $x = 0.5$ ?
11.
Data values lie between 2 and 5. What is the largest value linear interpolation can give between the points?

How to type answers

Powers with ^: x^7, 2^(1/3), a^-2.
Multiply with * or by writing side by side: 3x^2, 2(x+1).
ln(x) is the natural log, log(x, 2) is log base 2, e is Euler's number, sqrt(x) is a square root.
sin(x), cos(x) and tan(x) take radians; write pi for π and sin(x)^2 for sin²x.
Your first answer to each question is saved: it counts towards your score and brings the question back for review.
Exact answers like ln(20)/ln(3) are best. Decimals need about 4 significant figures.

Summary card One polynomial through many points wiggles near the ends, and more points make it worse. A cubic spline uses one cubic per gap, through both of its points, with no jump in slope or bend at any join. The slope is how steeply the curve climbs; the bend is how fast the slope changes. Two end conditions complete a spline. Natural: bend 0 at both ends. SciPy's default, not-a-knot: the first two pieces are one cubic, and so are the last two. CubicSpline(xs, ys) builds it; call the result like a function, with bc_type="natural" for a natural spline. Splines are usually far more accurate on smooth data. But they can overshoot, every point affects every piece, and they run away outside the data.

Corners

One polynomial through every point

Cubic splines

Rule: Cubic spline

What "smooth" means at a join

End conditions

Rule: Natural spline

Splines against straight lines

Overshoot

Common mistakes

Summary card

Free further reading