Simplify matrix into an upper triangular matrix
I'll illustrate on a smallish matrix, using exact arithmetic. That should make it relatively easy to verify correctness.
First we create a tridiagonal 5×5 matrix.
n = 5;
SeedRandom[33333];
mat = RandomInteger[{-100, 100}, {n, n}];
Do[mat[[i, j]] = 0; mat[[j, i]] = 0, {i, 3, n}, {j, 1, i - 2}];
mat
(*Out[68]= {{-30, -98, 0, 0, 0}, {12, 72, 29, 0, 0}, {0, -9, -49, 63,
0}, {0, 0, -21, 88, -16}, {0, 0, 0, -16, -98}} *)
To make it upper triangular we simply clear below pivots. I'm sure this can be done more slickly with Fold
, but I'm more used to procedural methods for this sort of thing.
Do[mat[[i]] -= mat[[i, i - 1]]/mat[[i - 1, i - 1]]*mat[[i - 1]], {i, 2, n}]
mat
(* Out[70]= {
{-30, -98, 0, 0, 0},
{0, 164/5, 29, 0, 0},
{0, 0, -(6731/164), 63, 0},
{0, 0, 0, 375356/6731, -16},
{0, 0, 0, 0, -(9627006/93839)}} *)
It is not always advised to perform LU decomposition without pivoting, for tridiagonal matrices or otherwise. Nevertheless, whenever it is stable (e.g. diagonally dominant cases), there is a simple formula for the $\mathbf L$ and $\mathbf U$ factors, based on the usual three-term recurrence for the determinant of a tridiagonal matrix (see e.g. this).
In particular, consider the tridiagonal matrix
MatrixForm[trid = With[{n = 5},
SparseArray[{Band[{2, 1}] -> Array[l, n - 1],
Band[{1, 1}] -> Array[d, n],
Band[{1, 2}] -> Array[u, n - 1]}, {n, n}]]]
$$\begin{pmatrix} d[1] & u[1] & 0 & 0 & 0 \\ l[1] & d[2] & u[2] & 0 & 0 \\ 0 & l[2] & d[3] & u[3] & 0 \\ 0 & 0 & l[3] & d[4] & u[4] \\ 0 & 0 & 0 & l[4] & d[5] \\ \end{pmatrix}$$
We can use RecurrenceTable[]
to evaluate the three-term recurrence for the tridiagonal determinant:
With[{n = 5},
rr = RecurrenceTable[{y[n] == d[n] y[n - 1] - l[n - 1] u[n - 1] y[n - 2],
y[-1] == 0, y[0] == 1}, y[n], {n, 0, n}]];
Check:
Last[rr] == Det[trid] // Simplify
True
Here, however, we are interested in using rr
to build the desired decomposition. In particular, the required factors are:
With[{n = 5},
lf = SparseArray[{Band[{2, 1}] -> Array[l, n - 1]/Ratios[Most[rr]],
Band[{1, 1}] -> 1}, {n, n}];
uf = SparseArray[{Band[{1, 1}] -> Ratios[rr],
Band[{1, 2}] -> Array[u, n - 1]}, {n, n}];]
Check:
lf.uf == trid // Simplify
True
Let's apply this to Daniel's example:
la = {12, -9, -21, -16};
da = {-30, 72, -49, 88, -98};
ua = {-98, 29, 63, -16};
RecurrenceTable[]
does not play well when the diagonals are specified as lists, so I will demonstrate a technique for evaluating the three-term recurrence, based on repeated multiplication of $2\times 2$ matrices:
n = Length[da];
rr = FoldList[Dot, IdentityMatrix[2],
Transpose[{{ConstantArray[0, n], ConstantArray[1, n]},
{-Append[la ua, 0], da}}, {2, 3, 1}]][[All, 2, 2]]
{1, -30, -984, 40386, 2252136, -231048144}
(I have previously used this in this blog post.)
The required factors are then
MatrixForm[lf = SparseArray[{Band[{2, 1}] -> la/Ratios[Most[rr]],
Band[{1, 1}] -> 1}, {n, n}]]
$$\begin{pmatrix} 1 & 0 & 0 & 0 & 0 \\ -\frac{2}{5} & 1 & 0 & 0 & 0 \\ 0 & -\frac{45}{164} & 1 & 0 & 0 \\ 0 & 0 & \frac{3444}{6731} & 1 & 0 \\ 0 & 0 & 0 & -\frac{26924}{93839} & 1 \\ \end{pmatrix}$$
MatrixForm[uf = SparseArray[{Band[{1, 1}] -> Ratios[rr],
Band[{1, 2}] -> ua}, {n, n}]]
$$\begin{pmatrix} -30 & -98 & 0 & 0 & 0 \\ 0 & \frac{164}{5} & 29 & 0 & 0 \\ 0 & 0 & -\frac{6731}{164} & 63 & 0 \\ 0 & 0 & 0 & \frac{375356}{6731} & -16 \\ 0 & 0 & 0 & 0 & -\frac{9627006}{93839} \\ \end{pmatrix}$$
Note that uf
is the matrix Daniel obtained in his answer.
A final check:
lf.uf == SparseArray[{Band[{2, 1}] -> la, Band[{1, 1}] -> da,
Band[{1, 2}] -> ua}, {n, n}] // Simplify
True