Using LU factorization

using FundamentalsNumericalComputation
A = [2 0 4 3; -4 5 -7 -10; 1 15 2 -4.5; -2 0 2 -13];
L,U = FNC.lufact(A);
L
4×4 Array{Float64,2}:
  1.0  0.0   0.0  0.0
 -2.0  1.0   0.0  0.0
  0.5  3.0   1.0  0.0
 -1.0  0.0  -2.0  1.0
U
4×4 Array{Float64,2}:
 2.0  0.0   4.0   3.0
 0.0  5.0   1.0  -4.0
 0.0  0.0  -3.0   6.0
 0.0  0.0   0.0   2.0
LtimesU = L*U
4×4 Array{Float64,2}:
  2.0   0.0   4.0    3.0
 -4.0   5.0  -7.0  -10.0
  1.0  15.0   2.0   -4.5
 -2.0   0.0   2.0  -13.0

It’s best to compare two floating-point quantities by taking their difference.

A - LtimesU
4×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

(Usually we can expect “zero” only up to machine epsilon at best. However, all the exact numbers in this example are also floating-point numbers.)

To solve a linear system, we no longer need the matrix \(A\).

b = [4,9,29,40]
z = FNC.forwardsub(L,b)
x = FNC.backsub(U,z)
4-element Array{Float64,1}:
 -3.0
  1.0
  4.0
 -2.0
b - A*x
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0