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
 (34)
Introduce new unknowns ej and fj, along with an equation
 (35)
Write (35) with shifted index:
 (36)
Insert (36) into (14):
 (37)
Now rearrange (37) to resemble (15):
 (38)
Compare (38) to (35) to see recursions for the new unknowns ej and fj:
 (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, .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
 (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.
 (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 etc. The subroutine rtris() solves equation (33) for q where n=5, endl,endr,a=c, and b.

# 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: Finite-differencing in the time Up: FINITE DIFFERENCING IN (omega,x)-SPACE Previous: The Crank-Nicolson method
Stanford Exploration Project
12/26/2000