Tridiagonal systems, determinants, and 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 d0 = 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, t0, t1, … tn and output values y0, y1, … yn. For the spline to interpolate the data, its value at ti must be yi.

A cubic spline then is a set of n cubic polynomials, one for each interval [ti, ti+1]. A cubic polynomial has four coefficients, so we have 4n coefficients in total. At each interior knot, t1 through tn-1, we have four constraints. Both polynomials that meet at ti must take on the value yi 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 t0 the value of the last polynomial is specified at the right end at tn. This gives us 4n – 2 equations in total.

We need two more equations. A clamped cubic spline species 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 (xti) and (ti+1 – x) because these expressions are 0 on different ends of the interval [ti, ti+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 hi = ti+1ti, 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 (xti) and (ti+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 ti the first term is zero and the second is zi, and when we stick in ti+1 the first term is zi+1 and the second is zero.

We can now write down the system of equations for the z‘s. We have z0 = zn from the natural cubic spline condition, and for 1 ≤ in-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.

Related posts

Leave a Reply

Your email address will not be published. Required fields are marked *