MATH2071: LAB 3: Solving linear systems

Introduction Exercise 1
Test matrices Exercise 2
The linear system ``problem'' Exercise 3
Gauss factorization Exercise 4
Permutation matrices Exercise 5
PLU factorization Exercise 6
PLU solution Exercise 7
A system of ODEs Exercise 8
  Exercise 9
  Exercise 10
  Exercise 11


In many numerical applications, notably the ordinary and partial differential equations this semester, the single most time-consuming step is the solution of a system of linear equations. We now begin a major topic in Numerical Analysis, Numerical Linear Algebra. We assume you are familiar with the concepts of vectors and matrices, matrix multiplication, vector dot products, and the theory of the solution of systems of linear equations, the matrix inverse, eigenvalues and so on. We now concentrate on the practical aspects of carrying out the computations required to do linear algebra in a numerical setting.

This lab will take four sessions. If you print this lab, you may prefer to use the pdf version.

We begin by noting three special test matrices that we will refer to when trying out algorithms. Then we state the linear system problem and consider three methods of solution, using the determinant, the inverse matrix, or Gauss factorization. We will find that factorization into three simply-solved factors is the best way to go, and we will write a m-file to perform the factorization and another to solve systems given the factorization. We will end up with an example using our m-files as part of the numerical solution of a partial differential equation.

Test matrices

We will now consider a few simple test matrices that we can use to study the behavior of the algorithms for solution of a linear system. We can use such problems as a quick test that we haven't made a serious programming error, but more importantly, we also use them to study the accuracy of the algorithms. We want test problems which can be made any size, and for which the exact values of the determinant, inverse, or eigenvalues can be determined. Matlab itself offers a ``rogues' gallery'' of special test matrices through the Matlab function gallery, but we will mainly be using our own m-files for these test matrices.

The three test we will be using are:

The second difference matrix

The second difference matrix arises in approximating the second derivative on equally spaced data. The formula is:

$\displaystyle {\bf A}_{i,j}=\left\{\begin{array}{cl}
-2 & \mbox{ for }j=i 1 &...
...{ for }\vert j-i\vert=1\\
0 & \mbox{ for }\vert j-i\vert>1 \end{array}\right.

So that the matrix looks like:

$\displaystyle \left( \begin{array}{rrrrrr}
-2 & 1 & & && \\
1 & -2 & 1 & && ...
...& 1 & & & & \ddots & &&\\
& & 1 &-2 &1 \\
&& & 1 &-2 &
\end{array} \right) $

where blanks represent zero values. The matrix is negative definite ($ -A$ is positive definite), symmetric, and tridiagonal. The determinant is plus or minus $ n+1$ where $ n$ denotes the size ($ A$ is $ n\times n$). Matlab provides a function to generate this matrix called gallery('tridiag',n), but we will use the m-file dif2.m. (Do not use the gallery('tridiag',n) form in these labs. Matlab has two forms of storage for matrices, ``sparse'' and ``full.'' Everything we have done so far has used the ``full'' form, but gallery('tridiag',n) returns the ``sparse'' form. Because ``sparse'' matrix form requires some special handling, we will not be using it for this lab.)

The Frank matrix

Row $ i$ of the $ n\times n$ Frank matrix has the formula:

$\displaystyle {\bf A}_{i,j}=\left\{ \begin{array}{cl}
0 & \mbox{ for }j<i-2\\
n+1-i & \mbox{ for }j=i-1\\
n+1-j & \mbox{ for }j\ge i

The Frank matrix for $ n=5$ looks like:

$\displaystyle \left(\begin{array}{ccccc}

The determinant of the Frank matrix is 1, but is difficult to compute numerically. This matrix has a special form called Hessenberg form wherein all elements below the first subdiagonal are zero. Matlab provides the Frank matrix in its ``gallery'' of matrices, gallery('frank',n), but we will use an m-file frank.m. It turns out that the inverse of the Frank matrix also has a simple formula.

The Hilbert matrix

The Hilbert matrix is related to the interpolation problem on the interval [0,1]. The matrix is given by the formula $ {\bf A}_{i,j}=1/(i+j-1)$. For example, with $ n=5$ the Hilbert matrix is:

$\displaystyle \left(\begin{array}{ccccc}
\frac11& \frac12& \frac13& \frac14& \f...
...\frac18 [5pt]
\frac15& \frac16& \frac17& \frac18& \frac19

The Hilbert matrix arises in interpolation and approximation contexts because it happens that $ {\bf A}_{i,j}=\int_0^1 (x^{i-1})(x^{j-1})dx$. The Hilbert matrix is at once nice because its inverse has integer elements and also not nice because it is extremely difficult to compute the inverse using the usual formulæ for matrix inverses.

Matlab has special functions for the Hilbert matrix and its inverse, called hilb(n) and invhilb(n), but we will use the m-file hilbert.m.

Download all three m-files, dif2.m, frank.m, and hilbert.m now.

The linear system ``problem''

The linear system ``problem'' can be posed in the following way. Find a $ n$-vector $ \bf x$ that satisfies the matrix equation

$\displaystyle {\bf A}{\bf x}={\bf b}$ (1)

where $ \bf A$ is a $ m\times n$ matrix and $ \bf b$ is a $ m$-vector.

You probably already know that there is ``usually" a solution if the matrix is square, (that is, if $ m=n$). We will concentrate on this case for now. But you may wonder if there is an intelligent response to this problem for the cases of a square but singular matrix, or a rectangular system, whether overdetermined or underdetermined.

At one time or another, you have probably been introduced to several algorithms for producing a solution to the linear system problem, including Cramer's rule (using determinants), constructing the inverse matrix, Gauss-Jordan elimination and Gauss factorization. We will see that it is usually not a good idea to construct the inverse matrix and we will focus on Gauss factorization for solution of linear systems.

We will be concerned about several topics:

Efficiency: what algorithms produce a result with less work?
Accuracy: what algorithms produce an answer that is likely to be more accurate?
Difficulty: what makes a problem difficult or impossible to solve?
Special cases: how do we solve problems that are big? symmetric? banded? singular? rectangular?

The inverse matrix

A classic result from linear algebra is the following:
Theorem The linear system problem (1) is uniquely solvable for arbitrary $ \bf b$ if and only if the inverse matrix $ {\bf A}^{-1}$ exists. In this case the solution $ \bf x$ can be expressed as $ {\bf x}={\bf A}^{-1}{\bf b}$.

So what's the catch? There are a few:

In the following exercises, you will see that constructing the inverse matrix takes time proportional to $ n^3$ (for an $ n\times n$ matrix), as does simply solving the system, but that solving the system is several times faster than constructing the inverse matrix. You will also see that matrices with special formats, such as tridiagonal matrices, can be solved very efficiently. And you will see even simple matrices can be difficult to numerically invert.

You will be measuring elapsed time in order to see how long these calculations take. The problem sizes are designed to take a modest amount of time on the computers in GSCC.

Exercise 1: The Matlab command inv(A) computes an approximation for the inverse of a matrix A. Do not worry for now just how it does it, but be assured that it uses one of the most efficient and reliable methods available.

Matlab provides the commands tic and toc to measure computational time. The tic command starts the timer, and the toc command stops it, either printing the elapsed time or returning its value as in the expression elapsedTime=toc;. The times are in seconds. (Note: tic and toc measure elapsed time. When a computer is doing more than one thing, the cputime function can be used to measure how much computer time is being used for this task alone.)

  1. Copy the following code to a Matlab function m-file named exer1.m and modify it to produce information for the table below. Be sure to add comments and your name and the date. Include the file with your summary.
    function elapsedTime=exer1(n)
    % elapsedTime=exer1(n)
    % comments
    % your name and the date
    A = dif2(n);
    b = ones(n,1);   % the right side vector doesn't change the time
    Ainv =  ???      % compute the inverse matrix
    xSolution = ???  % compute the solution

  2. You have seen in lecture that the time required to invert an $ n\times n$ matrix should be proportional to $ n^3$. Fill in the following table, where the column entitled ``Ratio'' should contain the ratio of the time for n divided by the time for the following value of n. (Note: timings are not precise! Furthermore, the first time a function is used results in Matlab reading the m-file and ``compiling it,'' and that takes considerable time.)
      Time to compute inverse matrices
       n     time      ratio
      160  _________  _________
      320  _________  _________
      640  _________  _________
     1280  _________  _________
     2560  _________
  3. Are these solution times roughly proportional to $ n^3$?

Exercise 2: Matlab provides a special operator, the backslash (\) operator, that is designed to solve a linear system without computing the inverse. It is used in the following way, for a matrix A and column vector b.
It may help you to remember this operator if you think of the A as being ``underneath'' the slash. The effect of this operator is to find the solution of the system of equations A*x=b.
  1. Copy exer1.m to exer2.m and replace the inverse matrix computation and solution with an expression using the Matlab ``\'' command, and fill in the following table
      Time to compute solutions
       n     time      ratio
      160  _________  _________
      320  _________  _________
      640  _________  _________
     1280  _________  _________
     2560  _________
  2. Are these solution times roughly proportional to $ n^3$?
  3. Compare the times for the inverse and solution and fill in the following table
      Comparison of times
       n   (time for inverse)/(time for solution)
      160               _________
      320               _________
      640               _________
     1280               _________
     2560               _________
  4. Theory shows that computation of a matrix inverse should take approximately three times as long as computation of the solution. Are your results consistent with theory?

You should be getting the message: ``You should never compute an inverse when all you want to do is solve a system.''

Warning: the ``\'' symbol in Matlab will work when the matrix A is not square. In fact, sometimes it will work when A actually is a vector. The results are not usually what you expect and no error message is given. Be careful of this potential for error when using ``\''. A similar warning could be made for the / (division) symbol because Matlab will try to ``interpret'' it if you use it with matrices or vectors and you will get ``answers'' that are not what you intend.

Exercise 3: Some matrices consist of mostly zero entries. These so-called ``sparse'' matrices sometimes admit extremely efficient methods for solving them. Matlab has a special storage scheme for sparse matrices, and the Matlab function gallery('tridiag',N) yields a sparse matrix that is the negative of the result from dif2(N).
  1. For N=10 and b=ones(N,1), use the ``\'' command to solve the system $ Ax=b$ once using the matrix from dif2(N) and once using the matrix from gallery('tridiag',N). Confirm that the solutions are negatives of each other. The easiest way to confirm that two vectors are the same is to compute the norm of their difference.
  2. Using the same methodology as in Exercise 2, fill in the following table using the matrix from gallery('tridiag',N). (Note that sparse matrix storage allows much larger matrices.)
      Time to compute solutions
      Size    time      ratio
      2560  _________  _________
      5120  _________  _________
     10240  _________  _________
     20480  _________  _________
     40960  _________  _________
     81920  _________
  3. Because we are now solving a sparse system, solution time is not necessarily $ O(n^3)$. Estimate $ p$ so that sparse solution time is $ O(n^p)$.
  4. Compare the times for N=2560 for the sparse matrix computation and for the full matrix computation in Exercise 2 above. You should see a big advantage for sparse storage. One method of sparse storage will be discussed in a later lab.

Exercise 4: Some nonsingular matrices result in systems that are very difficult to solve numerically. You should recall that the condition number, defined as the size of the largest eigenvalue over the smallest eigenvalue is a measure of the liklihood of roundoff errors in solving a matrix equation. In this exercise you will construct a system by picking a solution vector $ x_{\text{known}}$ and then constructing the right side $ b=Ax_{\text{known}}$. Then the known solution can be compared with the solution computed using the ``\'' command to see how bad a particular solution can be. Consider the following Matlab code:
A = dif2(N);
xKnown = ones(N,1);
b = ???          % compute the right side corresponding to xKnown
xSolution = ???  % compute the solution using "\"

  1. For the matrix defined by dif2, fill in the following table. (Recall that the Matlab command ``cond'' yields the condition number and the command ``det'' yields the determinant.)
            dif2 matrix
     size    error   determinant cond. no.
       10  _________ __________  __________
       40  _________ __________  __________
      160  _________ __________  __________
    The dif2 matrix is about as nicely behaved as a matrix can be. It is neither close to singular nor difficult to invert.
  2. For the matrix defined by hilbert, fill in the following table.
           Hilbert matrix
     size    error   determinant cond. no.
       10  _________ __________  __________
       15  _________ __________  __________
       20  _________ __________  __________
    The Hilbert matrix is so close to being singular that numerical approximations cannot distinguish it from a singular matrix.
  3. For the matrix defined by frank, fill in the following table.
              Frank matrix
     size    error   determinant cond. no.
       10  _________ __________  __________
       15  _________ __________  __________
       20  _________ __________  __________
    The Frank matrix can be proved to have a determinant equal to 1.0 because its eigenvalues come in $ (x,1/x)$ pairs. Thus, it is not close to singular. Nonetheless, roundoff errors conspire to make it appear to be singular on a computer. The condition number can alert you to this possibility. You can find some information about the Frank matrix from the Matrix Market, and the references therein.

In the following sections, you will write your own routine that does the same task as the Matlab ``\'' command.

Gauss factorization

The standard method for computing the inverse of a matrix is Gaussian elimination. You already know this algorithm, perhaps by another name such as row reduction. You have seen it discussed in class and you can find a discussion on the web at, for example, or at

First, we will prefer to think of the process as having two separable parts, the ``factorization'' or ``forward substitution'' and ``back substitution'' steps. The important things to recall include:

  1. In the $ k^{\mbox{th}}$ step of factorization, one of the remaining equations (rows) is chosen to become the ``pivot,'' and to be associated with $ x_k$. This equation (row) is swapped with the existing $ k^{\mbox{th}}$ equation (row).
  2. The pivot equation is used to eliminate references to $ x_k$ in the remaining equations (i.e., put zeros in position $ k$ in the rows below the $ k^{\mbox{th}}$ row). After this step, the matrix will have zeros in all positions below and to the left of the diagonal $ (k,k)$ position.
  3. After $ n-1$ steps of factorization, back substitution begins. The last equation involves a single variable, and can be solved for $ x_n$. But then equation $ n-1$ can be solved for $ x_{n-1}$, because $ x_n$ is now known, and similarly, the other equations are solved in backwards order.

The simplest version of Gauss factorization for a matrix $ A$ with entries $ A_{k,\ell}$ is called Gauss factorization with no pivoting, because it has a very simple way of choosing the pivot equation. The pivot on for the $ k^{\mbox{th}}$ row is the diagonal entry $ A_{k,k}$. If on step $ k$, the diagonal is zero the method will fail, even though a solution may actually exist. The forward substitution steps begin at the upper left of the matrix and proceed down the diagonal, converting all the entries below the diagonal to zero. The backward substitution steps begin at the lower right and work their way back up the diagonal.

The method of ``Gauss factorization with partial pivoting'' chooses as pivot the largest value $ A_{k,\ell}$ for $ \ell\ge k$. To keep the calculation orderly, this pivot row is actually moved into the $ k^{\mbox{th}}$ row of the matrix. In this case, the matrix gradually is transformed into upper triangular shape. For exact arithmetic, if there is a solution, this method is guaranteed to reach it. Moreover, for inexact arithmetic, this method has better accuracy properties than no pivoting.

The method of ``Gauss factorization with complete pivoting'' chooses as pivot for the $ k^{\mbox{th}}$ row the coefficient $ A_{i,j}$ of largest magniture for $ i\ge k$ and $ j\ge k$, This method has more bookkeeping than the partial pivoting method, and doesn't produce much improvement, so it is little used.

We now turn to writing code for the Gauss factorization. The discussion in this lab assumes that you are familiar with Gaussian elimination (sometimes called ``row reduction''), and is entirely self-contained. Because, however, interchanging rows is a source of confusion, we will assume at first that no row interchanges are necessary.

The idea is that we are going to use the row reduction operations to turn a given matrix $ \bf A$ into an upper triangular matrix (all of whose elements below the main diagonal are zero) $ \bf U$ and at the same time keep track of the multipliers in the lower triangular matrix $ \bf L$. These matrices will be built in such a way that $ \bf A$ could be reconstructed from the product $ \bf LU$ at any time during the factorization procedure.

Alert: The fact that $ {\bf A}={\bf LU}$ even before the computation is completed can be a powerful debugging technique. If you get the wrong answer at the end of a computation, you can test each step of the computation to see where it has gone wrong. Usually you will find an error in the first step.

Let's do a simple example first. Consider the matrix

$\displaystyle {\bf A}=\left(\begin{array}{cc}2&4 1&9\end{array}\right).

There is only one step to perform in the row reduction process: turn the 1 in the lower left into a 0 by subtracting half the first row from the second. If $ {\bf A}={\bf L}{\bf U}$, convince yourself that

$\displaystyle {\bf L}=\left(\begin{array}{cc}1&0 \frac12&1\end{array}\right)


$\displaystyle {\bf U}=\left(\begin{array}{cc}2&4 0&7\end{array}\right)

are matrices that describe the row reduction step ($ \bf L$) and its result ($ \bf U$) in such a way that the original matrix $ \bf A$ can be recovered as $ {\bf A}={\bf L}{\bf U}$.

Notice that the matrix $ \bf U$ is what you would obtain from a row reduction operation applied to $ \bf A$. Also notice that

$\displaystyle \mathbf{L}^{-1}=\left(\begin{array}{cc}1&0 -\frac12&1\end{array}\right)$

so that the product $ \mathbf{L}^{-1}\mathbf{A}$ can be described as ``take half the first row of $ \bf A$ and subtract it from the second row of $ \bf A$.''

Now, suppose you have an 5 by 5 matrix (switching to Matlab notation) A that is the transpose of the Frank matrix, and you wish to perform row reduction steps to turn all the entries in the first column below the first row into zero, and keep track of the operations in the matrix L and the result in U. Convince yourself that the following code does the trick.

A=frank(5)';  % TRANSPOSE: so A is not too easy
L=eye(n);     % square n by n identity matrix
for Irow=Jcol+1:n
  % compute Irow multiplier and save in L(Irow,Jcol)

  % multiply row "Jcol" by L(Irow,Jcol) and subtract from row "Irow"


Exercise 5:
  1. Use the above code to compute the first row reduction steps for the matrix A. Print the resulting matrix U and include it in the summary report. Check that all entries in the first column in the second through last rows are zero. (If any are non-zero, you have an error somewhere. Find the error before proceeding.)
  2. Multiply L*U and confirm for yourself that it equals A. An easy way to do this is to compute norm(L*U-A,'fro'). (The Frobenius norm is a little faster to compute than the default 2-norm, especially for large matrices.)
  3. Use the code given above as a model for an m-file called gauss_lu.m that performs Gaussian reduction without pivoting. The function should start out
    function [L,U]=gauss_lu(A)
    % function [L,U]=gauss_lu(A)
    % performs an LU factorization of the matrix A using
    % Gaussian reduction.
    % A is the matrix to be factored.
    % L is a lower triangular factor with 1's on the diagonal
    % U is an upper triangular factor.
    % A = L * U
    % your name and the date
  4. Use your routine to find L and U for the transpose of the Frank matrix of order 5. Include L and U in your summary.
  5. Verify that L is lower triangular (has zeros above the diagonal) and U is upper triangular (has zeros below the diagonal).
  6. Confirm that L*U recovers the original matrix.
  7. The Matlab expression R=rand(50,50) will generate a $ 50\times50$ matrix of random entries. Use gauss_lu to find its factors LR and UR.
    • Use norms to confirm that LR*UR=R.
    • Use tril and triu to confirm that LR is lower triangular and UR is upper triangular.
    • Do not send me copies of the matrices. (You should not even print the matrices out!)

It turns out that omitting pivoting leaves the function you just wrote vulnerable to failure.

Exercise 6:
  1. Compute L1 and U1 for the matrix
      A1 = [ -2   1   0   0   0
              1   0   1  -2   0
              0   0   0   1  -2
              1  -2   1   0   0
              0   1  -2   1   0]
    The method should fail for this matrix.
  2. On which step of the decomposition (values of Irow and Jcol) does the method fail? Why does it fail?
  3. What is the determinant of A1? What is the condition number (cond(A1))? Is the matrix singular or close to singular?

In Gaussian elimination it is entirely possible to end up with a zero on the diagonal, as the previous exercise demonstrates. Since you cannot divide by zero, the procedure fails. It is also possible to end up with small numbers on the diagonal and big numbers off the diagonal. In this case, the algorithm doesn't fail, but roundoff errors can overwhelm the procedure and result in wrong answers. One solution to these difficulties is to switch rows and columns around so that the largest remaining entry in the matrix ends up on the diagonal. Instead of this ``full pivoting'' strategy, it is almost as effective to switch the rows only, instead of both rows and columns. This is called ``partial pivoting.'' In the following section, permutation matrices are introduced and in the subsequent section they are used to express Gaussian elimination with partial pivoting as a ``PLU'' factorization.

Permutation matrices

In general, a matrix that looks like the identity matrix but with two rows interchanged is a matrix that causes the interchange of those two rows when multiplied by another matrix. A permutation matrix can actually be a little more complicated than that, but all permutation matrices are products of matrices with a single interchanged pair of rows.

Exercise 7: Consider the two permutation matrices
  P1=[1 0 0 0
      0 0 1 0
      0 1 0 0
      0 0 0 1];
  P2=[1 0 0 0
      0 0 0 1
      0 0 1 0
      0 1 0 0];
  1. Compute the product of P1 with the magic square matrix of order 4, A=magic(4), A2=P1*A. If the result is not the matrix A with rows two and three interchanged, then you made an error. Please include the matrix A2 in your summary.
  2. What permutation matrix would interchange the first and fourth rows of A?
  3. Compute the product P=P1*P2. The product of two permutation matrices is again a permutation matrix, so P is also a permutation matrix, but it is slightly more complicated than before. Please include the matrix P in your summary.
  4. Compute the product A3=P*A and notice that it is the original A with mixed-up rows. In general, a permutation matrix is one in which each row and each column has exactly one 1 in it and all other entries are zero. Please include the matrix A3 with your summary.

PLU factorization

You already have written gauss_lu.m to do Gauss eliminiation without pivoting. Now, you are going to add pivoting to the process. The previous section should have convinced you that pivoting is the same thing as multiplying by a permutation matrix. Just as in gauss_lu.m, where you started out with U equal to the identity matrix and saved the row operations in it, you will start off with a permutation matrix P=eye(n) and save the row interchanges in it.

Exercise 8:
  1. Copy your file gauss_lu.m to a file called gauss_plu.m. Change the signature line to
    function [P,L,U]=gauss_plu(A)
    and change the comments to describe the new output matrix P.
  2. Add the initialization statement P = eye(n) near the beginning.
  3. Add the following code to the beginning of the loop on Jcol, (i.e., just before the for Irow= ... statement. This code will be performed before any other processing inside that loop.
      % first, choose a pivot row by finding the largest entry
      % on or below position k on the diagonal
      [ colMax, pivotShifted ] = max ( abs ( U(Jcol:n, Jcol ) ) );
      % The value of pivotShifted from max needs to be adjusted.
      % See help max for more information.
      pivotRow = Jcol+pivotShifted-1;
      if pivotRow ~= Jcol     % no pivoting if pivotRow==Jcol
        U([Jcol, pivotRow], :)       = U([pivotRow, Jcol], :);
        L([Jcol, pivotRow], 1:Jcol-1)= L([pivotRow, Jcol], 1:Jcol-1);
        P(:,[Jcol, pivotRow])        = P(:,[pivotRow, Jcol]);
    The reason that pivotShifted must be adjusted to get pivotRow is because, for example, if J(jcol+1,jcol) happens to be the largest entry in the column, then the max function will return pivotShifted=2), not jcol+1. This is because (jcol+1) is the second entry in the vector (Jcol:n).
  4. Temporarily add the following statement at the end of the loop on Jcol (not the loop on Irow)
      if norm(P*L*U-A,'fro')> 1.e-12 * norm(A,'fro')
        error('If you see this message, there is a bug!')
  5. Test gauss_plu.m on the matrix A=magic(5), and check that P*L*U=A. Please include P, L, and U in your summary.
  6. Test gauss_plu.m on the matrix A1 from Exercise 6 above, and check that P1*L1*U1=A1. Recall that this was for matrix for which gauss_lu.m failed. You should not get any error messages.
  7. Remove the ``TEMPORARY DEBUGGING CHECK'' that you inserted above. This check involves too much computation and will mess up the timings in later exercises.

PLU solution

We have seen how to factor a matrix $ \bf A$ into the form:

$\displaystyle {\bf A} = {\bf P} {\bf L} {\bf U}

where $ \bf P$ is a permutation matrix, $ \bf L$ is a unit lower triangular matrix, and $ \bf U$ is an upper triangular matrix. To solve a matrix equation involving $ \bf A$, we use this factor information to ``peel away" the multiplication a step at a time:
$\displaystyle {\bf Ax}$ $\displaystyle =$ $\displaystyle {\bf b}$  
$\displaystyle {\bf P} {\bf L} {\bf U} {\bf x}$ $\displaystyle =$ $\displaystyle {\bf b}$  
$\displaystyle {\bf L} {\bf U} {\bf x}$ $\displaystyle =$ $\displaystyle {\bf P}^{-1} {\bf b}$  
$\displaystyle {\bf U} {\bf x}$ $\displaystyle =$ $\displaystyle {\bf L}^{-1}({\bf P}^{-1} {\bf b})$  
$\displaystyle {\bf x}$ $\displaystyle =$ $\displaystyle {\bf U}^{-1}({\bf L}^{-1}{\bf P}^{-1} {\bf b})$  

However, instead of explicitly computing the inverses of the matrices, we use special facts about their form. It's not the inverse matrix itself that we want, but the inverse times the right hand side.

Let's consider the factored linear system once more. It helps to make up some names for variables and look at the problem in a different way:

$\displaystyle {\bf P}({\bf L} {\bf U} {\bf x})$ $\displaystyle =$ $\displaystyle {\bf b}$  
$\displaystyle {\bf P}({\bf z})$ $\displaystyle =$ $\displaystyle {\bf b}$  
$\displaystyle {\bf L}({\bf U} {\bf x})$ $\displaystyle =$ $\displaystyle {\bf z}$  
$\displaystyle {\bf L}({\bf y})$ $\displaystyle =$ $\displaystyle {\bf z}$  
$\displaystyle {\bf U} {\bf x}$ $\displaystyle =$ $\displaystyle {\bf y}$  
$\displaystyle {\bf U} {\bf x}$ $\displaystyle =$ $\displaystyle {\bf y}$  
$\displaystyle {\bf x}$ $\displaystyle =$ $\displaystyle {\bf U}^{-1}{\bf y}$  
$\displaystyle {\bf x}$ $\displaystyle =$ $\displaystyle {\bf U}^{-1}{\bf L}^{-1}{\bf P}^{-1} {\bf b}$  

Notice that all I have done is to use parentheses to group factors, and name them. But you should now be able to see what an algorithm for solving this problem might look like.

In summary, the PLU Solution Algorithm is: To solve $ {\bf A}{\bf x}={\bf b}$ after finding the factors $ \bf P$, $ \bf L$, $ \bf U$ and the right hand side $ \bf b$,

  1. Solve $ {\bf P}{\bf z}={\bf b}$;
  2. Solve $ {\bf L}{\bf y}={\bf z}$; and,
  3. Solve $ {\bf U}{\bf x}={\bf y}$.

Ouch! I promised a simple algorithm, but now instead of having to solve one system, we have to solve three, and we have three different solution vectors running around. But things actually have gotten better, because these are really simple systems to solve:

  1. The solution of $ {\bf P}{\bf z}={\bf b}$ is $ {\bf z}={\bf P}^T{\bf b}$;
  2. It is easy to solve $ {\bf L}{\bf y}={\bf z}$ because $ \bf L$ is lower triangular and you can start at the top left; and
  3. It is easy to solve $ {\bf U}{\bf x}={\bf y}$ because $ \bf U$ is upper triangular and you can start at the bottom right.
The next exercise illustrates this process on a simple matrix and the following exercise has you write a Matlab routine to carry out the process in general. It also provides a worked-out example you can use to check Matlab code.

Exercise 9: Consider the simple matrix factorization
$\displaystyle {\bf A}$ $\displaystyle =$ $\displaystyle {\bf P}{\bf L}{\bf U}$  
$\displaystyle \left(\begin{array}{cc}1&9 2&4\end{array}\right)$ $\displaystyle =$ $\displaystyle \left(\begin{array}{cc}0&1 1&0\end{array}\right)
...&0 0.5&1\end{array}\right)
\left(\begin{array}{cc}2&4 0&7\end{array}\right)$ (2)

  1. Using the following matrices,
    P=[0  1
       1  0];
    L=[1  0
      0.5 1];
    U=[2  4
       0  7];
    confirm that $ {\bf A}={\bf PLU}$.
  2. Using pencil and paper, go through the steps of solving the system, with right hand side b=[19; 10] (note the semicolon to make b a column vector), using the PLU factorization. This will serve as a test problem for the code.
        0    b = [          19;         10 ]
        1    z = [  __________; __________ ]
        2    y = [  __________; __________ ]
        3    x = [  __________; __________ ]

In the following exercise, you will be writing a routine plu_solve.m, to solve the linear system (1) by solving the permutation, lower and upper triangular systems. Because this is a three stage process, we will divide it into four separate m-files: u_solve.m, l_solve.m, p_solve.m, and plu_solve.m. The plu_solve routine will look like:

function x = plu_solve ( P, L, U, b )
% function x = plu_solve ( P, L, U, b )
% solve a factored system 
% P=permutation matrix
% L=lower-triangular
% U=upper-triangular
% b=right side
% x=solution

% your name and the date

% Solve P*z=b for z using p_solve
z = p_solve(P,b);

% Solve L*y=z for y using l_solve
y = l_solve(L,z);

% Solve U*x=y for x using u_solve
x = u_solve(U,y);

Exercise 10:
  1. The easiest of these routines by far is p_solve because the matrix $ \bf P$ is orthogonal so its solution is $ {\bf z}={\bf P}^T{\bf b}$. This is a one-liner in Matlab. Replace the question marks with code in the following outline.
    function z=p_solve(P,b)
    % z=p_solve(P,b)
    % P is an orthogonal matrix
    % b is the right hand side
    % z is the solution of P*z=b
    % your name and the date
    z= ???
  2. Test your routine using the $ \bf P$ matrix and $ \bf b$ vector from the previous exercise. Do not continue until you are confident your file is correct.
  3. The second routine is still not so hard because all you have to do is start at the top and run down the rows. Replace the question marks with code in the following outline.
    function y=l_solve(L,z)
    % y=l_solve(L,z)
    % L is a lower-triangular matrix whose diagonal entries are 1
    % z is the right hand side
    % y is the solution of L*y=z
    % your name and the date
    % set n for convenience and simplicity
    % initialize y to zero and make sure it is a column vector
    % first row is really an easy one, especially since the diagonal
    % entries of L are equal to 1
    for Irow=2:n
      rhs = z(Irow);
      for Jcol=1:Irow-1
        rhs = rhs - ???
      y(Irow) = ???
  4. Test this routine using the $ \bf L$ matrix and $ \bf z$ vector from the previous exercise. Do not continue until you are confident your file is correct.
  5. The third routine a little harder because you have to start at the bottom and run up the rows. Replace the question marks with code in the following outline.
    function x = u_solve(U,y)
    % x=u_solve(U,y)
    % U is an upper-triangular matrix
    % y is the right hand side
    % x is the solution of U*x=y
    % your name and the date
    % set n for convenience and simplicity
    % initialize y to zero and make sure it is a column vector
    % last row is the easy one
    % the -1 in the middle means it is going up from the bottom
    for Irow=n-1:-1:1
      rhs = y(Irow);
      % the following loop could also be written as a single
      % vector statement.
      for Jcol=Irow+1:n
        rhs = rhs - ???
      x(Irow) = ???
  6. Test this routine using the $ \bf U$ matrix and $ \bf y$ vector from previous exercise. Do not continue until you are confident your file is correct.
  7. Now put them all together and verify your solution by calculating
      x = plu_solve(P,L,U,b)
    for the matrices and right side vector in the previous exercise.
  8. You have written and checked your code for a $ 2\times2$ system. Be sure it works for a $ 3\times3$ system using the following code.
      A = dif2(3);
      x = [1;1;1];
      b = A*x;
      [P L U]=gauss_plu(A);
      xSolution = plu_solve ( ??? );
      relErr = norm(x-xSolution)/norm(x)
    What are P, L, U, and xSolution? You should see the error is roundoff-sized.
  9. As a final test, solve a large system with randomly generated coefficients and right side. What is relErr?
      A = rand(100,100);
      x = rand(100,1);
      b = A*x;
    and the rest of the code as in the previous part of this Exercise. The error relErr will usually be a roundoff sized number. There is always some chance that the randomly generated matrix A will end up poorly conditioned or singular, though. If that happens, generate A and x again: they will be different this time.

You might wonder why we don't combine the m-files gauss_plu.m and plu_solve.m into a single routine. One reason is that there are other things you might want to do with the factorization instead of solving a system. For example, one of the best ways of computing a determinant is to factor the matrix. (The determinant of any permutation matrix is $ \pm1$, the determinant of L is 1 because L has all diagonal entries equal to 1, and the determinant of U is the product of its diagonal entries.) Another reason is that the factorization is much more expensive (i.e., time-consuming) than the solve is, and you might be able to solve the system many times using the same factorization. This latter case is presented in the following section.

A system of ODEs

You saw how to numerically solve certain ordinary differential equations in Lab 1. In this section you will see that the implicit Euler method can be used to solve a system of ordinary differential equations (or, what is the same thing, an ordinary differential equation involving vectors and matrices). At each time step in that solution, you need to solve a linear system of equations, and you will use plu_factor and plu_solve for this task. You will also see the advantage of splitting the solution into a ``factor'' step and a ``solve'' step.

Consider a system of ordinary differential equations.

$\displaystyle \frac{d\mathbf{x}}{dt}=\mathbf{Ax}$ (3)

where $ \mathbf{A}$ is an $ n\times n$ matrix and $ \mathbf{x}$ is an $ n$-vector. Since $ \mathbf{x}$ is a vector, it has subscripts, but since we are solving an ODE numerically, there will be time steps as well. We will write time levels using a superscript here so that it is a little clearer. Applying the implicit Euler method to (3) yields the system of equations

\begin{displaymath}\begin{array}{rcl} \mathbf{x}^{k+1}&=&\mathbf{x}^k + \Delta t...
...^k &=&(\mathbf{I}+\Delta t\mathbf{A})\mathbf{x}^k \end{array}\end{displaymath} (4)

where $ \mathbf{I}$ denotes the identity matrix.

In Matlab, we will represent the vector $ \mathbf{x}^k$ as x(:,k), where the ``:'' means ``all the subscripts'' and is equivalent to x(1:end,k). In addition, the Matlab function eye (that is a pun-get it?) generates the identity matrix.

The following code embodies the procedure described above. Copy it into a script m-file called back_euler_system.m. This particular system, using the dif2 matrix, turns out to be a discretization of the time-dependent heat equation. The variable $ \mathbf{x}$ represents temperature and is a function of time and also of an (here unspecified) one-dimensional spatial variable that is a multiple of its subscript. The plot is a compact way to illustrate the vector $ \mathbf{x}$ with its subscript measured along the horizontal axis.

% ntimes = number of temporal steps from t=0 to t=1 .
% dt = time increment

% N is the dimension of the space
% initial values (sawtooth)
x(1:N/2,1) = (1:N/2)';
x(N/2+1:N,1) = (N/2+1:N)'-N;
% discretization matrix is dif2 multiplied
% by a constant to make the solution interesting

for k = 1 : ntimes
  t(1,k+1) = t(1,k) + dt;
  x(:,k+1) = plu_solve(P,L,U,x(:,k));

% plot the solution at sequence of times
hold on
for k=2:2:ntimes
title 'Solution at selected times'
hold off

Exercise 11:
  1. Use the above routine to solve the ODE system. Record the time this solution takes. Save your solution for comparison by copying it (xsave=x). Note: The size of this problem is chosen to take a minute or two on the computers in GSCC. It will take longer if you are using a slower computer.
  2. Next, change the m-file so that you use the routine gauss_plu.m only once outside the loop to factor the matrix EulerMatrix into P, L, and U and then use plu_solve.m to solve the system inside the loop. Do the problem again and record the time. The time spent for this case should be a fraction of the time with an inversion each time step because factorization takes a lot more arithmetic than solution.
  3. Compare the two solutions by computing the norm of the difference (norm(xsave-x,inf)) to show they agree. Include your revised version of back_euler_system.m with your summary.

You should be aware that the system matrix in this exercise is sparse, but we have not taken advantage of that fact here. A lot of time could be saved if we modified the code to take advantage of sparseness.

If you happened to try Exercise 11 using the ``\'' form, you would find that it is substantially faster than our routines. It is faster because it takes advantage of programming efficiencies that are beyond the scope of the course, not because it uses different methods. One such efficiency would be the replacement of as many loops as possible with vector statements.

Back to MATH2071 page.

Mike Sussman 2009-01-07