Tridiagonal systems, determinants & natural cubic splines

Tridiagonal matrices

A tridiagonal matrix is a matrix that has nonzero entries only on the main diagonal and on the adjacent off-diagonals. This special structure comes up frequently in applications. For example, the finite difference numerical solution to the heat equation leads to a tridiagonal system. Another application, the one we’ll look at in detail here, is natural cubic splines. And we’ll mention an interesting result with Fibonacci numbers in passing.

Because these matrices have a special structure, the corresponding linear systems are quick and easy to solve. Also, we can find a simple recurrence relation for their determinants.

Determinants

Let’s label the component of a tridiagonal matrix as below

$M_n = \left( \begin{array}{lllll} a_1 & b_1 & & & \\ c_1 & a_2 & b_2 & & \\ & c_2 & \ddots & \ddots & \\ & & \ddots & \ddots & b_{n-1} \\ & & & c_{n-1} & a_n \\ \end{array} \right)$

where every component not shown is implicitly zero. We can expand the determinant of the matrix above using minors along the last row. This gives a recursive expression for the determinant

$d_n = a_n d_{n-1} - c_{n-1} b_{n-1} d_{n-2}$

with initial conditions d₀ = 1 and d_-1 = 0.

Note that if all the a‘s and b‘s are 1 and all the c‘s are -1, then you get the recurrence relation that defines the Fibonacci numbers. That is, the Fibonacci numbers are given by the determinant

$F_{n+1} = \left| \begin{array}{rrrrr} 1 & 1 & & & \ldots \\ -1 & 1 & 1 & & \ldots \\ & -1 & \ddots & \ddots & \\ & & \ddots & \ddots & 1 \\ & & & -1 & 1 \\ \end{array} \right|$

Natural cubic splines

A cubic spline interpolates a set of data points with piecewise cubic polynomials. There’s a (potentially) different cubic polynomial over each interval between input values, all fitted together so that the resulting function, its derivative, and its second derivative are all continuous.

Suppose you have input points, called knots in this context, t₀, t₁, … t_n and output values y₀, y₁, … y_n. For the spline to interpolate the data, its value at t_i must be y_i.

A cubic spline then is a set of n cubic polynomials, one for each interval [t_i, t_i+1]. A cubic polynomial has four coefficients, so we have 4n coefficients in total. At each interior knot, t₁ through t_n-1, we have four constraints. Both polynomials that meet at t_i must take on the value y_i at that point, and the two polynomials must have the same first and second derivative at that point. That gives us 4(n-1) equations. The value of the first polynomial is specified on the left end at t₀ the value of the last polynomial is specified at the right end at t_n. This gives us 4n – 2 equations in total.

We need two more equations. A clamped cubic spline specifies the derivatives at each end point. The natural cubic spline specifies instead that the second derivatives at each end are zero. What is natural about a natural cubic spline? In a certain sense it is the smoothest curve interpolating the specified points. With these boundary conditions we now have as many constraints as degrees of freedom.

So how would we go about finding the coefficients of each polynomial? Our task will be much easier if we parameterize the polynomials cleverly to start with. Instead of powers of x, we want powers of (x – t_i) and (t_i+1 – x) because these expressions are 0 on different ends of the interval [t_i, t_i+1]. It turns out we parameterize the spline over the ith interval as

$\begin{eqnarray*} S_i(x) &=& \frac{z_{i+1}}{6h_i}(x-t_i)^3 + \frac{z_i}{6h_i}(t_{i+1} - x)^3 \\ && + \left(\frac{y_{i+1}}{h_i} - \frac{h_iz_{i+1}}{6}\right) (x - t_i) \\ && + \left(\frac{y_{i }}{h_i} - \frac{h_iz_i }{6}\right) (t_{i+1} - x) \end{eqnarray*}$

where h_i = t_i+1 – t_i, the length of the ith interval.

This may seem unmotivated, and no doubt it is cleaner than the first thing someone probably tried, but it’s the kind of thing you’re lead to when you try to make the derivation work out smoothly.

The basic form is powers of (x – t_i) and (t_i+1 – x), each to the first and third powers, for reasons given above. Why the 6’s in the denominators? They’re not strictly necessary, but they cancel out when we take second derivatives. Let’s look at the second derivative.

$S_i''(x) = \frac{z_{i+1}(x - t_i)}{h_i} - \frac{z_{i}(t_{i+1} - x)}{h_i}$

Note how when we stick in t_i the first term is zero and the second is z_i, and when we stick in t_i+1 the first term is z_i+1 and the second is zero.

We can now write down the system of equations for the z‘s. We have z₀ = z_n from the natural cubic spline condition, and for 1 ≤ i ≤ n-1 we have

$h_{i-1}z_{i-1} + u_i z_i + h_i z_{i+1} = v_i$

where

$\begin{eqnarray*}b_i &=& (y_{i+1} - y_i)/h_i \\ u_i &=& 2(h_{i-1} + h_i) \\ v_i &=& 6(b_i - b_{i-1}) \end{eqnarray*}$
Note that this is a tridiagonal system because the ith equation only involves z‘s with subscripts i-1, i, and i+1.

Because of its tridiagonal structure, these equations can be solved simply and efficiently, much more efficiently than a general system of equations.

Tridiagonal systems, determinants, and natural cubic splines

Tridiagonal matrices

Determinants

Natural cubic splines

Related math posts

2 thoughts on “Tridiagonal systems, determinants, and natural cubic splines”