\ Linear equations with tridiagonal matrices by LU method
\ reference: Press, et al., "Numerical Recipes" (Cambridge
U. Press, 1986)
\ ---------------------------------------------------
\ (c) Copyright 1999 Julian V. Noble. \
\ Permission is granted by the author to \
\ use this software for any application pro- \
\ vided this copyright notice is preserved. \
\ ---------------------------------------------------
\ This is an ANS Forth program requiring the
\ FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets.
\
\ Environmental dependences:
\ Assumes independent floating point stack
FALSE [IF]
Algorithm:
Ax = r, A is tridiagonal
b(k), k=1,...,n are the diagonal elements
a(k), k=2,...,n are the lower subdiagonal elements
c(k), k=1,...,n-1 are the upper subdiagonal elements
Write A = LU where L is lower-triangular and U upper triangular
If A were a general matrix this would be possible also but the
decomposition would require O(n^3) steps. For the tridiagonal
case, however, the decomposition requires only O(n) steps.
Can assume L and U are bi-diagonal. In the case of L the lower
subdiagonal is just a(k). In the case of U we can choose the
diagonal elements to be 1 (unity). Thus we need to determine
the diagonal elements of L and the upper subdiagonal of U.
Call them L(k) and U(k) respectively. Then
L(1) = b(1)
U(k) = c(k)/L(k), k=1,...,n-1
L(k) = b(k) - a(k)*U(k-1), k=2,...,n
Finally let Ux = y and solve first
Ly = r
via
y(1) = r(1)/L(1)
y(2) = [r(2) - a(2)*y(1)] / L(2) etc.
then
x(n) = y(n)
x(n-1) = y(n-1) - x(n)*U(n-1) etc.
Usage:
Say a{ b{ c{ n }triangularize r{ x{ n }backsolve
[THEN]
MARKER -tridiag
BL PARSE undefined DUP PAD C! PAD CHAR+ SWAP CHARS MOVE PAD FIND NIP 0=
[IF] : undefined BL WORD FIND NIP 0= ; [THEN]
include arrays.f
include ftran111.f
20 VALUE Nmax
Nmax long 1 FLOATS 1array a{ \ input array in vector form
Nmax long 1 FLOATS 1array b{
Nmax long 1 FLOATS 1array c{
0 VALUE aa{ 0 VALUE bb{ 0 VALUE cc{ 0 VALUE NN
Nmax long 1 FLOATS 1array r{ \ inhomogeneous term
Nmax long 1 FLOATS 1array L{ \ diagonal of lower-triangular matrix
Nmax long 1 FLOATS 1array U{ \ subdiagonal of upper-triangular matrix
Nmax long 1 FLOATS 1array x{ \ solution vector
: }triangularize ( a{ b{ c{ n --)
TO NN TO cc{ TO bb{ TO aa{
f" bb{ 0 }" FDUP F0= ABORT" Reduce # of equations by 1"
L{ 0 } F!
f" U{ 0 } = cc{ 0 } / L{ 0 }"
NN 1- 0 DO f" U{ I } = cc{ I } / L{ I }"
f" L{I_1+} = bb{I_1+} - aa{I_1+} * U{I}"
LOOP
;
: }backsolve ( r{ x{ n --)
TO NN TO aa{ TO bb{
f" bb{0} = bb{0} / L{0}"
NN 1 DO f" bb{ I } = (bb{ I } - a{I}*bb{I_1-}) / L{I}"
LOOP
f" aa{ NN_1- } = bb{ NN_1- }"
0 NN 2 - DO f" aa{I} = bb{I} - U{I}*aa{I_1+}"
LOOP
;