Suppose that we are solving a system of *n* linear equations
**U x** = **b** where **U** = (u_{ij}) is an upper-triangular
matrix with no zero entries on the diagonal. The steps in
solving for **x** = (*x*_{i}) are:

For *i* = *n*, *n* - 1, ..., 2, 1, in that order, let:

This is called backward substitution because we are starting
with the last unknown *x*_{n} and, having solved for
it, using it to solve for the previous unknown *x*_{n - 1}, and so on.
In Matlab, this may be done by the (poorly implemented) code:

n = length( b );
x = zeros( n, 1 );
for i=n:-1:1
x(i) = b(i);
for j=(i + 1):n
x(i) = x(i) - U(i, j)*x(j);
end
x(i) = x(i)/U(i, i);
end

Note that if we initialize the vector **x** to be the zero vector,
then the sum is simply the dot product of the *i*th row of A
dotted with *x*. Thus, we may code this algorithm in Matlab
as:

n = length( b );
x = zeros( n, 1 );
for i=n:-1:1
x(i) = ( b(i) - U(i, :)*x )/U(i, i);
end

Note: this is **significantly faster** than the previous implementation
(and is also easier to remember).

Copyright ©2005 by Douglas Wilhelm Harder. All rights reserved.