Back to tutorial index

Solving non-singular linear systems


Topics in this lab

Introduction

In this lab, we will learn how to use Matlab to solve linear systems of the form $$Ax = b$$, where the matrix $A$ is an $N \times N$ matrix, and $b$ is an $N \times 1$ vector. The matrix $A$ is assumed to be non-singular and so the linear system has exactly one solution $x$.

clear all;
format short g;

Back to the top

Backslash operator : x = A\b

Linear systems arise in all areas of engineering, science and computational mathematics, and as a result, significant research has gone into finding the most computationally efficient algorithms for solving linear systems. In this tutorial, we will investigate some standard solvers for square, non-singular systems. In a later tutorial, we will investigate the solution to over-determined and under-determined systems.

The first linear system that most students encounter is the one whose solution describes the point at the intersection of two lines. For example, suppose we have the equations for two lines, given by $$ \begin{aligned} 3x - y = 1 \\ 2x + y = 7 \end{aligned} $$ Using matrix notation, we can write this as the $2 \times 2$ linear system given by $$ \left[\begin{array}{rr} 3 & -1 \\ 2 & 1 \end{array}\right] \left[\begin{array}{r} x \\ y \end{array}\right] = \left[\begin{array}{r} 1 \\ 7 \end{array}\right]. $$ More generally, a linear system is written as $$A{\bf x} = {\bf b}$$ where $A$ is an $n \times n$ matrix, and ${\bf x}$ and ${\bf b}$ are both $n \times 1$ vectors. Formally, the solution to a linear system can be written as $${\bf x} = A^{-1}{\bf b}$$ In Matlab, we can solve this system using the "backslash" operator.

A = [3 -1; -2 -1]
A =

     3    -1
    -2    -1

b = [1; -7];
x = A\b
x =

          1.6
          3.8

Checking this answer, we see that the residual error for this problem is essentially zero

r = b - A*x
r =

  -4.4409e-16
   8.8818e-16

Back to the top

The LU Decomposition

Behind the scenes, what is the "backslash" operator doing? If $A$ is a square, non-singular matrix, the backslash operator uses Gaussian Elimination, a procedure very similar to what you might do by hand to solve a linear system. Basically, the computer eliminates variables in order, until only a single equation for the last variable remains. Then, working backwards, each variable is solved in order. This procedure is carried out using the $LU$ decomposition.

Essentially, the backslash operator is doing the following behind the scenes :

format short

A = [1,2,3; 4,-1,5; 3,1,0];
b = [1;1;1];
[L,U,P] = lu(A);

Here, $L$ is a lower triangular matrix, $U$ is an upper triangular matrix and $P$ is a permutation matrix. These three matrices are related as $$ PA = LU $$
These matrices are

L
L =

    1.0000         0         0
    0.2500    1.0000         0
    0.7500    0.7778    1.0000

U
U =

    4.0000   -1.0000    5.0000
         0    2.2500    1.7500
         0         0   -5.1111

P*A
ans =

     4    -1     5
     1     2     3
     3     1     0

L*U
ans =

    4.0000   -1.0000    5.0000
    1.0000    2.0000    3.0000
    3.0000    1.0000         0

We can then use these to solve for $x$ as follows.

format long
y = L\(P*b);    % "Forward substitution" - P*A = P*b
x = U\y         % "Backward substitution"
x =

   0.239130434782609
   0.282608695652174
   0.065217391304348

We can compare this solution to the solution using the backslash operator.

A\b
ans =

   0.239130434782609
   0.282608695652174
   0.065217391304348

The solutions are exactly the same.

Back to the top

Using linsolve and mldivide

Most Matlab operators have equivalent functions that can be used when using an operator ("backslash" in this case) is inconvenient. For example, you cannot pass a "backslash" operator as an argument to a function. In this case, the linsolve or the mldivide functions may be more convenient ways to obtain essentially the same solution as that obtained using backslash.

A = [1,2,3; 4,-1,5; 3,1,0];
b = [1;1;1];

x = linsolve(A,b)  % Return the same solution as backslash
x =

   0.239130434782609
   0.282608695652174
   0.065217391304348

A second function, mldivide (for "matrix left divide") is exactly the functional equivalent of backslash

A\b - mldivide(A,b)
ans =

     0
     0
     0

Back to the top

Solving sparse linear systems

If the matrix $A$ is stored in sparse matrix format, then the backslash operator will take advantage of the sparse format and solve the system much more quickly that the standard full matrix solve. To see this, will solve the trivial linear system $Ix = b$ where $I$ is the indentity matrix.

I = eye(5000);
b = ones(5000,1);

We will time both solves using tic and toc and see which solve is faster.

% Time the solve using full matrix.
tic
x = I\b;
toc
Elapsed time is 0.081546 seconds.

Create a sparse matrix that stores only the non-zero entries.

I = sparse(eye(5000));

tic
x = I\b;
toc
Elapsed time is 0.000076 seconds.

The sparse format is orders of magnitude times faster than the full format. This is due to the fact that the sparse format is able to avoid most multiplication by zeros.

An easy way to create banded sparse matrices is to use the spdiags command.

Back to the top

Why don't we just use inv?

It is very tempting to use the inv function to solve a linear system. After all, when we write down the solution to a linear system using mathematical notation, we write $${\bf x} = A^{-1}{\bf b}$$ So, it seems natural to use the inv function.

A = [1,2,3; 4,-1,5; 3,1,0];
b = [1;1;1];

x = inv(A)*b
x =

   0.239130434782609
   0.282608695652174
   0.065217391304348

For this simple problem, the solution is identical to what we got above using backslash. But there are examples of matrices for which the solution using inv will be much worse than that using the backslash operator. For this reason, the inv should be avoided.

Here is an example of the famous Hilbert matrix, introduced by the German mathematician David Hilbert in 1894. This matrix is an example of particular ill-conditioned matrix, meaning that it is very hard to solve a linear systems accurately that involves this matrix.

In the following example, we compute the solution to a Hilbert system using both backslash and the inverse operator, and compare the results to that obtained using the exact solution.

format short
A = hilb(10)
b = ones(10,1);
A =

  Columns 1 through 7

    1.0000    0.5000    0.3333    0.2500    0.2000    0.1667    0.1429
    0.5000    0.3333    0.2500    0.2000    0.1667    0.1429    0.1250
    0.3333    0.2500    0.2000    0.1667    0.1429    0.1250    0.1111
    0.2500    0.2000    0.1667    0.1429    0.1250    0.1111    0.1000
    0.2000    0.1667    0.1429    0.1250    0.1111    0.1000    0.0909
    0.1667    0.1429    0.1250    0.1111    0.1000    0.0909    0.0833
    0.1429    0.1250    0.1111    0.1000    0.0909    0.0833    0.0769
    0.1250    0.1111    0.1000    0.0909    0.0833    0.0769    0.0714
    0.1111    0.1000    0.0909    0.0833    0.0769    0.0714    0.0667
    0.1000    0.0909    0.0833    0.0769    0.0714    0.0667    0.0625

  Columns 8 through 10

    0.1250    0.1111    0.1000
    0.1111    0.1000    0.0909
    0.1000    0.0909    0.0833
    0.0909    0.0833    0.0769
    0.0833    0.0769    0.0714
    0.0769    0.0714    0.0667
    0.0714    0.0667    0.0625
    0.0667    0.0625    0.0588
    0.0625    0.0588    0.0556
    0.0588    0.0556    0.0526

We can compare the exact solution (which we can obtain using the analytically known exact inverse)to the solutions obtained using backslash and inv. To compare the differences between solutions, we use the norm function norm(e,Inf), defined as $$ {\tt norm(e,Inf)} = \parallel \! u \! \parallel_\infty = \max_{i} |e_i| $$

x_true = invhilb(10)*b;
error_backslash = norm(A\b-x_true,Inf)/norm(x_true,Inf);
error_inv = norm(inv(A)*b - invhilb(10)*b,Inf)/norm(x_true,Inf);

Here are the results.

fprintf('Error using backslash : %12.4e\n',error_backslash);
fprintf('Error using inv       : %12.4e\n',error_inv);
Error using backslash :   9.1573e-05
Error using inv       :   5.3014e-02

The solution obtained from the backslash operator is several order of magnitude better than the the solution obtained using the inverse operator.

Using x = A^(-1)*b results in the same error as inv operator, and should also be avoided.

Back to the top

Lab exercises

Solve for the unknowns in the following problems.
  1. $$\begin{eqnarray} x - y & = & 5 \\ 3x + 5y & = & 1 \end{eqnarray} $$
  2. $$\begin{eqnarray} x - y + 3z & = & 5 \\ 3x + 5y - z & = & 1 \\ 5y - 3z & = & -22 \\ \end{eqnarray} $$
  3. $$\begin{eqnarray} a & = & 4(b - c) \\ \frac{b-a}{2c} & = & 1 \\ a + b & = & c-2 \end{eqnarray} $$
  4. Solve the linear system $$\begin{eqnarray} ax_1^2 + b x_1 + c & = & y_1 \\ ax_2^2 + b x_2 + c & = & y_2 \\ ax_3^2 + b x_3 + c & = & y_3 \end{eqnarray} $$ for unknown parameters $a$, $b$ and $c$, where points $(x_1,y_1)$, $(x_2,y_2)$, $(x_3,y_3)$ are known and are given by $(1,3)$, $(2,5)$, $(3,-1)$, respectively.
  5. Solve the linear system $A{\bf x} = {\bf b}$, where $A$ is a $5 \times 5$ matrix with ones on the sub- and super-diagonal, and -2 on the diagonal, i.e. $$ A = \left[\begin{array}{rrrrr} -2 & 1 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 1 & -2 \\ \end{array}\right]. $$ The right hand side vector ${\bf b}$ is a $5 \times 1$ column vector of ones. Create the system using Matlab functions (1) diag and (2) spdiag.
  6. Often, linear systems arise in the context of solving recurrence relations. For example, the famous Fibonacci sequence, named after the 13th century Italian mathematician Fibonacci, is generated from the recurrence relation $$ x_{j+1} = x_j + x_{j-1}, \qquad j = 1,2,3,4,... $$ where $x_0 = 0$ and $x_1 = 1$. The first few numbers in the sequence are $[0,1,1,2,3,5,8,13,21,34,...]$. Construct a linear system that relates values $x_j$, $j = 1,2,3,4,5...$ to the known starting values $x_0$ and $x_1$. Solve the linear system to generate the first 10 Fibonacci numbers. Try this using both diag and spdiag.

    Hint: Write out the equations for enough values of $j$ so that you can see a pattern in the linear system. Then construct the system.
  7. Suppose you only knew $x_0$ and $x_{11}$. Could you still generate the Fibonacci sequence? Modify the system you found above to solve for the first 10 values of the sequence, using only the known values $x_0 = 0$ and $x_{11} = 89$.

Compare your answers with the solutions.

Back to the top

Get the code

Do you want to try the above code fragments on your own? Download the Matlab script that produces this page here. (lab_7.m)

Powered by MathJax

Published with MATLAB® 8.5