Simplify matrix into an upper triangular matrix

Question

I have a tridiagonal matrix that I am trying to simplify into an upper triangular matrix using Mathematica, so I can use back substitution and solve my linear system. I have found the command RowReduce, but that simplifies the matrix too far, and I get an identity matrix, since the matrix I am trying to simplify is a square matrix. Is there a way in Mathematica to simplify a matrix into an upper triangular matrix?

Note: I want to solve the linear system by hand. I just want Mathematica to simplify the matrix.

Update: Specifically, I'd like to get an upper triangular matrix by using Gaussian elimination.

Answer 1

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)}} *)

Answer 2

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