\documentclass[12pt]{article}
\usepackage[]{amssymb,amsmath}

\title{Tridiagonal Systems}

\author{Charles Kankelborg}

\date{2001-01-21}

\begin{document}
% Definition of title page:

\maketitle 

\section{The Problem}\label{sec:problem}

A matrix equation $M\mathbf{u}=\mathbf{r}$ is just a way of writing N 
linear, algebraic equations. The $i$th equation is:
\begin{equation}
\sum_j{M_{ij}\,u_j}=r_i.
\end{equation}
The standard problem is to solve for $\mathbf{u}$ given $M$ and 
$\mathbf{b}$. A square matrix is said to be {\em tridiagonal} when
\begin{equation}
M_{ij}=0,\quad \left| i-j \right| > 1.
\end{equation}
For example,
\begin{equation}
\begin{pmatrix}\label{eq:4x4}
b_1 & c_1 & 0 & 0 \\
a_2 & b_2 & c_2 & 0 \\
0 & a_3 & b_3 & c_3 \\
0 & 0 & a_4 & b_4
\end{pmatrix}
\begin{pmatrix}
u_1 \\ u_2 \\ u_3 \\ u_4
\end{pmatrix}
=
\begin{pmatrix}
r_1 \\ r_2 \\ r_3 \\ r_4
\end{pmatrix}.
\end{equation}

\section{Algorithm}\label{sec:algorithm}

Have a look {\it Numerical Recipes} Section 2.6. The {\it Numerical Recipes}
function {\tt tridag} is easily understood without LU decomposition.

To solve the $4\times 4$ system in equation \ref{eq:4x4} for
$\mathbf{u}$, we begin by using row 1 to eliminate $a_2$ from row 2. 
Similarly, we eliminate all the $a_i$ row by row. 

Here is how the {\tt 
tridag} implementation works. The first row stands for the equation:

\begin{equation*}
b_1\,u_1+c_1\,u_2=r_1.
\end{equation*}

\begin{align}
\mbox{Let \quad} \beta_1 &= b_1,\\
v_1 &= \frac{r_1}{\beta_1},\\
\gamma_2 &= \frac{c_1}{\beta_1}
\end{align}
If we divide row 1 by $\beta_1$, its equation becomes
\begin{equation}\label{eq:1b}
u_1+\gamma_2\,u_2 = v_1.
\end{equation}

The equation for the second row is:
\begin{equation*}
a_2\,u_1+b_2\,u_2+c_2\,u_3 = r_2.
\end{equation*}

\begin{align*}
\mbox{Let\quad}
\beta_2 &= b_2 - a_2\,\gamma_2,\\
v_2 &= \frac{(r_2 - a_2\,v_1)}{\beta_2}, \\
\gamma_3 &= \frac{c_2}{\beta_2}.
\end{align*}
We now subtract $a_2$ times row 1 from row 2. Dividing the result by 
$\beta_2$, our new row 2 is written:
\begin{equation}
u_2 + \gamma_3\,u_3 = v_2.
\end{equation}

This process may be repeated over the range $i=2,3,\ldots N$, using 
the recursion formulas
\begin{align}
\gamma_i &= \frac{c_{i-1}}{\beta_{i-1}}, \\
\beta_i &= b_i - a_i\,\gamma_i, \\
v_i &= \frac{r_i - a_i\,v_{i-1}}{\beta_i}.
\end{align}
The result is that the matrix equation is rewritten in uper triangular 
form. For our $4 \times 4$ example,
\begin{equation}\label{eq:4x4b}
\begin{pmatrix}
1 & \gamma_2 & 0 & 0 \\
0 & 1 & \gamma_3 & 0 \\
0 & 0 & 1 & \gamma_4 \\
0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
u_1 \\ u_2 \\ u_3 \\ u_4
\end{pmatrix}
=
\begin{pmatrix}
v_1 \\ v_2 \\ v_3 \\ v_4
\end{pmatrix}.
\end{equation}
This problem is trivially solved by starting from the bottom and back 
substituting as you go up:
\begin{align}
u_N &= v_N,\\
u_i &= v_i - \gamma_{i+1}\,u_{i+1}.
\end{align}

\section{How to Break the Algorithm}\label{sec:break}

If $\beta_i=0$, for some $i$, the program quits to avoid a division by 
zero. The authors of {\it Numerical Recipes} assure us that this 
doesn't happen very often. 

It can be shown that $\beta_i = 0$ cannot occur in the calculation of
$y_i^{\prime\prime}$ for the cubic spline.

\end{document}
