5
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ Check UpperTriangularize? $\endgroup$ – Öskå Apr 7 '14 at 18:45
  • $\begingroup$ That just sets all the elements below the diagonal to 0. I actually want the first steps of Gaussian elimination to be applied. $\endgroup$ – rioneye Apr 7 '14 at 18:49
  • 5
    $\begingroup$ LUDecomposition? $\endgroup$ – Öskå Apr 7 '14 at 19:13
  • $\begingroup$ LU decomposition is a bit different than what I am hoping for. I would like to get an upper triangular matrix using Gaussian elimination (mathworld.wolfram.com/GaussianElimination.html). $\endgroup$ – rioneye Apr 7 '14 at 20:24
  • 1
    $\begingroup$ Look at RowReduce... $\endgroup$ – ciao Apr 7 '14 at 22:02
7
$\begingroup$

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)}} *)
| improve this answer | |
$\endgroup$
1
$\begingroup$

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
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.