next up previous print clean
Next: Finite-differencing in the time Up: FINITE DIFFERENCING IN (omega,x)-SPACE Previous: The Crank-Nicolson method

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  
 \begin{displaymath}
a_j \, q_{j+1} \ +\ b_j \, q_j \ +\ c_j \, q_{j-1}
 \eq d_j\end{displaymath} (34)
Introduce new unknowns ej and fj, along with an equation  
 \begin{displaymath}
q_j \eq e_j \, q_{j+1} \ +\ f_j\end{displaymath} (35)
Write (35) with shifted index:  
 \begin{displaymath}
q_{j-1} \eq e_{j-1} \, q_j \ +\ f_{j-1}\end{displaymath} (36)
Insert (36) into (14):  
 \begin{displaymath}
a_j \, q_{j+1} \ +\ b_j \, q_j \ +\ 
c_j \, ( e_{j-1} \, q_j \ +\ f_{j-1} ) \eq d_j\end{displaymath} (37)
Now rearrange (37) to resemble (15):  
 \begin{displaymath}
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} }\end{displaymath} (38)
Compare (38) to (35) to see recursions for the new unknowns ej and fj:
      \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 e0 and f0. With e0 and all the aj , bj , cj, we can use (39) to compute all the ej.

On the right-hand boundary we need a boundary condition. The general two-point boundary condition is  
 \begin{displaymath}
c_{n-1} \, q_{n-1} \ +\ e_{\rm right} \, q_n \eq d_n\end{displaymath} (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.  
 \begin{displaymath}
q_{n-1} \eq e_{n-1} \, q_n \ +\ f_{n-1}\end{displaymath} (42)
Both qn and qn-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 qn 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 ej depends on the medium through aj, bj, cj, but it does not depend on the solution qj (even through dj). This means that it may be possible to save and reuse ej. (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).


next up previous print clean
Next: Finite-differencing in the time Up: FINITE DIFFERENCING IN (omega,x)-SPACE Previous: The Crank-Nicolson method
Stanford Exploration Project
12/26/2000