Thomas Algorithm

As alluded to in the Gaussian Elimination chapter, the Thomas Algorithm (or TDMA, Tri-Diagonal Matrix Algorithm) allows for programmers to massively cut the computational cost of their code from to in certain cases! This is done by exploiting a particular case of Gaussian Elimination where the matrix looks like this:

This matrix shape is called Tri-Diagonal (excluding the right-hand side of our system of equations, of course!). Now, at first, it might not be obvious how this helps. Well, firstly, it makes the system easier to encode: we may divide it into four separate vectors corresponding to , , , and (in some implementations, you will see the missing and set to zero to get four vectors of the same size). Secondly, and most importantly, equations this short and regular are easy to solve analytically.

We'll start by applying mechanisms familiar to those who have read the Gaussian Elimination chapter. Our first goal is to eliminate the terms and set the diagonal values to . The and terms will be transformed into and . The first row is particularly easy to transform since there is no , we simply need to divide the row by :

Let's assume that we found a way to transform the first rows. How would we transform the next one? We have

Let's transform row in two steps.

Step one: eliminate with the transformation :

Step two: get with the transformation :

Brilliant! With the last two formula, we can calculate all the and in a single pass, starting from row , since we already know the values of and .

Of course, what we really need are the solutions . It's back substitution time!

If we express our system in terms of equations instead of a matrix, we get

plus the last row that is even simpler: . One solution for free! Maybe we can backtrack from the last solution? Let's (barely) transform the above equation:

and that's all there is to it. We can calculate all the in a single pass starting from the end.

Overall, we only need two passes, and that's why our algorithm is ! The transformations are quite easy too, isn't that neat?

Example Code

function(a::Vector{Float64}, b::Vector{Float64}, c::Vector{Float64},
         d::Vector{Float64}, soln::Vector{Float64})
    # Setting initial elements
    c[0] = c[0] / b[0]
    d[0] = d[0] / b[0]

    for i = 1:n
        # Scale factor is for c and d
        scale = 1 / (b[i] - c[i-1]*a[i])
        c[i] = c[i] * scale
        d[i] = (d[i] - a[i] * d[i-1]) * scale
    end

    # Set the last solution for back-substitution
    soln[n-1] = d[n-1]

    # Back-substitution
    for i = n-2:0
        soln[i] = d[i] - c[i] * soln[i+1]
    end

end
#include <stdio.h>
#include <string.h>

void thomas(double * const a, double * const b, double * const c,
            double * const x, const size_t size) {

    double y[size];
    memset(y, 0, size * sizeof(double));

    y[0] = c[0] / b[0];
    x[0] = x[0] / b[0];

    for (size_t i = 1; i < size; ++i) {
        double scale = 1.0 / (b[i] - a[i] * y[i - 1]);
        y[i] = c[i] * scale;
        x[i] = (x[i] - a[i] * x[i - 1]) * scale;
    }

    for (int i = size - 2; i >= 0; --i) {
        x[i] -= y[i] * x[i + 1];
    }
}

int main() {
    double a[] = {0.0, 2.0, 3.0};
    double b[] = {1.0, 3.0, 6.0};
    double c[] = {4.0, 5.0, 0.0};
    double x[] = {7.0, 5.0, 3.0};

    printf("The system,\n");
    printf("[1.0  4.0  0.0][x] = [7.0]\n");
    printf("[2.0  3.0  5.0][y] = [5.0]\n");
    printf("[0.0  3.0  6.0][z] = [3.0]\n");
    printf("has the solution:\n");

    thomas(a, b, c, x, 3);

    for (size_t i = 0; i < 3; ++i) {
        printf("[%f]\n", x[i]);
    }

    return 0;
}
# Author: gammison

# note this example is inplace and destructive
def thomas(a, b, c, d):

    #set the initial elements
    c[0] = c[0] / b[0]
    d[0] = d[0] / b[0]

    n = len(d) #number of equations to solve
    for i in range(1,n):
        #scale factor for c and d
        scale = 1 / (b[i] - c[i-1]*a[i])

        c[i] = c[i] * scale
        d[i] = (d[i] -a[i] * d[i-1]) * scale


    # do the back substitution
    for i in range(n-2,-1,-1):
        d[i] = d[i] - c[i]*d[i+1]

    return d

def main():
    # example for matrix
    # [1  4  0][x]   [7]
    # [2  3  5][y] = [5]
    # [0  3  6][z]   [3]
    #                 [.8666]
    # soln will equal [1.533]
    #                 [-.266]
    # note we index a from 1 and c from 0
    a = [0, 2, 3]
    b = [1, 3, 6]
    c = [4, 5, 0]
    d = [7, 5, 3]
    soln = thomas(a, b, c, d)
    print(soln)

if __name__ == '__main__':
    main()
public class thomas {
    private static void thomasAlgorithm(double a[], double b[], double c[], double x[], int size) {

        double y[] = new double[size];

        y[0] = c[0] / b[0];
        x[0] = x[0] / b[0];

        for (int i = 1; i < size; ++i) {
            double scale = 1.0 / (b[i] - a[i] * y[i - 1]);
            y[i] = c[i] * scale;
            x[i] = (x[i] - a[i] * x[i - 1]) * scale;
        }

        for (int i = size - 2; i >= 0; --i) {
            x[i] -= y[i] * x[i + 1];
        }
    }

    public static void main(String[] args) {
        double a[] = {0.0, 2.0, 3.0};
        double b[] = {1.0, 3.0, 6.0};
        double c[] = {4.0, 5.0, 0.0};
        double x[] = {7.0, 5.0, 3.0};

        System.out.println("The system,\n");
        System.out.println("[1.0  4.0  0.0][x] = [7.0]\n");
        System.out.println("[2.0  3.0  5.0][y] = [5.0]\n");
        System.out.println("[0.0  3.0  6.0][z] = [3.0]\n");
        System.out.println("has the solution:\n");

        thomasAlgorithm(a, b, c, x, 3);

        for (int i = 0; i < 3; ++i)
            System.out.println("[" + x[i] + "]\n");
    }
}
import Data.List (zip4)
import Data.Ratio

thomas :: Fractional a => [a] -> [a] -> [a] -> [a] -> [a]
thomas a b c = init . scanr back 0 . tail . scanl forward (0, 0) . zip4 a b c
  where
    forward (c', d') (a, b, c, d) =
      let denominator = b - a * c'
       in (c / denominator, (d - a * d') / denominator)
    back (c, d) x = d - c * x

main :: IO ()
main = do
  let a = [0, 2, 3] :: [Ratio Int]
      b = [1, 3, 6]
      c = [4, 5, 0]
      d = [7, 5, 3]
  print $ thomas a b c d

results matching ""

    No results matching ""