实现高斯赛德尔的矩阵项版本

Question

我正在尝试实现第 1 章中的以下描述。 Heath 的“科学计算入门概述”中的第 11 部分，介绍了求解线性方程组的高斯-赛德尔迭代法

但是，我一定遗漏了一些东西，因为完全按照原样编写方程无法产生正确的结果。请帮忙！

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]

Answer 1

根据 @Ian Bush 的评论，我误解了方程中的 L、U 和 D。 L 和 U 是 A 的严格下三角矩阵和上三角矩阵，而 D 是 A 的对角线项。这是与所描述的方程相匹配的校正函数（即使在实践中可能无法计算逆 D + L）。

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)                                                      
    U[diagind(U)] .= 0                                                          
                                                                                
    # 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]