How to implement tridiagonal matrix algorithm?

Question

I have created the tridiagonal matrix below:

LTable={{-10.4151, 9.4098, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {105.526, -104.76, 94.4741, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 106.107, -104.511, 93.893, 0,0, 0, 0, 0, 0, 0}, {0, 0, 106.749, -104.396, 93.2507, 0, 0, 0, 0, 0,0}, {0, 0, 0, 107.459, -104.413, 92.5409, 0, 0, 0, 0, 0}, {0, 0, 0,0, 108.244, -104.564, 91.7564, 0, 0, 0, 0}, {0, 0, 0, 0, 0,109.111, -104.854, 90.8894, 0, 0, 0}, {0, 0, 0, 0, 0, 0,110.069, -105.29, 89.9312, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 111.128, -105.883, 88.8723, 0}, {0, 0, 0, 0, 0, 0, 0, 0,112.298, -106.643, 87.702}, {0, 0, 0, 0, 0, 0, 0, 0, 0, -11.1525,3.02486}};
RTable={1., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.};
Print[LTable // MatrixForm,Table[Subscript[y, i], {i, 1, n + 1}] // MatrixForm, "=", RTable // MatrixForm]

I know I can solve this using the built-in LinearSolve function, but I want to see too how can I implement the calculation of tridiagonal matrix in Mathematica in terms of For loops.

Answer 1

Using the Gauss algorithm for $LU$-factorization on banded matrices will result in banded $L$ and $U$ factors with the same bandwith. Thus, it suffices to work on the diagonals. If there are only 3 diagonals, then the algorithm reduces to the Thomas algorithm.

Here is a Mathematica implementation that uses Compile to produce runtime optimized libraries for machined real-valued tridiagonal matrices:

cThomasLUDecomposition = Compile[{{l0, _Real, 1}, {d0, _Real, 1}, {u0, _Real, 1}},
   Block[{n, u, uk, d, dold, l, lk},
    n = Length[d0];
    d = Table[0., {n}];
    u = Table[0., {n}];
    l = Table[0., {n}];
    
    d[[1]] = dold = Compile`GetElement[d0, 1];
    
    Do[
     u[[k]] = uk = Compile`GetElement[u0, k];
     l[[k]] = lk = Compile`GetElement[l0, k] / dold;
     d[[k + 1]] = dold = Compile`GetElement[d0, k + 1] - lk uk;
     ,
     {k, 1, n - 1}];
    
    {l, d, u}
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

cThomasLUSolve = Compile[{{l, _Real, 1}, {d, _Real, 1}, {u, _Real, 1}, {y, _Real, 1}},
   Block[{n, x, xold},
    n = Length[d];
    x = Table[0., {n}];
    
    (*Solving x = LinearSolve[L,y];*)
    x[[1]] = xold = Compile`GetElement[y, 1];
    Do[
     x[[k]] = xold = Compile`GetElement[y, k] - Compile`GetElement[l, k - 1] xold
     , {k, 2, n}];
    
    (*Solving x = LinearSolve[U,x];*)
    x[[n]] = xold = Compile`GetElement[x, n]/Compile`GetElement[d, n];
    Do[
     x[[k]] = xold = (Compile`GetElement[x, k] - Compile`GetElement[u, k] xold)/Compile`GetElement[d, k]
     , {k, n - 1, 1, -1}];
    x
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Here a usage example:

A = {{-10.4151, 9.4098, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {105.526, -104.76,
     94.4741, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 106.107, -104.511, 93.893, 
    0, 0, 0, 0, 0, 0, 0}, {0, 0, 106.749, -104.396, 93.2507, 0, 0, 0, 
    0, 0, 0}, {0, 0, 0, 107.459, -104.413, 92.5409, 0, 0, 0, 0, 
    0}, {0, 0, 0, 0, 108.244, -104.564, 91.7564, 0, 0, 0, 0}, {0, 0, 
    0, 0, 0, 109.111, -104.854, 90.8894, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 
    110.069, -105.29, 89.9312, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 
    111.128, -105.883, 88.8723, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 
    112.298, -106.643, 87.702}, {0, 0, 0, 0, 0, 0, 0, 0, 0, -11.1525, 
    3.02486}};

(*Extracting the diagonals.*)
a = Diagonal[A, -1];
b = Diagonal[A, 0];
c = Diagonal[A, +1];

(*Compute the factorization only one.*)
{l, d, u} = cThomasLUDecomposition[a, b, c];

(*Now you can use that factorization to solve with as many right-hand sides as you like.*)

y = RandomReal[{-1, 1}, Length[b]];
x = cThomasLUSolve[l, d, u, y];

Max[Abs[A . x - y]]

2.84217*10^-14

By the way, the factors $L$ and $U$ can be obtained as follows:

U = DiagonalMatrix[SparseArray[d]] + DiagonalMatrix[SparseArray[Most[u]], +1];
L = DiagonalMatrix[SparseArray[ConstantArray[1., Length[d]]]] + DiagonalMatrix[SparseArray[Most[l]], -1];
Max[Abs[L . U - A]]

2.84217*10^-14

Answer 2

I found a good reference for the algorithm. It used two approaches for the input: vector and matrix.

https://tamaskis.github.io/files/Tridiagonal_Matrix_Algorithm.pdf

Answer 3

Here is a quick and dirty solution.

We solve the last equation or row from hand. Then we recursively solve the second to last, Third to last.. until the second row. The first row we again solve by hand.

ars = Table[Symbol["y" <> ToString[i]], {i, 11}];
sol = {y10 -> -LTable[[11, 11]]/LTable[[11, 10]] y11};
Do[
 eq = (LTable[[i]] );
 eq = (eq[[i - 1 ;; i + 1]] . vars[[i - 1 ;; i + 1]]) /. sol;
 sol = Join[sol, Solve[eq == 0, vars[[i - 1]]][[1]]];
 , {i, 10, 2, -1}]
sol0 = Solve[(LTable[[1, 1 ;; 2]] . {y1, y2} /. sol) == 1, y11][[1]];
sol = sol /. sol0;
sol = PrependTo[sol, sol0[[1]]];
res = vars /. so

{0.120071, 0.239171, 0.131094, -0.124365, -0.289298, -0.182, \
0.133879, 0.372935, 0.272769, -0.141348, -0.521142}

To check if the solution: res is correct:

LTable . res // Chop

{1., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}