Solving tridiagonal simultaneous equations

Much of the world's scientific computing power gets used up solving tridiagonal simultaneous equations. For reference and completeness the algorithm is included here.

Let the simultaneous equations be written as a difference equation  

$$a_j \, q_{j+1} \ +\ b_j \, q_j \ +\ c_j \, q_{j-1} \eq d_j$$ (34)

Introduce new unknowns e_j and f_j, along with an equation  

$$q_j \eq e_j \, q_{j+1} \ +\ f_j$$ (35)

Write (35) with shifted index:  

$$q_{j-1} \eq e_{j-1} \, q_j \ +\ f_{j-1}$$ (36)

Insert (36) into (14):  

$$a_j \, q_{j+1} \ +\ b_j \, q_j \ +\ c_j \, ( e_{j-1} \, q_j \ +\ f_{j-1} ) \eq d_j$$ (37)

Now rearrange (37) to resemble (15):  

$$q_j \eq { - \ a_j \over b_j \ +\ c_j \, e_{j-1} } \ \ q_{j+... ... +\ \ { d_j \ -\ c_j \, f_{j-1} \over b_j \ +\ c_j \, e_{j-1} }$$ (38)

Compare (38) to (35) to see recursions for the new unknowns e_j and f_j:

$$\begin{eqnarray} e_j \ \ \ & = & \ \ \ { - \ a_j \over b_j \ +\ c_j \, e_{j-1} }... ... \ \ \ { d_j \ -\ c_j \, f_{j-1} \over b_j \ +\ c_j \, e_{j-1} }\end{eqnarray}$$ (39) (40)

First a boundary condition for the left-hand side must be given. This may involve one or two points. The most general possible end condition is a linear relation like equation (35) at j=0, namely, $q_0 \ =\ e_0 q_1 + f_0$.Thus, the boundary condition must give us both e₀ and f₀. With e₀ and all the a_j , b_j , c_j, we can use (39) to compute all the e_j.

On the right-hand boundary we need a boundary condition. The general two-point boundary condition is  

$$c_{n-1} \, q_{n-1} \ +\ e_{\rm right} \, q_n \eq d_n$$ (41)

Equation (41) includes as special cases the zero-value and zero-slope boundary conditions. Equation (41) can be compared to equation (36) at its end.  

$$q_{n-1} \eq e_{n-1} \, q_n \ +\ f_{n-1}$$ (42)

Both q_n and q_n-1 are unknown, but in equations (41) and (42) we have two equations, so the solution is easy. The final step is to take the value of q_n and use it in (36) to compute $q_{n-1} , \ q_{n-2} , \ q_{n-3} , \$ etc. The subroutine rtris() solves equation (33) for q where n=5, endl$=e_{\rm left}$,endr$=e_{\rm right}$,a=c$=-\alpha$, and b$=1-2\alpha$. 

# real tridiagonal equation solver
subroutine rtris( n, endl, a, b, c, endr, d, q)
integer i,        n
real q(n), d(n), a, b, c, den, endl, endr
temporary real f(n), e(n)
e(1) = -a/endl;   f(1) = d(1)/endl
do i= 2, n-1 {
        den = b+c*e(i-1);   e(i) = -a/den;   f(i) = (d(i)-c*f(i-1))/den
        }
q(n) = (d(n)-c*f(n-1)) / (endr+c*e(n-1))
do i= n-1, 1, -1
        q(i) = e(i) * q(i+1) + f(i)
return;         end

If you wish to squeeze every last ounce of power from your computer, note some facts about this algorithm. (1) The calculation of e_j depends on the medium through a_j, b_j, c_j, but it does not depend on the solution q_j (even through d_j). This means that it may be possible to save and reuse e_j. (2) In many computers, division is much slower than multiplication. Thus, the divisor in (19a,b) can be inverted once (and perhaps stored for reuse).