Untitled

It's easy to get just the lower triangular part of any matrix using the tril command.

A = magic(5)
A = 
  24     1     8    15
   5     7    14    16
   6    13    20    22
  12    19    21     3
  18    25     2     9
L = tril(A)
L = 
   0     0     0     0
   5     0     0     0
   6    13     0     0
  12    19    21     0
  18    25     2     9

We'll set up and solve a linear system with this matrix.

b = ones(5,1);
x = forwardsub(L,b)
x = 
   5.8824e-02
  -7.0588e-02
   9.1403e-02
  -2.2754e-02
  -6.8448e-02

It's not clear what the error in this answer is. However, the residual, while not zero, is virtually

in size.

b - L*x
ans = 
            0
            0
   1.1102e-16
  -2.2204e-16
   4.4409e-16

Next we'll engineer a problem to which we know the exact answer. You should be able to convince yourself that for any α and β,

(missing equation)

alpha = 0.3;
beta = 2.2;
U = eye(5) + diag([-1 -1 -1 -1],1);
U(1,[4 5]) = [ alpha-beta, beta ];
x_exact = ones(5,1);
b = [alpha;0;0;0;1];

x = backsub(U,b);
err = x - x_exact
err = 
     0
     0
     0
     0
     0

Everything seems OK here. But another example, with a different value for β, is more troubling.

alpha = 0.3;
beta = 1e12;
U = eye(5) + diag([-1 -1 -1 -1],1);
U(1,[4 5]) = [ alpha-beta, beta ];
b = [alpha;0;0;0;1];
x = backsub(U,b);
err = x - x_exact
err = 
  -4.8828e-05
            0
            0
            0
            0

It's not so good to get four digits of accuracy after starting with sixteen! But the source of the error is not hard to track down. Solving for

performs

in the first row. Since

is so much smaller than

, this a recipe for losing digits to subtractive cancellation.