# A problem with Poisson equation

I'm computing the Hartree potentials of atoms by solving the Poisson equation and I use hydrogen atom as a test case. The Poisson equation for hydrogen atom in atomic units is given by $$\nabla^2 V_H = -4 \exp(-2 r)$$ where $$r = \sqrt{x^2+y^2+z^2}$$. The numerical solution of $$V_H(x,y,z)$$ with $$z=0$$ is illustrated in the following figure:

The numerical solution is computed with the Conjugate Gradient method so that the Laplacian is computed with stencil $$\nabla^2 f(x,y,z) \approx \frac{a+b+c-6 f(x,y,z)}{h^2}$$ where $$a = f(x+h,y,z)+f(x-h,y,z),$$ $$b = f(x,y+h,z)+f(x,y-h,z),$$ and $$c = f(x,y,z+h)+f(x,y,z-h).$$

As the right-hand side of the Poisson equation is spherically symmetric the Poisson equation takes the form $$\frac{d^2 v_H}{dr^2} + \frac{2}{r} \frac{d v_H}{dr} = -4 \exp(-2 r)$$ where $$v_H$$ is the Hartree potential as a function of $$r$$. The solution of this equation is $$v_H(r) = - \frac{r + 1}{r} \exp(-2r),$$ which is illustrated in the following figure: (here $$r=\sqrt{x^2+y^2+z^2}$$ and $$z=0$$).

The Octave code computing the Laplacian is here:

function aLap = StencilLap3d( a, rStep )
rH2 = rStep * rStep;
vSize = size( a );
nXSize = vSize( 1 );
nYSize = vSize( 2 );
nZSize = vSize( 3 );
nXDim = ( nXSize - 1 ) / 2;
nYDim = ( nYSize - 1 ) / 2;
nZDim = ( nZSize - 1 ) / 2;
assert( nXDim == round( nXDim ) );
assert( nYDim == round( nYDim ) );
assert( nZDim == round( nZDim ) );
aLap = zeros( nXSize, nYSize, nZSize );
rXP = 0;
rXM = 0;
rYP = 0;
rYM = 0;
rZP = 0;
rZM = 0;
for nX = (-nXDim):nXDim
for nY = (-nYDim):nYDim
for nZ = (-nZDim):nZDim
if ( nX < nXDim )
rXP = a( nXDim + 1 + nX + 1, nYDim + 1 + nY, nZDim + 1 + nZ );
else
rXP = 0.0;
endif
if ( nX > -nXDim )
rXM = a( nXDim + 1 + nX - 1, nYDim + 1 + nY, nZDim + 1 + nZ );
else
rXM = 0.0;
endif
if ( nY < nYDim )
rYP = a( nXDim + 1 + nX, nYDim + 1 + nY + 1, nZDim + 1 + nZ );
else
rYP = 0.0;
endif
if ( nY > -nYDim )
rYM = a( nXDim + 1 + nX, nYDim + 1 + nY - 1, nZDim + 1 + nZ );
else
rYM = 0.0;
endif
if ( nZ < nZDim )
rZP = a( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ + 1 );
else
rZP = 0.0;
endif
if ( nZ > -nZDim )
rZM = a( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ - 1 );
else
rZM = 0.0;
endif
n0 = a( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ );
aLap( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ ) = ...
( rXP + rXM + rYP + rYM + rZP + rZM - 6 * n0 ) / rH2;
endfor
endfor
endfor
endfunction


Obviously the solution of the differential equation should be the same as the numerical solution, so there is something wrong here. Based on my previous computations I suspect that the error is in the differential equation. Can somebody tell what is wrong?

• May I ask which tool are you using for the plots? – Algo Aug 13 at 14:52
• I use Octave. IMHO it is very good software. It contains a symbolic math package, too. – Tommi Höynälänmaa Aug 13 at 15:21
• The sign of the numerical and exact solutions are opposite. I would suspect a bug in your code. – cfdlab Aug 14 at 5:53
• You are solving the Poisson eqn in $R^3$ and the exact solution is in $L^2(R^3)$ but not in $H^1(R^3)$. – cfdlab Aug 14 at 6:19
• I have checked the numerical result by substituting it in the equation $\nabla^2 V_H = -4 \exp(-2r)$. What is the significance of the exact solution not being in $H^1(R^3)$? – Tommi Höynälänmaa Aug 14 at 8:56

The solution to the differential equation is $$V_H(r) = -\frac{K}{r} - \frac{r+1}{r} \exp(-2r).$$ By requiring that $$\lim_{r \to 0} V_H(r)$$ is real we get $$K=-1$$ and $$V_H(r) = \frac{1}{r} - \frac{r+1}{r} \exp(-2r).$$