Solving nonsingular 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 nonsingular 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, nonsingular systems. In a later tutorial, we will
investigate the solution to overdetermined and underdetermined 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.4409e16 8.8818e16
Back to the top
The LU Decomposition
Behind the scenes, what is the "backslash" operator doing?
If $A$
is a square, nonsingular 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 nonzero 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 illconditioned 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\bx_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.1573e05 Error using inv : 5.3014e02
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
 $$\begin{eqnarray} x  y & = & 5 \\ 3x + 5y & = & 1 \end{eqnarray} $$
 $$\begin{eqnarray} x  y + 3z & = & 5 \\ 3x + 5y  z & = & 1 \\ 5y  3z & = & 22 \\ \end{eqnarray} $$
 $$\begin{eqnarray} a & = & 4(b  c) \\ \frac{ba}{2c} & = & 1 \\ a + b & = & c2 \end{eqnarray} $$
 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.
 Solve the linear system $A{\bf x} = {\bf b}$, where $A$ is a $5 \times 5$ matrix with ones on the sub and superdiagonal, 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.

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_{j1}, \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.  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)
Published with MATLAB® 8.5