Solving triangular systems

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

using FundamentalsNumericalComputation
A = rand(1.:9.,5,5)
L = tril(A)
5×5 Array{Float64,2}:
 8.0  0.0  0.0  0.0  0.0
 4.0  4.0  0.0  0.0  0.0
 7.0  3.0  5.0  0.0  0.0
 5.0  8.0  2.0  3.0  0.0
 4.0  8.0  2.0  7.0  8.0

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

b = ones(5)
x = FNC.forwardsub(L,b)
5-element Array{Float64,1}:
  0.125
  0.125
 -0.05
 -0.17499999999999996
  0.10312499999999997

It’s not clear what the error in this answer is. However, the residual, while not zero, is comparable to \(\varepsilon_\text{mach}\) in size.

b - L*x
5-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0

Next we’ll engineer a problem to which we know the exact answer.

alpha = 0.3;
beta = 2.2;
U = diagm(0=>ones(5),1=>[-1,-1,-1,-1])
U[1,[4,5]] = [ alpha-beta, beta ]
U
5×5 Array{Float64,2}:
 1.0  -1.0   0.0  -1.9   2.2
 0.0   1.0  -1.0   0.0   0.0
 0.0   0.0   1.0  -1.0   0.0
 0.0   0.0   0.0   1.0  -1.0
 0.0   0.0   0.0   0.0   1.0
x_exact = ones(5)
b = [alpha,0,0,0,1]

x = FNC.backsub(U,b)
err = x - x_exact
5-element Array{Float64,1}:
 2.220446049250313e-16
 0.0
 0.0
 0.0
 0.0

Everything seems OK here. But another example, with a different value for \(\beta\), is more troubling.

alpha = 0.3;
beta = 1e12;
U = diagm(0=>ones(5),1=>[-1,-1,-1,-1])
U[1,[4,5]] = [ alpha-beta, beta ]
b = [alpha,0,0,0,1]

x = FNC.backsub(U,b)
err = x - x_exact
5-element Array{Float64,1}:
 -4.882812499995559e-5
  0.0
  0.0
  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 \(x_1\) performs \((\alpha-\beta)+\beta\) in the first row. Since \(|\alpha|\) is so much smaller than \(|\beta|\), this a recipe for losing digits to subtractive cancellation.