coulomb_2D_A with a different parametrization of the surfaces.¶where ${\bf s}$ is Griffiths' script-r vector from source point $dq=\sigma\, da^\prime$ at position ${\bf r}^\prime$ to the observation point at ${\bf r}$, $\sigma$ is the surface charge density of the source object, and $da^\prime$ is the infinitesimal area on the surface.
Two dimensional surfaces can be parametrized in terms of two parameters --
lets call them $p_1$ and $p_2$. These parameters might by $x$ and $y$ for planar area at constant $z$, or
$\theta$ and $\phi$ for a spherical surface. (This is just a generalization of the parametrization
of a one-dimensional source in terms of the parameter $t$ in the notebook coulomb2.) The eventual
integrations will be double integrations over the parameters $p_1$, $p_2$ that are chosen for the
surface in question. We will need several functions that will depend on the surface charge:
$da$ will be of the form $da = \mbox{factor}\times dp_1 dp_2$. We need a function that gives this factor (that will depend on $p_1$ and $p_2$)
we need a function that will convert the parameters $p_1$ and $p_2$ into cartesian coord-nates
It is convenient to work with dimensionless parameters defined in terms of some characteristic length $R$ and charge $Q$ characterizing the source:
$$ {\bf r}^\ast \equiv \frac{\bf r}{R} \quad\quad {\bf s}^\ast \equiv \frac{\bf s}{R}\quad\quad \sigma^\ast = \frac{\sigma}{Q/R^2}\quad\quad \mbox{and} \quad\quad{\bf E}^\ast \equiv \frac{{\bf E}}{Q/(4\pi \epsilon_0 R^2)}. $$Coulomb's law for 2-d sources written in terms of these parameters is
$$ {\bf E}^\ast = \iint \frac{{\bf s}^\ast}{{s^\ast}^3}\sigma^\ast\, d{a^\prime}^\ast. $$This is really three separate double integrals, one for each component. Limits of integration must be determined from analysis of the source.
All calculations in what follow will be in terms of dimensionless variables; asterisks will
be implicit. The dimensional electric
field can be recovered by multiplying the numerical results below by the factor $Q/(4\pi\epsilon_0 R^2)$.
This notebook uses the quad method from the integrate submodule of scipy
for numerical integration. There are other choices, but this works well (and I like to use scipy
tools when they're available).
Marty Ligare, November 2022
import numpy as np
from scipy import integrate
import numdifftools as nd
import matplotlib as mpl
import matplotlib.pyplot as plt
# Following is an Ipython magic command that puts figures in the notebook.
# For figures in separate windows, comment out following line and uncomment
# the next line
# Must come before defaults are changed.
%matplotlib notebook
#%matplotlib
# As of Aug. 2017 reverting to 1.x defaults.
# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended,
# which don't seem to be universal
# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?
mpl.style.use('classic')
# M.L. modifications of matplotlib defaults using syntax of v.2.0
# More info at http://matplotlib.org/2.0.0/users/deflt_style_changes.html
# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]
plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size
Disk with radius $R$ and total charge $Q$ in $x$-$y$ plane centered on the origin. The obvious choice for length and charge scales are $R$ and $Q$. The dimensionless surface charge density is then
\begin{eqnarray*} \sigma^\ast &\equiv& \frac{\sigma}{Q/R^2} \\ &=& \frac{Q}{\pi R^2} \frac{1}{Q/R^2} \\ &=& \frac{1}{\pi} \end{eqnarray*}In this example, the parameters $p_1$ and $p_2$ could be chosen to be $x$ and $y$, or else $r$ and $\phi$, where $r$ is the radial position in the plane of the disk. In this notebook I'll use cartesian coordinates:
$$ p_1 \longrightarrow x\quad\quad\mbox{and}\quad\quad p_2\longrightarrow y. $$Note that in the cell Definition of functions describing source the arguments are written as cartesian variables to make interpretation convenient. The definitions in this cell will depend on the source and its pararmetrization under consideration. The cell below that, Source-independent functions, should not have to be modified from problem to problem, so it is written in terms the generic parameters $p_1$ and $p_2$.
def circle(x):
'''Function that provides upper limit for integration over y'''
return np.sqrt(1-x**2)
def ncircle(x):
'''Function that provides lower limit for integration over y'''
return -np.sqrt(1-x**2)
p1_1 = -1. # Lower limit for integration over x
p1_2 = 1. # Upper limit for integration over x
p2_1 = ncircle
p2_2 = circle
def sigma(x, y):
'''Charge density for disk'''
return 1/np.pi
def da(x, y):
'''The factor that multiplies dx*dy in da = dx*dy'''
return 1
def rp_x(x, y):
'''Conversion of (r,phi) to cartesian x-component of r-prime vector'''
return x
def rp_y(x, y):
'''Conversion of (r,phi) to cartesian y-component of r-prime vector'''
return y
def rp_z(x, y):
'''Conversion of (r,phi) to cartesian z-component of r-prime vector'''
return 0
def integrand_x(p1,p2,x,y,z):
'''x-component of field at (x,y,z) due to infinitesimal source dq at p1, p2
sx, sy, sz are components of vector separation between observation point and
source point'''
sx = x - rp_x(p1, p2)
sy = y - rp_y(p1, p2)
sz = z - rp_z(p1, p2)
return sx/(sx**2 + sy**2 + sz**2)**(3/2)*sigma(p1,p2)*da(p1,p2)
def integrand_y(p1,p2,x,y,z):
'''y-component of field at (x,y,z) due to infinitesimal source dq at p1, p2.
sx, sy, sz are components of vector separation between observation point and
source point'''
sx = x - rp_x(p1, p2)
sy = y - rp_y(p1, p2)
sz = z - rp_z(p1, p2)
return sy/(sx**2 + sy**2 + sz**2)**(3/2)*sigma(p1,p2)*da(p1,p2)
def integrand_z(p1,p2,x,y,z):
'''z-component of field at (x,y,z) due to infinitesimal source dq at p1, p2.
sx, sy, sz are components of vector separation between observation point and
source point'''
sx = x - rp_x(p1, p2)
sy = y - rp_y(p1, p2)
sz = z - rp_z(p1, p2)
return sz/(sx**2 + sy**2 + sz**2)**(3/2)*sigma(p1,p2)*da(p1,p2)
def ex(x,y,z):
'''Perform double integration for x-component'''
a = integrate.dblquad(integrand_x, p1_1, p1_2, p2_1, p2_2, args=(x,y,z))
return a[0]
def ey(x,y,z):
'''Perform double integration for y-component'''
a = integrate.dblquad(integrand_y, p1_1, p1_2, p2_1, p2_2, args=(x,y,z))
return a[0]
def ez(x,y,z):
'''Perform double integration for z-component'''
a = integrate.dblquad(integrand_z, p1_1, p1_2, p2_1, p2_2, args=(x,y,z))
return a[0]
ex(0,0,1), ey(0,0,1), ez(0,0,1)
(0.0, 0.0, 0.5857864376269019)
Analytical result for points on the $z$-axis:
$$ {\bf E}(0,0,z) =\frac{Q}{2\pi\epsilon_0 R^2} \left(1 - \frac{1}{\sqrt{1+\frac{R^2}{z^2}}}\right)\, \hat{\bf z} $$giving
\begin{eqnarray*} E_z(0,0,R) &=& 2\left(1 - \frac{1}{\sqrt{2}}\right)\frac{Q}{4\pi\epsilon_0 R^2}\\ &=& 0.5858 \times \frac{𝑄}{4\pi\epsilon_0 R^2} \end{eqnarray*}We get the same on-axis results as the analytical technique, but it's now trivial to get field at other points.
This surface can be parametrized in terms of the coordinates $x$ and $y$ along with the relationship
$$ \phi(x, y, z) = x^2 + y^2 + z^2 = 1. $$When using $x$ and $y$ we must be careful. In determing $dq = \sigma\, dA$ we must be using the infintesimal area on the surface, and not the area $dx\, dy$ which is an infintesimal area in the $x$-$y$ plane beneath the surface. As is demonstrated in Chapter 5 of Mathematical Methods in the Physical Sciences by Mary Boas, the correct area is
$$ dA = \frac{\left|\vec{\nabla}\phi\right|}{\frac{d\phi}{dz}}\, dx\, dy. $$It's straightforward to show that for the hemispherical surface
$$ dA = \frac{1}{\sqrt{1 - x^2 - y^2}}\, dx\, dy $$It would also be possible to use tools from the numdifftools package to evaluate numerically
the factor in the expression for $dA$.
def circle(x):
'''Function that provides upper limit for integration over y'''
return np.sqrt(1-x**2)
def ncircle(x):
'''Function that provides lower limit for integration over y'''
return -np.sqrt(1-x**2)
p1_1 = -1. # Lower limit for integration over x
p1_2 = 1. # Upper limit for integration over x
p2_1 = ncircle
p2_2 = circle
def sigma(x, y):
'''Charge density on hemisphere'''
return 1/np.pi/2
def da(x, y):
'''
The factor that multiplies dx*dy in the expression for da
'''
return 1/np.sqrt(1 - x**2 - y **2)
def rp_x(x, y):
'''Conversion of (theta,phi) to cartesian x-component of r-prime vector'''
return x
def rp_y(x, y):
'''Conversion of (theta,phi) to cartesian y-component of r-prime vector'''
return y
def rp_z(x, y):
'''Conversion of (theta,phi) to cartesian z-component of r-prime vector'''
return np.sqrt(1 - x**2 - y**2)
ex(0,0,2),ey(0,0,2),ez(0,0,2)
(0.0, 0.0, 0.36180339887479274)
I don't have a formula at hand, but I can think of one consistency check. Considering adding a second inverted hemisphere below the $x$-$y$ plane, making a complete hemispherical shell with total charge $2Q$. And note that by symmetry the fields due the bottom and top hemispheres are related:
$$ {\bf E}_{\rm bottom}(0,0,z) = -{\bf E}_{\rm top}(0,0,-z) $$The total field at $(0,0,2)$ should thus be
$$ {\bf E}_{\rm sphere}(0,0,2) = {\bf E}_{\rm top}(0,0,2) - {\bf E}_{\rm top}(0,0,-2) $$ez(0,0,2) - ez(0,0,-2)
0.4999999999995622
The field of the spherical shell should be that of a point charge $2Q$ located at the origin, which this is.
ez(0,0,0.9) - ez(0,0,-0.9)
9.827139102469573e-13
version_information is from J.R. Johansson (jrjohansson at gmail.com)
See Introduction to scientific computing with Python:
http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-0-Scientific-Computing-with-Python.ipynb
for more information and instructions for package installation.
If version_information has been installed system wide (as it has been on Bucknell linux computers with shared file systems), continue with next cell as written. If not, comment out top line in next cell and uncomment the second line.
%load_ext version_information
#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
version_information numpy, scipy, matplotlib
| Software | Version |
|---|---|
| Python | 3.7.13 64bit [GCC 7.5.0] |
| IPython | 7.31.1 |
| OS | Linux 4.9.0 9 amd64 x86_64 with debian 9.13 |
| numpy | 1.21.5 |
| scipy | 1.7.3 |
| matplotlib | 3.5.2 |
| Mon Nov 14 14:01:04 2022 EST | |