julia

Implementing matrix term version of Gauss-seidel


I am trying to implement the below description from Ch. 11 of Heath's "Scientific Computing An Introductory Survey" of the Gauss-Seidel iterative method for solving a system of linear equations

enter image description here

however, I must be missing something because writing the equation exactly as is fails to produce the correct results. Please help!

using LinearAlgebra: lu
# From dicretization of Laplace equation
A = [
  4.0  -1.0  -1.0   0.0
 -1.0   4.0   0.0  -1.0
 -1.0   0.0   4.0  -1.0
  0.0  -1.0  -1.0   4.0]

# block diagonal matrix of A
D = [
  4.0  -1.0   0.0   0.0
 -1.0   4.0   0.0   0.0
  0.0   0.0   4.0  -1.0
  0.0   0.0  -1.0   4.0]

b = [0, 0, 1, 1] # rhs for laplace w/ Dirichlet boundary conditions
L, U = lu(A)

# gauss seidel using matrix terms
niters = 100
D_plus_L_inv = inv(D + L)
x = zeros(length(b))
for k in 1:niters
  x = D_plus_L_inv*(b - U*x)
end 

# clearly not equal
@show A \ b
# A \ b = [0.12499999999999999, 0.12499999999999999, #0.37499999999999994, 0.37499999999999994]

@show x
# x = [0.021795262674411238, 0.023610129063946845, 0.14893710589872386, 0.14211027447080438]

Solution

  • Per @Ian Bush comments, I was misinterpreting what L, U, and D were in the equations. L and U are strict lower and upper triangular matrices of A, while D is the diagonal entries of A. Here is the corrected function matching the described equation (even though computing the inverse D + L may not be done in practice).

    using LinearAlgebra
    
    """                                                                             
        gauss_seidel_matrix_form_solver(                                            
            A_::Matrix, b::Vector, x0::Vector, niters::Int)                          
                                                                                    
    Return solution to linear system `Ax = b` via matrix Gauss-Seidel method.       
                                                                                    
    # References                                                                    
    [1] : Ch. 11.5.3 Heath.                                                         
    """                                                                             
    function gauss_seidel_matrix_form_solver(                                       
        A_::Matrix, b::Vector, x0::Vector, niters::Int)  
        A = copy(A_)
                                
        # Get diagonal entries of A, store in matrix                                
        D_vec = diag(A)                                                             
        D = diagm(D_vec)                                                            
                                                                                    
        # Get the strict lower and upper triangular matrices                        
        L = LowerTriangular(A)                                                      
        L[diagind(L)] .= 0                                                          
        U = UpperTriangular(A)                                                                                                              
                                                                                    
        # Compute the inverse of D + L ahead of time                                
        D_plus_L_inv = inv(D + L)                                                   
                                                                                    
        # initialize the guess                                                      
        x = x0                                                                      
                                                                                    
        # gauss seidel iterations                                                   
        for k in 1:niters                                                           
            x = D_plus_L_inv*(b - U*x)                                              
        end                                                                         
                                                                                    
        return x                                                                    
    end  
    
    A = [
      4.0  -1.0  -1.0   0.0
     -1.0   4.0   0.0  -1.0
     -1.0   0.0   4.0  -1.0
      0.0  -1.0  -1.0   4.0]
    
    b = [0, 0, 1, 1]
    
    u_gs = gauss_seidel_matrix_form_solver(A, b, zeros(length(b)), 100)
    u_builtin = A \ b
    
    @show u_gs
    # u_gs = [0.125, 0.125, 0.375, 0.375]
    @show u_builtin
    # u_builtin = [0.12499999999999999, 0.12499999999999999, 0.37499999999999994, 0.37499999999999994]