Initialize printing:



In [1]:

    
%pylab inline
from sympy.interactive.printing import init_printing
init_printing()









    



Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

Poisson Equation

Radial Poisson equation is $$ \phi''(r)+{2\over r} \phi'(r) = -4\pi\rho(r). $$ Alternatively, this can also be written as: $$ {1\over r}(r \phi(r))'' = -4\pi\rho(r). $$

Example I

Positive Gaussian charge density:



In [2]:

    
from sympy import var, exp, pi, sqrt, integrate, refine, Q, oo, Symbol, DiracDelta
var("r alpha Z")
alpha = Symbol("alpha", positive=True)
rho = Z * (alpha / sqrt(pi))**3 * exp(-alpha**2 * r**2)
rho









    Out[2]:





$$\frac{Z \alpha^{3} e^{- \alpha^{2} r^{2}}}{\pi^{\frac{3}{2}}}$$

The total charge is $Z$:



In [3]:

    
integrate(4*pi*rho*r**2, (r, 0, oo))









    Out[3]:





$$Z$$

Solve for $\phi(r)$ from the Poisson equation:



In [4]:

    
phi = integrate(-4*pi*rho*r, r, r)/r
phi









    Out[4]:





$$\frac{Z \operatorname{erf}{\left (\alpha r \right )}}{r}$$

Tell SymPy that $\alpha$ is positive:



In [5]:

    
phi = refine(phi, Q.positive(alpha))
phi









    Out[5]:





$$\frac{Z \operatorname{erf}{\left (\alpha r \right )}}{r}$$

Plot the charge:



In [6]:

    
from sympy import plot



In [7]:

    
plot(rho.subs({Z: 1, alpha: 0.1}), (r, 0, 100));

Plot the potential



In [8]:

    
plot(phi.subs({Z: 1, alpha: 0.1}), 1/r, (r, 0, 100), ylim=[0, 0.12]);



In [8]: