Gaussian elimination WITHOUT pivoting

You need to download the m-file lunp.m and put it in the same folder as your other m-files.

Solving a linear system using Gaussian elimination WITHOUT pivoting

We want to solve the linear system

clearvars                         % Always clear variables at the beginning
format compact                    % omit blank lines between results
format short g                    % show each entry with 5 digits accuracy
A = [4 3 3; 8 7 9; -8 -1 3]
A = 3×3
     4     3     3
     8     7     9
    -8    -1     3
b = [1;4;2]
b = 3×1
     1
     4
     2

We use Gaussian elimination without pivoting:

elimination in column 1: For

: subtract

times row 1 from row j

elimination in column 2: For

: subtract

times row 2 from row j

...

This transforms the matrix A to an upper triangular matrix U.

The lower triangular matrix L has 1's on the diagonal and the multipliers

below the diagonal.

[L,U] = lunp(A)
L = 3×3
     1     0     0
     2     1     0
    -2     5     1
U = 3×3
     4     3     3
     0     1     3
     0     0    -6

For a given right hand side vector b:

Solve

by forward substitution: (Matlab's \ command checks whether the matrix is triangular and then performs substitution in the correct order)

y = L\b
y = 3×1
     1
     2
    -6

Solve

by back substitution to obtain the solution of the linear system

Problems with this algorithm:

(1) Algorithm can fail even for nonsingular matrix A when zero pivot is encountered:

We want to solve the linear system

. The algorithm stops since it encounters the pivot

A = [0 1; 2 0]
A = 2×2
     0     1
     2     0
[L,U] = lunp(A)
Warning: LU-decomposition does not exist
L =
     []
U =
     []

(2) Algorithm can give very bad results due to tiny pivots:

We want to solve the linear system

A = [.1 .3 0; .3 .9 1; 0 1 1]
A = 3×3
          0.1          0.3            0
          0.3          0.9            1
            0            1            1
b = [1;1;1]
b = 3×1
     1
     1
     1

The exact solution is

xex = A\b
xex = 3×1
            1
            3
           -2

We now use Gaussian elimination without pivoting:

[L,U] = lunp(A)
L = 3×3
            1            0            0
            3            1            0
            0   4.5036e+15            1
U = 3×3
          0.1          0.3            0
            0   2.2204e-16            1
            0            0  -4.5036e+15

Note that we obtain a tiny pivot element

. In exact arithmetic we would have obtained a zero pivot.

y = L\b
y = 3×1
            1
           -2
   9.0072e+15
xhat = U\y
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  4.930381e-32.
xhat = 3×1
           -2
            4
           -2

There is no division by zero, but Matlab alerts us that there may be loss of accuracy.

Note that the computed solution

is completely wrong: 100% relative error for x, and 58% relative error for

relerr_x = norm(xhat-xex,Inf)/norm(xex,Inf)  % relative error of xhat
relerr_x =             1
bhat = A*xhat;
relerr_b = norm(bhat-b,Inf)/norm(b)          % backward error: relative error of bhat = A*xhat
relerr_b =       0.57735