# Iteratively solving 3D Poisson equation in MATLAB

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!

• You need to provide more information about how you are using the matlab iterative solvers. Are you using pcg? With or without preconditioning from, e.g. ichol? – Bill Greene Feb 14 '17 at 19:13
• For what $Nx, Ny, Nz$ do you get this error message? If you got this error for small dimensions, then it is something wrong in the calls to matlab routines. As for preconditioner, one of the choices probably can be a discrete laplacian with Dirichlet boundary conditions ignoring the structure inside the domain. For cubic domains this preconditioner can be easily implemented using discrete Fourier transform. – VorKir Feb 14 '17 at 23:37
• @BillGreene like I mentioned, my knowledge on solving such systems is (at this moment) very limited. So when I was referring to use built-in iterative solvers "out-of-the-box", I literally meant running e.g. bicgstab(A,b) (whose documentation merely says The n-by-n coefficient matrix A must be square and should be large and sparse. The column vector b must have length n.) or pcg(A,b) without any preconditioning. I did try to use ichol, but there I am stuck with Encountered nonpositive pivot.... – slvrbld Feb 15 '17 at 10:12
• @VorKir I have 2 particular test cases which I have successfully run with a direct solver and which I can't get to run with any iterative solver other than a (way too slow) SOR routine. Those cases are 80x50x100 (where A is sparse and real with ~2.6e6 nonzero terms) and a more recent case that also takes into account conductivity and permittivity where I have a 201x46x46 domain (where A is sparse and complex with ~2.8e6 nonzero terms). Those cases are pretty much the heaviest I can run with a direct solver. – slvrbld Feb 15 '17 at 10:22
• @slvrbld There is (to my knowledge) no convergence theory (and hence no convergence guarantee) for BiCGStab -- if it works, fine, if it doesn't, there's not much to go on. Sounds like you encountered the second case. If your matrix is not symmetric (and you can't make it so by changing the way the boundary conditions are incorporated -- which should be possible), GMRES or QMR would be the better option. – Christian Clason Feb 15 '17 at 15:08

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.

• This works as long as K is SPD, but my symmetry is broken at the moment I impose mixed boundary conditions to my nodes. ichol then refuses to factorize and pcg can (according to its description) no longer be used. For my Dirichlet BCs, I set the contribution of neighboring nodes to 0 and fix the value in my b vector. For Neumann conditions, I only consider the point itself and its neighbor perpendicular to the face I apply the condition on. Do you have a suggestion how to retain my symmetry and still be able to apply my BCs? – slvrbld Feb 15 '17 at 14:44
• If your approach to applying BCs doesn't maintain the symmetry, you can't use the conjugate gradient method but instead must use a different Krylov solver such as gmres. You might try incomplete LU factorization (ilu) as a preconditioner. – Bill Greene Feb 15 '17 at 15:57
• Thanks. Now we're getting somewhere! I should've mentioned explicitly that my matrix was non-symmetric. For my test-case, I've achieved convergence to acceptable residuals in times very comparable to the direct approach (~120s for my 201x46x46 case with complex A) by using ilu with default properties for preconditioning. When using gmres, I noticed that the rate of convergence is very sensitive to the RESTART parameter. Could you modify your answer (if possible with some info on how to optimize my ilu and gmres settings) to better reflect the question? Then I'll mark it accepted. – slvrbld Feb 16 '17 at 9:30
• For questions on iterative solvers, in general, and gmres, in particular, I highly recommend Saad's book on this topic which he has conveniently made available online: www-users.cs.umn.edu/~saad/books.html. In my previous comment I should have also mentioned that I believe it is possible to apply both Dirichlet and Neumann constraints in a way that the matrix will be symmetric. This post might be of some help in that regard: scicomp.stackexchange.com/questions/5355. IMHO, this is the right approach even if it requires a bit more work. – Bill Greene Feb 16 '17 at 15:45
• Thanks for the link to Saad's book, seems like a very useful resource for me. Christian Clason's comment to my original post actually made me look into ways to get my BCs in by modifying the RHS as well, but meanwhile I have also been extending the model to include dielectrics and conductors. Unfortunately, retaining symmetry when I also include non-constant relative permittivity seems not possible (please correct me if I'm wrong here :)). However, it will be a nice exercise for me to get my old model to produce symmetric matrices. – slvrbld Feb 16 '17 at 16:59