Solving a tridiagonal system of equations
|Simple Examples||Subroutines Examples||Dynamic Examples||Steady State Overall Mass Balance||Array and Matrix Examples|
|Multi-Step Trigger||Checking |
|Belt Filter |
PGM FILE - Solve a Tridiagonal System of equations
For PSD calculations and elsewhere, we often get tridiagonal systems of equations. We could solve these using matrix functions, but it is much more efficient to solve the system directly. Instead of having to store a full nxn matrix, we use three nx1 arrays a, b, c(and we don't actually use the first element of a or the last element of c. For an 8x8 system we use the following elements of the arrays:
[b0 c0 0 ... ] [a1 b1 c1 ] [0 a2 b2 c2 ... ] [ ... ] [ ... a7 b7]
For larger arrays, this saves a lot of space: if we have a tridiagonal system of size 100, we only need to store 300 elements, as opposed to storing 10000 elements for the full system. Solving a tridiagonal system is also very efficient, requiring just O(n) operations.
Unfortunately we can't pass array or matrix classes to functions, but we can still define a function to do our work, and reuse it in different PGM's. We have to define the arrays we use globally and then pas s the solution back in a global array.
We then need five arrays: a, b, c are used to define the diagonal and off-diagonal elements (b contains the elements on the diagonal). r contains the right hand side, and x is used to store the solution of the equation A x = r.
In this example we solve the nxn system:
[4 1 0 ... ]  [1 4 1 ]  [0 1 4 1 ... ] x =  [ ... ] . [ ... 1 4 1] . [ ... 1 4] [n]
The code here doesn't do any checks for divide by zero - for a lot of cases the matrix is diagonally dominant and we don't anticipate any problems with this simplified approach. YMMV.
;declares array instances. array a,b,c,x,r long n, sz*<<5>> ;define the function: function tridiag() long n, j real z array w ;; workspace n = a.GetLen() j = 1 w.SetLen(n) z = b ;; Assume b != 0, maybe check x = r/z while (j<n) ;; Forward Reduction w[j] = c[j-1] / z z = b[j] - a[j] * w[j] ;; check that z!=0 ! x[j] = (r[j] - a[j] * x[j-1]) / z j = j+1 endwhile j = n-2 while (j>=0) ;; Back Substitution x[j] = x[j] - w[j+1] * x[j+1] j = j-1 endwhile return 0 ;; if problems, return error code endfunct
You can save the TRIDIAG code in a PGM file (eg: tridiag.pgm) and then include it whenever you need it. See Include Files. To use the function, just set up the arrays and then call:
; >>tridiag.pgm ; use this if the function is written as an include file ; Set up arrays a.setlen(sz) b.setlen(sz) c.setlen(sz) r.setlen(sz) x.setlen(sz) a.setall(1.0) b.setall(4.0) c.setall(1.0) n = 0 while (n<sz) r[n] = n+1 n = n+1 endwhile ; Call the tridiagonal function we have set up above. tridiag() watch x[All, 10] $