Table of contents

1. Poisson’s equation
2. Solving Poisson’s equation numerically
   1. Visualization of iterative solutions
3. Class material
   1. Additional resources
   2. Footnotes

Poisson’s equation

The electrostatic potential \(\Phi(\mathbf{r})\) is related to the charge density \(\rho(\mathbf{r})\) through Poisson’s equation

\[\boldsymbol{\nabla}^2 \Phi = -4\pi\rho\]

a hyperbolic PDE. In Cartesian coordinates it reads

\[\frac{\partial^2 \Phi(x, y, z)}{\partial x^2} + \frac{\partial^2 \Phi(x, y, z)}{\partial y^2} + \frac{\partial^2 \Phi(x, y, z)}{\partial z^2} = -4 \pi \rho(x,y,z).\]

If no charges are inside the domain boundaries, Poisson’s equation reduces to Laplace’s equation

\[\frac{\partial^2 \Phi}{\partial x^2} + \frac{\partial^2 \Phi}{\partial y^2} + \frac{\partial^2 \Phi}{\partial z^2} = 0.\]

We can solve either equation if we are either given the value of the potential on the boundary or the electric field on the boundary normal to the surface.

Solving Poisson’s equation numerically

We will look at various finite difference schemes (both explicit and implicit) to solve the different classes of PDEs on discretized grids in space and time. For Poisson’s equation, we will introduce a simple iterative scheme to successively improve an initial (poor) guess for the solution. The most simple of these iterative algorithms is the Jacobi method that updates the value of the potential (here in 2D) at grid point \((i, j)\) as the average of the neighbors,

\[\Phi_{i,j} = \frac{1}{4}\Big(\Phi_{i+1,j} + \Phi_{i-1,j} + \Phi_{i,j+1} + \Phi_{i,j-1}\Big) + \pi\rho_{i,j} \Delta^2\]

where \(\Delta\) is the spacing between grid points in the \(x\) and \(y\) direction.

We start with \(\Phi_{i,j} = 0\) inside the domain boundaries and the potential defined on the boundary itself.

In Python, the Jacobi update looks very similar

Phi[1:-1, 1:-1] = 0.25*(Phi[2:, 1:-1] + Phi[0:-2, 1:-1] + \
                        Phi[1:-1, 2:] + Phi[1:-1, 0:-2]) + \
                        np.pi * rho[1:-1, 1:-1] * \Delta**2

where we made use of NumPy slicing and element-wise operations (as explained in numpy_arrays.ipynb).

Visualization of iterative solutions

The movies show how the solutions of simple 2D electrostatic problems relax to the converged solution (using successive over-relaxation with the Gauss-Seidel algorithm), namely a box with one side at 100 V and three sides at 0 V (left) and the same box with an additional charge dipole included (right).

Class material

1. 14_PDEs-1.ipynb: Poisson’s and Laplace’s equation: Background (board notes on PDEs (PDF)), Gauss-Seidel algorithm, problems: 14_PDEs-1-Students.ipynb ¹
2. 14_PDEs-2.ipynb ¹: fast Jacobi, Gauss-Seidel and Successive Over Relaxation algorithms (using NumPy array operations, see board notes on numpyfied Poisson solvers (PDF)) and how to replace slow Python loops with fastnumpy array operations (numpy_arrays.ipynb), Poisson equation, charge density: 14_PDEs-2-Students.ipynb ¹
3. When plotting arrays of data with matplotlib, one sooner or later one needs to use the numpy.meshgrid() function which “returns coordinate matrices from coordinate vectors”—most users find this function mystifying. The notebook meshgrid.ipynb makes a bit clearer, what meshgrid does and why it is useful.

Additional resources

• Computational Physics, Ch 19–23
• Partial differential equation, Andrei D. Polyanin et al. (2008), Scholarpedia, 3(10):4605. doi: 10.4249/scholarpedia.4605
• Numerical Recipes in C, WH Press, SA Teukolsky, WT Vetterling, BP Flannery. 2nd ed, 2002. Cambridge University Press. Chapter 19.

---

Footnotes

1. As usual, git pull the resources repository PHY432-resources to get a local copy of the notebook. Then copy the notebook and all other code into your work directory in order to complete the exercises. ↩ ↩² ↩³