For a region of space without charges ($\rho = 0$) this reduces to Laplace's equation
$$ \nabla^2 \Phi(x, y, z) = 0 $$General solution by iteration: $$ \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 $$
Do not change $\Phi_{i,j}$ until a complete sweep has been completed.
def Laplace_Jacobi_slow(Phi):
# Don't use, very slow AND inefficient
Phi_new = Phi.copy()
Nx, Ny = Phi.shape
for xi in range(1, Nx-1):
for yj in range(1, Ny-1):
Phi_new[xi, yj] = 0.25*(Phi[xi+1, yj] + Phi[xi-1, yj]
+ Phi[xi, yj+1] + Phi[xi, yj-1])
Phi[:, :] = Phi_new
return Phi
Fast implementation using numpy array operations (vectorized, run at C speed, not Python speed) (see the board notes (PDF) or local PDF for an explanation for how to set up the numpy array operations that give you a 100-fold speed up over Laplace_Jacobi_slow()
):
def Laplace_Jacobi(Phi):
"""One update in the Jacobi algorithm"""
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])
return Phi
Immediately use updated new values for $\Phi_{i-1, j}$ and $\Phi_{i, j-1}$ (if starting from $\Phi_{1, 1}$).
Leads to accelerated convergence and therefore less round-off error (but distorts symmetry of boundary conditions... hopefully irrelevant when converged but check!)
def Laplace_Gauss_Seidel(Phi):
"""One update in the Gauss-Seidel algorithm"""
Nx, Ny = Phi.shape
for xi in range(1, Nx-1):
for yj in range(1, Ny-1):
Phi[xi, yj] = 0.25*(Phi[xi+1, yj] + Phi[xi-1, yj]
+ Phi[xi, yj+1] + Phi[xi, yj-1])
return Phi
Divide the lattice into a checkerboard of black and white cells. Update the odd cells first (like Jacobi) but then use the odd cells to update the even ones (Gauss-Seidel-like; see the board notes (PDF) or local PDF ). This leads to faster convergence and can be done with numpy array operations.
def Laplace_Gauss_Seidel_odd_even(Phi):
"""One update in the Gauss-Seidel algorithm on odd or even fields"""
# odd update (uses old even)
Phi[1:-2:2, 1:-2:2] = 0.25*(Phi[2::2, 1:-2:2] + Phi[0:-2:2, 1:-2:2] + Phi[1:-2:2, 2::2] + Phi[1:-2:2, 0:-2:2])
Phi[2:-1:2, 2:-1:2] = 0.25*(Phi[3::2, 2:-1:2] + Phi[1:-2:2, 2:-1:2] + Phi[2:-1:2, 3::2] + Phi[2:-1:2, 1:-2:2])
# even update (uses new odd)
Phi[1:-2:2, 2:-1:2] = 0.25*(Phi[2::2, 2:-1:2] + Phi[0:-2:2, 2:-1:2] + Phi[1:-2:2, 3::2] + Phi[1:-2:2, 1:-1:2])
Phi[2:-1:2, 1:-2:2] = 0.25*(Phi[3::2, 1:-2:2] + Phi[1:-2:2, 1:-2:2] + Phi[2:-1:2, 2::2] + Phi[2:-1:2, 0:-2:2])
return Phi
Solve the box-wire problem and make sure that the solution is converged to tol = 1e-3
.
Check convergence with the Frobenius Norm:
$$ ||\mathsf{\Phi}|| := \sqrt{\sum_{i,j}|\Phi_{ij}|^2} $$which is implemented as numpy.linalg.norm() when the argument is a matrix.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
Interactive plotting (with ipympl):
# for interactive work
%matplotlib widget
Only execute the next line if you don't want interactive plotting (e.g., when exporting to LaTeX/PDF or html):
%matplotlib inline
If you provide a filename then output is only written to a file and figures are close to conserve memory. This allows you to plot files in loops and later assemble them into movies using other programs such as ffmpeg, ImageMagick, mencoder, QuickTime 7, ...
def plot_contour(Phi, filename=None, zlabel=r"potential $\Phi$ (V)",
cmap=plt.cm.coolwarm):
"""Plot Phi as a contour plot.
Arguments
---------
Phi : 2D array
potential on lattice
filename : string or None, optional (default: None)
If `None` then show the figure and return the axes object.
If a string is given (like "contour.png") it will only plot
to the filename and close the figure but return the filename.
cmap : colormap
pick one from matplotlib.cm
"""
fig = plt.figure(figsize=(5,4))
ax = fig.add_subplot(111)
x = np.arange(Phi.shape[0])
y = np.arange(Phi.shape[1])
X, Y = np.meshgrid(x, y)
Z = Phi[X, Y]
cset = ax.contourf(X, Y, Z, 20, cmap=cmap)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_aspect(1)
cb = fig.colorbar(cset, shrink=0.5, aspect=5)
cb.set_label(zlabel)
if filename:
fig.savefig(filename)
plt.close(fig)
return filename
else:
return ax
def plot_surf(Phi, filename=None, offset=-20, zlabel=r'potential $\Phi$ (V)',
elevation=40, azimuth=-65, cmap=plt.cm.coolwarm):
"""Plot Phi as a 3D plot with contour plot underneath.
Arguments
---------
Phi : 2D array
potential on lattice
filename : string or None, optional (default: None)
If `None` then show the figure and return the axes object.
If a string is given (like "contour.png") it will only plot
to the filename and close the figure but return the filename.
offset : float, optional (default: 20)
position the 2D contour plot by offset along the Z direction
under the minimum Z value
zlabel : string, optional
label for the Z axis and color scale bar
elevation : float, optional
choose elevation for initial viewpoint
azimuth : float, optional
chooze azimuth angle for initial viewpoint
cmap : colormap
pick one from matplotlib.cm
"""
x = np.arange(Phi.shape[0])
y = np.arange(Phi.shape[1])
X, Y = np.meshgrid(x, y)
Z = Phi[X, Y]
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(X, Y, Z, rstride=2, cstride=2, linewidth=0.5, color="gray")
surf = ax.plot_surface(X, Y, Z, cmap=cmap, alpha=0.6)
cset = ax.contourf(X, Y, Z, 20, zdir='z', offset=offset+Z.min(), cmap=cmap)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel(zlabel)
ax.set_zlim(offset + Z.min(), Z.max())
ax.view_init(elev=elevation, azim=azimuth)
cb = fig.colorbar(surf, shrink=0.5, aspect=5)
cb.set_label(zlabel)
if filename:
fig.savefig(filename)
plt.close(fig)
return filename
else:
return ax
Use the fastest Poisson solver that we have at the moment: Laplace_Gauss_Seidel_odd_even()
Max_iter=30000
tol = 1e-3
Nmax = 100
Phi = np.zeros((Nmax, Nmax), dtype=np.float64)
Phi_old = np.zeros_like(Phi)
# initialize boundaries
# everything starts out zero so nothing special for the grounded wires
Phi[0, :] = 100 # wire at x=0 at 100 V
for n_iter in range(Max_iter):
Phi_old[:, :] = Phi
Phi = Laplace_Gauss_Seidel_odd_even(Phi)
DeltaPhi = np.linalg.norm(Phi - Phi_old)
if DeltaPhi < tol:
print("Laplace_Gauss_Seidel_odd_even converged in {0} iterations to {1}".format(n_iter+1, DeltaPhi))
break
else:
print("Laplace_Gauss_Seidel_odd_even did NOT converge in {0} iterations, DeltaPhi={1}".format(n_iter+1, DeltaPhi))
Laplace_Gauss_Seidel_odd_even converged in 7558 iterations to 0.0009995267090286367
Plot the result and visualy compare to what we had before:
plot_contour(Phi)
<matplotlib.axes._subplots.AxesSubplot at 0x11ba257b8>
plot_surf(Phi)
<matplotlib.axes._subplots.Axes3DSubplot at 0x11db512e8>
Accelerate convergence with the scheme
\begin{align} r_{i, j} &= \Phi_{i,j}^\text{new} - \Phi_{i, j}^\text{old}\\ \Phi_{i,j}^\text{new} &= \Phi_{i,j}^\text{old} + \omega r_{i,j} \end{align}where the new solution is computed with the Gauss-Seidel scheme.
Values of $1 \leq \omega \leq 2$ may work well, $\omega > 2$ can lead to numerical instabilities. Experiment!
Run to convergence (tol = 1e-03). Record the number of steps required (and visually check the solution).
Max_iter=10000
tol = 1e-3
Nmax = 100
omega = 1.99
Phi = np.zeros((Nmax, Nmax), dtype=np.float64)
Phi_old = np.zeros_like(Phi)
residual = np.zeros_like(Phi)
# initialize boundaries
# everything starts out zero so nothing special for the grounded wires
Phi[0, :] = 100 # wire at x=0 at 100 V
for n_iter in range(Max_iter):
Phi_old[:, :] = Phi
Phi = Laplace_Gauss_Seidel_odd_even(Phi)
residual[:, :] = Phi - Phi_old
DeltaPhi = np.linalg.norm(residual)
if DeltaPhi < tol:
print("SOR converged in {0} iterations to {1}".format(n_iter+1, DeltaPhi))
break
# SOR
Phi[:, :] = Phi_old + omega*residual # = omega*Phi + (1-omega)*Phi_old
else:
print("SOR did NOT converge in {0} iterations, DeltaPhi={1}".format(n_iter+1, DeltaPhi))
SOR converged in 3797 iterations to 0.0009987117944452034
plot_surf(Phi)
<matplotlib.axes._subplots.Axes3DSubplot at 0x121ff7d30>
Add a positive charge and a negative charge in the box.
Now we need to solve Poisson's equation.
$$ \nabla^2 \Phi(x, y, z) = -4\pi\rho(x, y, z)\\ $$Finite difference solution: $$ \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 $$
Modify the Laplace solvers to now solve Poisson's equation:
import numpy as np
def Poisson_Jacobi(Phi, rho, Delta=1.):
"""One update in the Jacobi algorithm for Poisson's equation"""
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 * Delta**2 * rho[1:-1, 1:-1]
return Phi
def Poisson_Gauss_Seidel(Phi, rho, Delta=1.):
"""One update in the Gauss-Seidel algorithm for Poisson's equation"""
Nx, Ny = Phi.shape
for xi in range(1, Nx-1):
for yj in range(1, Ny-1):
Phi[xi, yj] = 0.25*(Phi[xi+1, yj] + Phi[xi-1, yj]
+ Phi[xi, yj+1] + Phi[xi, yj-1]) \
+ np.pi * Delta**2 * rho[xi, yj]
return Phi
def Poisson_Gauss_Seidel_odd_even(Phi, rho, Delta=1.):
"""One update in the Gauss-Seidel algorithm on odd or even fields"""
a = np.pi * Delta**2
# odd update (uses old even)
Phi[1:-2:2, 1:-2:2] = 0.25*(Phi[2::2, 1:-2:2] + Phi[0:-2:2, 1:-2:2]
+ Phi[1:-2:2, 2::2] + Phi[1:-2:2, 0:-2:2]) + a * rho[1:-2:2, 1:-2:2]
Phi[2:-1:2, 2:-1:2] = 0.25*(Phi[3::2, 2:-1:2] + Phi[1:-2:2, 2:-1:2]
+ Phi[2:-1:2, 3::2] + Phi[2:-1:2, 1:-2:2]) + a * rho[2:-1:2, 2:-1:2]
# even update (uses new odd)
Phi[1:-2:2, 2:-1:2] = 0.25*(Phi[2::2, 2:-1:2] + Phi[0:-2:2, 2:-1:2]
+ Phi[1:-2:2, 3::2] + Phi[1:-2:2, 1:-1:2]) + a * rho[1:-2:2, 2:-1:2]
Phi[2:-1:2, 1:-2:2] = 0.25*(Phi[3::2, 1:-2:2] + Phi[1:-2:2, 1:-2:2]
+ Phi[2:-1:2, 2::2] + Phi[2:-1:2, 0:-2:2]) + a * rho[2:-1:2, 1:-2:2]
return Phi
Solve the Poisson problem with SOR (as above), but now employing the Poisson solvers (we are using the fastest one that we have):
Nmax = 100
Max_iter = 10000
omega = 1.99
Phi = np.zeros((Nmax, Nmax), dtype=np.float64)
Phi_old = np.zeros_like(Phi)
rho = np.zeros_like(Phi)
# initialize boundaries
# everything starts out zero so nothing special for the grounded wires
Phi[:, 0] = 100 # wire at y=0 at 100 V
rho[25:27, 39:41] = 5.0
rho[75:77, 39:41] = -5.0
Delta = 1.0
for n_iter in range(Max_iter):
Phi_old[:, :] = Phi
Phi = Poisson_Gauss_Seidel_odd_even(Phi, rho, Delta=Delta)
residual[:, :] = Phi - Phi_old
DeltaPhi = np.linalg.norm(residual)
if DeltaPhi < tol:
print("SOR converged in {0} iterations to {1}".format(n_iter+1, DeltaPhi))
break
# SOR
Phi[:, :] = Phi_old + omega*residual # = omega*Phi + (1-omega)*Phi_old
else:
print("SOR did NOT converge in {0} iterations, DeltaPhi={1}".format(n_iter+1, DeltaPhi))
SOR converged in 3823 iterations to 0.0009991605868815044
plot_surf(Phi, elevation=40, azimuth=20)
<matplotlib.axes._subplots.Axes3DSubplot at 0x1222e8cc0>
Given a potential $\Phi$, Poisson's equation gives the charge density: $$ \rho(\mathbf{x}) = -\frac{1}{4\pi} \nabla^2\Phi(\mathbf{x}) $$
Discretized Laplacian $$ \nabla^2 \Phi \approx \frac{\Phi(x+\Delta x,y) + \Phi(x-\Delta x,y) +\Phi(x,y+\Delta y) +\, \Phi(x,y-\Delta y) - 4\Phi(x, y)}{\Delta^2} $$
def laplacian2d(f, Delta=1):
"""Finite difference approximation of Del^2 f.
Arguments
---------
f : M x N matrix
Delta : float
Returns
-------
M x N matrix, boundaries set to 0
"""
L = np.zeros_like(f, dtype=np.float64)
L[1:-1, 1:-1] = f[2:, 1:-1] + f[:-2, 1:-1] + f[1:-1, 2:] + f[1:-1, :-2] - 4*f[1:-1, 1:-1]
return L/Delta**2
def laplacian2dsimple(f, Delta=1):
L = np.zeros_like(f, dtype=np.float64)
for i in range(1, L.shape[0]-1):
for j in range(1, L.shape[1]-1):
L[i, j] = f[i+1, j] + f[i-1, j] + f[i, j+1] + f[i, j-1] - 4*f[i, j]
return L/Delta**2
Making sure that our fancy numpy-version of the Laplacian is producing the same results as the naive implementation with Python loops:
import numpy as np
def test_laplacian2d():
ftest = np.random.random((200, 200))
assert np.allclose(laplacian2d(ftest), laplacian2dsimple(ftest))
test_laplacian2d()
Compute the charge density from the converged potential $$ \rho(x, y) = -\frac{1}{4\pi} \nabla^2 \Phi(x, y) $$
rhox = - laplacian2d(Phi)/(4*np.pi)
Does rhox
show the charges of $+5$ and $-5$ that we introduced with the charge density $\rho$?
print(rhox.min())
print(rhox.max())
-5.0000041867249845 5.000000000000001
Indeed, the point charges are recovered from the potential.
Plot the charge density:
plot_contour(rhox, zlabel=r"charge density $\rho$")
<matplotlib.axes._subplots.AxesSubplot at 0x1228fb898>
The position of the point charges are found correctly but the charge on the wire appear to be absent. This is expected because the Laplacian is not defined on the boundary. If you wanted to see the charge on the wire you would need to include it inside the integration domain and also assign a thickness (e.g., one or two grid cells).
The 3D plot shows the "point" charges as spikes of the size of one grid cell:
plot_surf(rhox, zlabel=r"charge density $\rho$", elevation=20, azimuth=20, offset=-3)
<matplotlib.axes._subplots.Axes3DSubplot at 0x1226e87b8>
The position of the point charges are found correctly but the charge on the wire appear to be absent. The problem is that the Laplacian is not defined on the boundary.
Also, if the charges are only defined on a single grid cell then the 3D visualization is sometimes not able to show it correctly.