I'm trying to write a function that can solve a tridiagonal system of linear equations using the Thomas algorithm. It basically solves the following equation. (Details can be found at the Wiki page here Tridiagonal matrix algorithm.)
$$ \begin{bmatrix} {b_ 1} & {c_ 1} & { } & { } & { 0 } \\ {a_ 2} & {b_ 2} & {c_ 2} & { } & { } \\ { } & {a_ 3} & {b_ 3} & \ddots & { } \\ { } & { } & \ddots & \ddots & {c_{n-1}}\\ { 0 } & { } & { } & {a_n} & {b_n}\\ \end{bmatrix} \cdot \begin{bmatrix} {x_ 1 } \\ {x_ 2 } \\ {x_ 3 } \\ \vdots \\ {x_n } \\ \end{bmatrix} = \begin{bmatrix} {d_ 1 } \\ {d_ 2 } \\ {d_ 3 } \\ \vdots \\ {d_n } \\ \end{bmatrix} $$
This can be done easily in various ways in Mathematica using the built-in functions such as Solve
, LinearSolve
, LUDecomposition
, etc. Instead of using Mathematica's wonderful black box solvers, I have decided to follow the procedures in the algorithm and write my own version which I can better control.
Here is the procedure in detail:
$$ c'_i= \begin{cases} \begin{array}{lcl} \cfrac{c_i}{b_i} &&; i = 1 \\ \cfrac{c_i}{b_i - c'_{i - 1} a_i} &&; i = 2, 3, \ \dots, n-1 \\ \end{array} \end{cases} $$
$$ d'_i= \begin{cases} \begin{array}{lcl} \cfrac{d_i}{b_i} &&; i = 1 \\ \cfrac{d_i - d'_{i - 1} a_i}{b_i - c'_{i - 1} a_i} &&; i = 2, 3, \dots, n. \\ \end{array} \end{cases} $$
$$ \begin{array}{lcl} x_n&=&d'_n\\ x_i&=&d'_i-c'_ix_{i+1} \end{array} $$
I tried to write the function in Mathematica like this:
tridag[{a_, b_, c_}, d_] :=
Module[{n = Length[d], c1 = Range[Length[d]], d1 = Range[Length[d]], x = Range[Length[d]]},
c1[[1]] = c[[1]]/b[[1]];
d1[[1]] = c1[[1]]/c[[1]]*d[[1]];
Do[c1[[i]] = c[[i]]/(b[[i]] - a[[i]]*c1[[i - 1]]);
d1[[i]] = c1[[i]]/c[[i]]*(d[[i]] - a[[i]] d1[[i - 1]]);, {i, 2, n - 1}];
x[[n]] = d1[[n]];
Do[x[[i]] = d1[[i]] - c1[[i]]*x[[i + 1]], {i, n - 1, 1, -1}];
x]
where a, b, c are the three diagonal rows in the matrix.
For a test case of
A = {{1, 2, 0, 0, 0}, {2, 2, 3, 0, 0}, {0, 3, 3, 4, 0}, {0, 0, 4, 4, 5}, {0, 0, 0, 5, 5}};
R = {5, 15, 31, 53, 45};
my function gives the same answer as Mathematica's linear solver:
LinearSolve[A, R]
{1, 2, 3, 4, 5}
tridag[{PadLeft[Diagonal[A, 1], 5], Diagonal[A], PadRight[Diagonal[A, -1], 5]}, R]
{1, 2, 3, 4, 5}
Here are my questions:
- Since I am just a novice, my function is simply a direct translation of the above procedurual equations. So I would like to learn a really functional Mathematica way of implementing this function (other than using the built-in solver)?
- A more general question is that, given an algorithm written in a procedure way, i.e., with implementation in a procedure language like C or Fortran in mind, how can it be translated into functional Mathematica code? Can you give some examples or suggestions or guidelines?
Edit
As belisarius points out, "'the real elegant Mathematica way' is always to use an already implemented function". To be clear, I'm not asking how to implement a solver that is superior to the built-in solver. I'm just asking how to implement a solver in the Mathematica way, i.e., in the functional way. I think knowing how to write one's own solver is helpful when using the built-in solvers. It gives the user confidence when the built-in solver does not work. An example may be the HilbertMatrix
:
LinearSolve[ HilbertMatrix [10], Table[Random[], {10}]]
LinearSolve::luc: "Result for LinearSolve of badly conditioned matrix {{1., 0.5, 0.333333, 0.25, 0.2, 0.166667, 0.142857, 0.125, 0.111111, 0.1}, <<8>>, {0. 1, 0.0909091, 0.0833333, 0. 0769231, 0.0714286, 0.0666667, 0.0625,0.0588235, 0.0555556, 0.0526316}} may contain significant numerical errors."
Besides, I think trying to implement a simple solver is a good way for me as a new hand to learn the Mathematica language.
m
be the matrix,Last /@ RowReduce[Join[m, {d}\[Transpose], 2]]
faithfully implements this algorithm. I suspect even @belisarius might not downvote this solution :-). If you insist the algorithm be expressed in terms of the vectorsa
,b
, andc
, then begin withm = SparseArray[{Band[{2, 1}] -> a, Band[{1, 1}] -> b, Band[{1, 2}] -> c}];
. $\endgroup$ – whuber Mar 9 '13 at 4:32RowReduce
is "black box." The purpose in offering that solution is to demonstrate that there's nothing special about the Thomas algorithm; it's just a special case of row reduction. At the end, whatever algorithm you code has to rely on some primitive operations, whether they are addition or row reduction or something in between. At this point it seems we're guessing a little concerning just how far to go and when to stop. $\endgroup$ – whuber Mar 9 '13 at 21:49