The source is here.

Points of interest

use of iterative equation solver
convolution with the Green's function
use of diffuse domain approximation to represent interior surfaces
line and volume plots
imposition of exterior Neumann boundary conditions

Background

We wish to solve Laplaces's equation within a non-homogeneous medium, i.e., in which the permittivity is not constant, with an external electric field.

$\nabla . \left( \epsilon(r) \nabla v(r) \right) = 0$

Expanding and rearranging yields,

$\nabla^2 v(r) = - \frac{1}{\epsilon} \nabla \epsilon(r) . \nabla v(r)$

We can interpret $\nabla \epsilon . \nabla u / \epsilon$ as the induced surface charge density.

MADNESS's default Coulomb Green's function ( ) corresponds to the free-space Dirichlet boundary condition $v(\infty)=0$ . But we wish to impose the Neumann condition $- \partial v / \partial z(\infty) = E$ . Thus, we write our solution as the sum of the asymptotic solution ( ) and a component that goes to zero at infinity so that our MADNESS Green's function can be used. Substituting

$v(r) = - E z + u(r)$

our equation becomes

$\nabla^2 u(r) = \frac{1}{\epsilon} \nabla \epsilon(r) . \left(E - \nabla u(r) \right)$

The quantity in parentheses is the total electric field — the sum of the applied field and that due to the induced surface charge.

The exact solution is readily derived by expanding in spherical harmonics and matching the potential and electric displacement at the surface, or by looking inside Jackson, or by a quick google to find this page.

(Note that while the Green's function to Poisson's equation is $-1 / 4 \pi |r-s|$ , MADNESS provides convolution with , so you have to keep track of the $-1/4\pi$ yourself.)

(Note that to obtain an accurate solution with $\sigma<0.2$ you will need to increase the order of the wavelet and decrease the threshold, e.g., $k=8, \mbox{thresh}=1e-5$ .