Solving a tridiagonal system of equations

From SysCAD Documentation
Jump to navigation Jump to search

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 (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]
$