Poisson 方程的有限差分矩阵是对称的和正定的。因此，预处理共轭梯度算法是该问题的首选迭代求解器。

预处理器的选择对方法的收敛性有很大的影响。众所周知，不完全 Cholesky 分解可以很好地解决这个问题。

以下 MATLAB 代码构造 3D Poisson 问题的有限差分矩阵，并求解所有 1 右侧的方程。

function laplaceEqnTestCSE
n = 50;
L=1;
h = L/(n+1);
dim=3;
K=laplaceEqn(dim, n);
neq = rows(K);
b = ones(neq,1);
ne2 = ceil(neq/2);
tic
disp('Begin solve.');
if(0)
u = K\b;
else
tol=1e-8;
maxIter=100;
L = ichol(K);
u = pcg(K, b, tol, maxIter, L, L');
end
toc
end

function K=laplaceEqn(dim, n)
h = 1/(n+1);

K1D = spdiags(ones(n,1)*[-1 2 -1],-1:1,n,n);  % 1d Poisson matrix
%subplot(2,3,4), spy(K1D)
if(dim==1)
K = K1D/h^2;
return;
end

I1D = speye(size(K1D));                       % 1d identity matrix
K2D = kron(K1D,I1D)+kron(I1D,K1D);            % 2d Poisson matrix
%subplot(2,3,5), spy(K2D)
if(dim==2)
K = K2D/h^2;
return;
end

I2D = speye(size(K2D));                       % 2d identity matrix
K3D = kron(K2D,I1D)+kron(I2D,K1D);            % 3d Poisson matrix
if(dim==3)
K = K3D/h^2;
return;
end

end

我在 Octave 中运行了这段代码，并在 54 次迭代中得到了一个收敛的解决方案。在我的桌面 Windows PC 上花了大约一秒钟。

一些评论指出，一些离散边界条件或引入空间变化系数的方法可能会导致不对称的有限差分矩阵。但从根本上说，拉普拉斯算子是自伴的，因此理想情况下，我们希望数值系数矩阵也是自伴的（即对称的）。如果不是这种情况，离散化方法可能不是最好的方法。Leveque 的有限差分法入门书（Leveque) 以前在这里推荐过。它解决了如何应用边界条件的问题，以及在示例 2.1 中如何处理空间变化的系数以保持系数矩阵的对称性。套用 Leveque 的话说，通过取导数然后用 FD 对其进行近似来处理空间变化系数的最“明显”方法并不是离散这种情况的最佳方法。他更详细地讨论了为什么原始 PDE 的重要基本性质应该反映在数值近似中。