Simple 1D and 2D Poisson solvers in Octave and Fortran

Hi there,
here are some MATLAB/Octave and Fortran90 scripts for solving the Poisson's equation. First, we will show the finite-difference formulation in 1D and 2D, then we discuss hot to use the FISHPACK library for solving 2D problems.

In 1D, using central second order finite differences, we write:

We use V for the unknown here. In case of the electric potential (V standing for "voltage"), the right hand side b will be equal to the electric charge "rho" divided by the vacuum constant "epsilon0". The discretized equation creates a linear system that is roughly tri-diagonal (depending on the boundary conditions), as will be discussed in the following.

Preliminary note on the numerical solution of linear systems

Dirichlet BCs: imposed V at the extremes

If we want to impose the potential V at the boundaries (Dirichlet Boundary Conditions), we use the notation V(0) = V0 and V(L) = V_N and we compose a linear system that happens to be tridiagonal. Since the values at the extremes are given, there are two ways to write this system: either we put these values at the right hand side and solve only for the internal points (Option 1 in the following picture), or we treat them as unknowns and prescribe an equation in the form "1*V0 = V0" (Option 2 in the following).

Here, we provide some scripts to solve this problem:

• Octave script;
• Fortran 90 file.

Periodic BCs

Periodic boundary conditions are obtained through the condition "V(0) = V(L)": the solution at the two extremes of the domain is the same. This translates into the discretized condition "V1 = VN", that is implemented as shown in the following image. Notice that the most naive implementation that we may think of is actually wrong, as it gives rise to a singular matrix (see the image). The reason lies into the mathematics of the problem: periodic BCs are not sufficient to fully determine the solution. Indeed, specifying that "V(0) = V(L)" imposes only one condition to the system, but we know that a second order differential equation (such as the Poisson's equation) requires two BCs. To put it in other terms: every solution that differs by a constant still satisfies the differential equation and the periodic BCs. Numerically, this translates into a singular discretization matrix.

This problem is easily solved though, by fixing the potential at one of the two extremes. Most likely, this value will not matter physically, as we will be interested in the derivative of the potential. What value you set does not matter much, you may set zero for instance. Notice that you could cook up a formulation where you set the value of some other node of the domain (like the center for example, to impose some further symmetry).

Ok, so we set "VN = Vfix = 0". This will be the last row of our matrix, that will read "[0 0 .... 0 1] [V1 V2 .... VN]^T = Vfix". The first point on the other hand will have the condition "V1 = VN", that will introduce an off-diagonal entry in the matrix. But if you think about that, imposing VN also means imposing V1. So, ultimately, fixing the value of one extreme, also means fixing the value for the other one! The system now has two Dirichlet BCs, and can be written in diagonal form, as follows:

And here are the scripts (pretty much identical to the Dirichlet versions):

• Octave script;
• Fortran 90 file.

2D Poisson equation with Octave

We now move to the 2D Poisson equation. The idea is analogous as before: all one has to do is discretizing both the x and y directions.

The final result is quite a large matrix of size M(Nx*Ny, Nx*Ny). The matrix is very sparse, and you should remember this if you want a decent numerical efficiency. For this problem, it is convenient to introduce an additional point numbering, let's call it "ID", that goes from 1 to Nx*Ny. Moving from this "ID" numbering to "(i,j)" is of course very simple (see the following Octave script). This numbering is convenient because the solution and the source term will be shaped as a large column vector, ranging from ID = 1 to ID = Nx*Ny.

Now, every entry in the matrix represents the discretized derivative in one single node "ID". You can build it in the following way:

1. You loop on the rows, using the index ID = 1, Nx*Ny; ID is the point (i,j);
2. You locate the IDs corresponding to the four points (i+1,j), (i-1,j), (i,j+1), (i,j-1);
3. At the matrix row "ID", you assign a value "-4" to the column "ID" and a value "+1" at the column relative to the other 4 points (same row);

Boundary conditions are imposed as discussed before for the 1D cases. Here is a simple Octave implementation for this problem.

2D Poisson equation with Fortran: FISHPACK90

A Fortran implementation could be easily prepared following the structure of the Octave file. However, we employ here the FISHPACK library, that is much more flexible. (In case the link is broken, here's a copy of the sources - all credits to the authors [link: click here]) Indeed, the 2D problem is not trivial in some cases: consider for instance a problem that is periodic in one direction, but you want Dirichlet BCs on the other direction. That's not trivial, but FISHPACK does it easily and quite efficiently. The FISHPACK library actually solves the inhomogeneous Helmholtz equation, that is a generalization of the Poisson equation. You just have to set the parameter "LAMBDA" equal to zero to recover Poisson. We use here the "HWSCRT" cartesian (CRT) solver: you need to specify the following parameters:

• xmin, xmax, Nx-1, BC x type, dummy vector(Nx), dummy vector(Nx);
• ymin, ymax, Ny-1, BC y type, dummy vector(Ny), dummy vector(Ny);
• LAMBDA, b(i,j), IDIMF, PERTURB, IERR.

The vectors that I have denoted as "dummy" are needeed only if you want to specify neumann BCs. The matrix b(i,j) contains the source terms for internal points, and the BCs in case of Dirichlet conditions. IDIMF, PERTURB and IERR are just some parameters and you should check the documentation. In the following tar.gz file you will find a Fortran implementation for a full-Dirichlet case and for a test case with Dirichlet along x and periodic BCs along y.

CLICK HERE TO DOWNLOAD THE TAR.GZ FILE.

ATTENTION! MEMORY LEAKS!

The FISHPACK90 library is known to suffer from memory leaks when compiled with gfortran. This happens to be a problem when you call a solver multiple times. Imagine that you are running some simulation where time evolves, and imagine that at every single time step you need to solve the Poisson equation. Well, memory leaks will accumulate, and you will see the system memory filling up, until your computer will start filling up the swap memory and everything will get so slow, until eventually the OS dies.

After some checks using valgrind, fishpack appears not to free some allocated memory indeed: to fix things in a dirty way you can take the solver that you are interested in and freed the memory manually:

      [...]
!     allocate required work space length (generous estimate)
      IRWK=4*(N+1)+(13+INT(ALOG(FLOAT(N+1)/ALOG(2.0))))*(M+1)
      ICWK = 0
      CALL ALLOCATFISH (IRWK, ICWK, W, IERROR)
!     check that allocation was successful
      IF (IERROR == 20) RETURN
      call hwscrtt(A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC,BDD,
     +             ELMBDA,F,IDIMF,PERTRB,IERROR,w%rew)

!     release allocated work space
      CALL FISHFIN (W)
 
! FISHFIN WILL NOT DO ITS JOB... JUST DEALLOCATE HERE
      deallocate(W%rew)

      RETURN
      END SUBROUTINE HWSCRT