实现Gauss-Seidel迭代法矩阵形式结果异常，请求技术帮助

我正尝试参考Heath所著《Scientific Computing: An Introductory Survey》第11章中求解线性方程组的Gauss-Seidel迭代法描述，实现其矩阵形式，但严格按照公式编写代码后无法得到正确结果。以下是我的代码及运行结果，恳请协助排查问题：

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]

错误原因分析 #

1. D矩阵定义错误：Gauss-Seidel迭代中的D是原矩阵A的纯对角矩阵，仅保留对角线元素，而非你定义的块对角矩阵。当前D包含了下三角的非对角元素，完全不符合迭代公式要求。
2. L、U矩阵拆分错误：你用lu(A)得到的L和U是LU分解的结果（L为单位下三角，U为上三角），但Gauss-Seidel迭代要求：
   • L：A的严格下三角矩阵（对角线元素为0，仅保留对角线下方元素）
   • U：A的严格上三角矩阵（对角线元素为0，仅保留对角线上方元素）

修正后的代码 #

using LinearAlgebra

# 原矩阵与右端项
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, 0.0, 1.0, 1.0]

# 正确拆分A为L（严格下三角） + D（对角） + U（严格上三角）
D = Diagonal(A)          # 提取A的对角元素构造对角矩阵
L = tril(A, -1)          # 严格下三角矩阵
U = triu(A, 1)           # 严格上三角矩阵

# Gauss-Seidel迭代
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

# 对比结果
@show A \ b
@show x

运行结果 #

修正后迭代100次的结果会与A \ b的精确解完全一致：

A \ b = [0.12499999999999999, 0.12499999999999999, 0.37499999999999994, 0.37499999999999994]
x = [0.125, 0.125, 0.375, 0.375]

内容的提问来源于stack exchange，提问作者Jared