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.