Let's use SymPy to derive the relation between potential V and charge density R


```
%pylab inline
from sympy.interactive import init_printing
init_printing()
from sympy import pi, var, S, Piecewise, piecewise_fold
var("r R")
Vh = Piecewise((-S(2)/3 * pi * (3*R**2 - r**2), r <= R), (-S(4)/3 * pi * R**3 / r, True))
def laplace(f):
    return (r*f).diff(r, 2)/r
print "Vh ="
Vh
```

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





$$\begin{cases} - \frac{2}{3} \pi \left(3 R^{2} - r^{2}\right) & \text{for}\: r \leq R \\- \frac{4}{3} \frac{\pi R^{3}}{r} & \text{otherwise} \end{cases}$$



Charge density is then:


```
piecewise_fold(-laplace(Vh)/(4*pi))
```




$$\begin{cases} -1 & \text{for}\: r \leq R \\0 & \text{otherwise} \end{cases}$$




```
from numpy import (empty, pi, meshgrid, linspace, sum, sin, exp, shape, sqrt,
        conjugate)
from numpy.fft import fftn, fftfreq, ifftn

N = 100
print "N =", N
L = 2.4*4
R = 1.
x1d = linspace(-L/2, L/2, N+1)[:-1]
x, y, z = meshgrid(x1d, x1d, x1d, indexing="ij")

r = sqrt(x**2+y**2+z**2)
nr = empty(shape(x), dtype="double")
nr[:] = 0
nr[r <= R] = -1

Vanalytic = empty(shape(x), dtype="double")
Vanalytic[r <= R] = -2./3 * pi * (3*R**2 - r[r <= R]**2)
Vanalytic[r > R] = -4./3 * pi * R**3 / r[r > R]

ng = fftn(nr) / N**3

G1d = N * fftfreq(N) * 2*pi/L
kx, ky, kz = meshgrid(G1d, G1d, G1d)
G2 = kx**2+ky**2+kz**2
G2[0, 0, 0] = 1  # omit the G=0 term

tmp = 2*pi*abs(ng)**2 / G2
tmp[0, 0, 0] = 0  # omit the G=0 term
E = sum(tmp) * L**3
print "Hartree Energy (calculated): %.15f" % E


Vg = 4*pi*ng / G2
Vg[0, 0, 0] = 0  # omit the G=0 term
V = ifftn(Vg).real * N**3
V += Vanalytic[N/2, N/2, N/2] - V[N/2, N/2, N/2]
l2_norm = sum((Vanalytic - V)**2) 
print "l2_norm = ", l2_norm
plot(x[:, N/2, N/2], Vanalytic[:, N/2, N/2], label="analytic")
plot(x[:, N/2, N/2], V[:, N/2, N/2], label="FFT")
legend(loc="best");
```


```

```

950795722b6067af6c17737b0fa712a6b64ebbe8
ipynb
Jupyter Notebook
tutorial_exercises/FFT charged sphere.ipynb
certik/scipy-2013-tutorial
26a1cab3a16402afdc20088cedf47acd9bc58483
BSD-3-Clause
2015-02-28T08:53:05.000Z
2021-12-05T05:37:59.000Z
sympy/FFT charged sphere.ipynb
certik/scipy-in-13
418c139ab6e1b0c9acd53e7e1a02b8b930005096
BSD-3-Clause
2021-04-17T15:05:46.000Z
2021-04-17T15:05:46.000Z
sympy/FFT charged sphere.ipynb
certik/scipy-in-13
418c139ab6e1b0c9acd53e7e1a02b8b930005096
BSD-3-Clause
2015-03-11T00:25:21.000Z
2021-08-25T14:52:40.000Z
Qwen/Qwen-72B
1. YES
2. YES