Chapter 2

f_02_01.jl

Gaussian elimination for linear simultaneous equations

Other versions: j_02_01.jl, nm21.f95



In [16]:

    
using NMfE

Gaussian elimination of an equation like:

$[\mathbf{A}]\{\mathbf{x}\}=\{\mathbf{b}\}$

First define matrix A and right hand side vector b:



In [17]:

    
A=[10. 1. -5.;-20. 3. 20.;5. 3. 5.]









    Out[17]:





3x3 Array{Float64,2}:
  10.0  1.0  -5.0
 -20.0  3.0  20.0
   5.0  3.0   5.0



In [18]:

    
b=[1., 2., 6.]









    Out[18]:





3-element Array{Float64,1}:
 1.0
 2.0
 6.0

The function lufac() can be found in the src directory of package NMfE. lufac() is a pretty 'direct' translation of lufac.f95 to Julia. Julia has a built-in function called lufact() which will be used below.



In [19]:

    
n = size(A, 1)
x = 0.0









    Out[19]:





0.0



In [20]:

    
# Convert to upper triang;e
for k in 1:n-1
  if abs(A[k, k]) > 1.0e-6
    for i in k+1:n
      x = A[i, k] / A[k, k]
      A[i, k] = 0.0
      for j in k+1:n
        A[i, j] -= A[k, j] * x
      end
      b[i] -= b[k] * x
    end
  else
    println("Zero pivot found in row $k")
  end
end



In [21]:

    
A |> display









    





3x3 Array{Float64,2}:
 10.0  1.0  -5.0
  0.0  5.0  10.0
  0.0  0.0   2.5



In [22]:

    
for i in n:-1:1
  x = b[i]
  if i < n
    for j in i+1:n
      x -= A[i, j] * b[j]
    end
  end
  b[i] = x / A[i, i]
end



In [23]:

    
println("Solution vector: $b")









    



Solution vector: [1.0,-2.0,1.4]



In [24]:

    
# Copy b to d for later use in @assert test. Notice the use of copy!
d = copy(b)









    Out[24]:





3-element Array{Float64,1}:
  1.0
 -2.0
  1.4



In [30]:

    
# Restore A and b
A = [10. 1. -5.; -20. 3. 20.; 5. 3. 5.]
b = [1., 2., 6.];



In [26]:

    
c = A\b









    Out[26]:





3-element Array{Float64,1}:
  1.0
 -2.0
  1.4



In [27]:

    
@assert round(A * c, 14) == b "A * c !== b"



In [28]:

    
@assert round(c, 14) == d "c !== d"



In [29]:

    
?@assert









    Out[29]:





@assert cond [text]
Throw an AssertionError if cond is false. Preferred syntax for writing assertions. Message text is optionally displayed upon assertion failure.