Incomplete – complete notebook 15_PDEs-2.ipynb will be available after the class.
Still solving Poisson's equation for the electric potential $\Phi(\mathbf{r})$ and the charge density $\rho(\mathbf{r})$:
$$ \nabla^2 \Phi(x, y, z) = -4\pi\rho(x, y, z)\\ $$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] = NotImplemented
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):
%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
Complete the code (use one of the fast numpy Laplace solvers) and get converged results (to tol = 1e-03
or better). Try
Laplace_Jacobi()
Laplace_Gauss_Seidel_odd_even()
Which one converges faster?
Max_iter=10000
tol = 1e-3
Nmax = 100
nsave = 100 # print a progress message every nsave steps
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
# Laplace solver
Phi =
# check convergence
DeltaPhi =
if DeltaPhi < tol:
print("Laplace solver converged in {0} iterations to {1}".format(n_iter+1, DeltaPhi))
break
if n_iter % nsave == 0:
print("Iteration {0}".format(n_iter), end="\r")
else:
print("Laplace solver did NOT converge in {0} iterations, DeltaPhi={1}".format(n_iter+1, DeltaPhi))
Do a quick wireframe plot:
# plot Phi(x,y)
x = np.arange(Nmax)
y = np.arange(Nmax)
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)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel(r'potential $\Phi$ (V)')
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, 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(r"potential $\Phi$ (V)")
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()
Plot the result $\Phi(x, y)$ (i.e., Phi
) and visualy compare to what we had before:
plot_contour(Phi)
plot_surf(Phi)
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.
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)
# implement SOR and convergence check
else:
print("SOR did NOT converge in {0} iterations, DeltaPhi={1}".format(n_iter+1, DeltaPhi))
plot_surf(Phi)
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 (replace the NotImplemented
with the appropriate term):
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]) \
+ NotImplemented
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]) \
+ NotImplemented
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]) + NotImplemented
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]) + NotImplemented
# 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]) + NotImplemented
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]) + NotImplemented
return Phi
Solve the Poisson problem with SOR (as above):
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))
plot_surf(Phi, elevation=40, azimuth=20)
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} $$
Implement the Laplacian operator (replace NotImplemented
):
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):
NotImplemented
return L/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] = NotImplemented
return L/Delta**2
rhox = - laplacian2dsimple(Phi)/(4*np.pi)
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())
Plot the charge density:
plot_surf(rhox, zlabel=r"charge density $\rho$", elevation=20, azimuth=20, offset=-3)
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 visualization is sometimes not able to show it correctlty.)
import numpy as np
def test_laplacian2d():
ftest = np.random.random((200, 200))
assert np.allclose(laplacian2d(ftest), laplacian2dsimple(ftest))
test_laplacian2d()