Solving a tridiagonal system of equations

Navigation: PGMs ➔ Example PGM Files ➔ Tridiagonal System

Simple Examples		Subroutines Examples					Dynamic Examples		Steady State Overall Mass Balance			Array and Matrix Examples
Basic Layout	Simple Calculations	Initialise PreStart	Multi-Step Trigger	Checking Project	Counter, While and Random	Belt Filter Wash Loss	Startup Actions	Mass Balance	Mass Balance	Species Balance	Elemental Balance	Lookup Value	Set Values	Tridiagonal System

Related Links: PGM Files using Class and Functions, Matrix Class, Array Class

---

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]
 [1 4 1         ]     [2]
 [0  1 4 1 ...  ] x = [3] 
 [        ...   ]      .
 [    ...  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[0]    ;; Assume b[0] != 0, maybe check
    x[0] = r[0]/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 (e.g.: 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]
$