Iteratively solving 3D Poisson equation in MATLAB

Question

I have written a function that sets up a sparse matrix A and RHS b for the 3D Poisson equation in a relatively efficient way. The set-up is nothing fancy: I have extended the 2D 5-point stencil to an equivalent 7-point stencil for 3D. Conductors are (at this moment) simply blocks of Dirichlet BCs and I am not (yet) taking dielectrics into account. Therefore, the constructed A is quite straightforward:

• Maximum number of non-zero terms is Nx*Ny*Nz*7
• In practice, my actual number of non-zero terms is about 90% of the maximum
• Due to the BCs, A is non-symmetric (I should have mentioned this from the start)

The correctness of my arrays has been verified by running some very simple problems with Nx, Ny and Nz between 50-100 and then using mldivide. However, my peak memory usage is already hitting ~20GB for these test problems, while I'll be using significantly larger arrays for my actual problem.

So now I'm looking into iterative approaches, but my knowledge currently does not extend beyond Gauss-Seidel, Jacobi and SOR. When I attempt to use any of the built-in iterative solvers out of the box, I typically end up with the following message:

METHOD stopped at iteration 2 without converging to the desired tolerance 1e-06
because a scalar quantity became too small or too large to continue computing.
The iterate returned (number 1) has relative residual 0.97.

As this is quite a common system to solve, I was expecting to easily find some resources on this, but I am finding significantly less useful information than I expected. Some advice and pointers would be appreciated!

Answer 1

The finite difference matrix for the Poisson equation is symmetric and positive definite. So the preconditioned conjugate gradient algorithm is the iterative solver of choice for this problem.

The choice of preconditioner has a big effect on the convergence of the method. Incomplete Cholesky factorization is known to work well for this problem.

The following MATLAB code constructs the finite difference matrix for the 3D Poisson problem and solves the equation for a right-hand-side of all ones.

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

I ran this code in Octave and got a converged solution in 54 iterations. It took about a second on my desktop Windows PC.

Several of the comments pointed out that some methods for discretizing the boundary conditions or the introduction of a spatially-varying coefficient can result in an unsymmetric finite difference matrix. But fundamentally, the Laplace operator is self-adjoint so, ideally, we would expect the numerical coefficient matrix to also be self-adjoint (i.e. symmetric). If that is not the case, the discretization approach is probably not the best one. Leveque's introductory book on finite difference methods (Leveque) has been previously recommended here. It addresses both the issue of how to apply boundary conditions and, in example 2.1, how to handle spatially-varying coefficients to maintain symmetry of the coefficient matrix. To paraphrase Leveque, the most "obvious" way to deal with spatially-varying coefficients by taking the derivative and then approximating it with FD is not the best way to discretize this case. He discusses in more detail why the important fundamental properties of the original PDE should be reflected in the numerical approximation.