\ 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 ;