Description
HONORS PHYSICS II
Project: Calculating Potentials Due to Charged Conductors: Numerical Solution of the Laplace’s
Equation Using the Relaxation Method1
1 Introduction
The techniques we have learned in class for calculating electric potential and electric fields are
useful only if the charge distribution is known. In many practical situations, however, the charge
distribution is not known. Instead, the value of the potential is known on the boundaries of a region.
As we know, in an electrostatic situation, the surface of a conductor is always an equipotential
surface, but the distribution of charge on the surface is in general not uniform and is not readily
calculated by the techniques we have learned (except for some simple problems that can be solved
with the method of images).
Consider a region of space that does not contain any electric charge in its interior and is enclosed
by one or more conductors maintained at fixed potentials (for example, by a battery). The question
that we are going to explore in this project is: How can we determine the potential as a function
of the position in this region? As we know, this questions is answered by solving the Laplace’s
equation ∇2V = 0 with the given boundary conditions. However, solving the equation analytically
is usually a difficult task; therefore we will do it numerically, with the help of a computer.
2 Method
The key to solving this problem numerically is to use the following fact about the electric potential:
In a region where there is no charge, the value of the electric potential at a given point is equal to
the average value of the potential at surrounding points. We will justify this statement2 by using
Gauss’s law in conjuction with the equation
E = −∇V, (1)
which gives the electric field components in terms of partial derivatives of the potential.
We will confine our discussion to situations, in which the potential depends only on two coordinates, say x and y. An example is the potential due to a long uniformly charged cylinder with
the z axis being the axis of symmetry. Here the potential at a point of space depends only on the
1This project description closely follows section 24-8 of H.D. Young, R.A. Freedman, University Physics (vol. 2),
10th edition.
2This is a general property of a class of functions called harmonic functions, which are solutions to the Laplace’s
equation. See Griffiths (section 3.1) for details.
1
point’s coordinates in a plane perpendicular to the axis of the cylinder, not on the coordinate z
along the axis. For such a two-dimensional situation, consider a point P with coordinates (x, y, z),
and enclose it by a Gaussian surface in the shape of a cubical box of side 2∆l, centered on P (see
Figure 1).
Figure 1: In a region where there is no charge, the value of the potential at a point P equals the
average of the potential values at points surrounding P.
If there is no charge in the volume enclosed by the box, the total electric flux ΦE through the
box is equal to zero. From Eq. (1) the z-component of the electric field Ez = −∂V /∂z equals zero
because the potential V is not a function of z. Hence, there is no flux through the two faces of the
Gaussian surface that are parallel to the xy-plane. Since the box is small, to a good approximation
the flux through each of the other faces of the box is equal to the product of the normal component
of E at the center of each face and the area (2∆l)
2 of each face. The total flux (equal to zero) can
then be expressed as
ΦE = Ex(x + ∆l, y, z)(2∆l)
2 + (−Ex(x − ∆l, y, z)) (2∆l)
2 +
+Ey(x, y + ∆l, z)(2∆l)
2 + (−Ey(x, y − ∆l, z)) (2∆l)
2 = 0. (2)
Using Eq. (1), we can write the electric-field components to the same approximation as
Ex(x + ∆l, y, z) ≈ −
V (x + ∆l, y) − V (x, y)
∆l
,
Ex(x − ∆l, y, z) ≈ −
V (x, y) − V (x − ∆l, y)
∆l
,
Ey(x, y + ∆l, z) ≈ −
V (x, y + ∆l) − V (x, y)
∆l
,
Ey(x, y − ∆l, z) ≈ −
V (x, y) − V (x, y − ∆l)
∆l
. (3)
Substituting Eqns. (3) into Eq. (2) and dividing through by 4∆l, we obtain
−[V (x + ∆l, y) − V (x, y)] + [V (x, y) − V (x − ∆l, y)] +
−[V (x, y + ∆l) − V (x, y)] + [V (x, y) − V (x, y − ∆l)] = 0.
If we solve this for V (x, y), the potential at point P, we find
V (x, y) = 1
4
[V (x + ∆l, y) + V (x − ∆l, y) + V (x, y + ∆l) + V (x, y − ∆l)] . (4)
2
In words, the value of the potential at P is the average of the potential values at the points
surrounding P. This statement becomes exact in the limit that ∆l becomes infinitesimally small.
To see how to use Eq. (4) to calculate the potential due to a set of charged conductors, let us
consider a specific situation. Figure 2a shows a hollow conducting box with a square cross section
Figure 2: (a) A conducting box whose z-dimension is much longer than the dimension L along the
x- and y-axes. (b) A view down the z-axis of the two segments into which the box is divided and
the potential difference between the segments.
and with a long axis parallel to the z-axis. The length of the box is very much greater than the
dimension L. The top of the box, labelled a, is insulated from the other three sides, collectively
labelled b; this is done by having the top be a separate piece of metal with a small gap between it
and the vertical sides of the box (Figure 2b). A fixed potential difference V0 is maintained between
segments a and b of the box. We choose the potential of the lower segments to be Vb = 0, so
the potential of the upper segment is Va = V0. As a result of the potential difference, there is a
positive charge on a (the higher-potential conductor) and a negative charge on b (the lower potential
conductor).
Our goal is to find the potential V at all points of the interior volume of the box. Because the
box is long, the potential inside the box is, to a good approximation, a function of x and y only. We
imagine making a rectangular grid of points inside the box, separated by a distance ∆l (Figure 3).
The outermost points of the grid are on the conductor surfaces themselves. Eq. (4) then relates the
potential at different grid points to each other; since the potentials of the conductors are specified,
we can determine the potential at each grid in the empty interior of the box.
Figure 3: A view down the long axis of the conducting box showing a rectangular grid of points
separated by ∆l.
3
Complications arise because Eq. (4) relates the potential at four different grid points. The value
of V is unknown at each grid point in the interior of the box, and each such point is surrounded by
two or three other interior points at which V is also unknown. So Eq. (4) cannot be used to solve
for the values of the potential at these interior points in a single step. Instead, we need to use an
iterative method: we will make a series of successive approximations to find a set of values of V at
the interior points such that Eq. (4) is satisfied at every point. The procedure that we follow to do
this is called the relaxation method.
Here is the algorithm to carry out this calculation
1. Choose a positive3 value of the potential difference V0.
2. Choose the number m of grid points4 across or down the region shown in Figure 3. The total
number of grid points is then m2
, and the total number of grid point in the interior of the box
is (m − 2)2
.
3. Let (j, k), where 1 ≤ j ≤ m and 1 ≤ k ≤ m be a pair of integer indexes that identify a
particular grid point and its location in the grid (the jth column and the kth row of the grid).
4. Assign an initial value to the potential V (j, k) for each grid. For k = 1 (the top row, corresponding to the surface of the upper conductor a in Figures 1a and 2), set V (j, k) = V0. If
j = 1, or j = m, or k = m, set V (j, k) equal to zero; these correspond to the left, right, and
bottom surfaces, respectively, of the lower conductor b in Figures 1a and 2.
For all other values of (j, k), corresponding to grid points in the empty interior of the box,
assign an arbitrary value of V (j, k). The closer this arbitrary value is to the actual value,
the fewer iterations will be required to obtain a good solution. But any value greater than
zero and less than V0 will work. A good choice for the arbitrary value might be V0/2, for all
interior grid points. Do not use V = 0, as this choice will cause problems in Step 7.
5. Specify the desired accuracy of the results. The algorithm will iterate to find a solution for
V (j, k) at grid points inside the box, and the iteration will stop when the relative change in
the values from one iteration to the next ε is less than the desired accuracy.5 The smaller the
value chosen, the more iterations will be required.
6. Use Eq. (4) to compute new values of the potential
Vnew(j, k) = 1
4
[V (j + 1, k) + V (j − 1, k) + V (j, k + 1) + V (j, k − 1)]
for all grid points in the interior.
7. Repeat Step 6 until the value
Vnew(j, k) − V (j, k)
V (j, k)
is less that the required accuracy ε for all grid points in the interior.
It is easy to modify this algorithm to treat other types of conducting boundaries. The values
of the potential on the conducting surfaces can be chosen at will. The box can be changed from a
3The procedure runs into a trouble if the potential difference is negative.
4Values of m between 20 and 500 work well. Large values of m require lengthy calculation; small values of m give
low resolution.
5Reasonable values of ε are from 0.01 to 0.001.
4
square to a rectangle by having a different number of rows than columns. The box can also have
shapes other than rectangular.
Once the values of the potential have been calculated for all grid points, the electric field components Ex and Ey can be calculated as well, using Eq. (3) and the central difference formula for
the derivative
Ex(x, y) = −
∂V
∂x
(x,y)
≈ −
V (x + ∆l, y) − V (x − ∆l, y)
2∆l
, (5)
Ey(x, y) = −
∂V
∂y
(x,y)
≈ −
V (x, y + ∆l) − V (x, y − ∆l)
2∆l
. (6)
(7)
3 Project tasks
For each of the following systems of conductors and two different choices of the grid density (give
the number of rows and columns) complete the tasks listed below. Formulate observations and
conclusions.
1. Specify the initial values of the potential inside the region.
2. Find the electric potential inside the region with the assumed accuracy ε, starting with the
initial configuration from the previous step. Plot the final result by visualizing the electric
potential and equipotential lines in the form of a density plot and contour plot.6
3. Find the corresponding electric field and visualize it (both magnitude and direction).7
4. Choose the dense grid and repeat your calculations (only for the electric potential) for three
different values of ε (differing by orders of magnitude). Plot the results in the form of density
plots and comment on them.
5. Include a table with the following data: grid size, accuracy, and the number of iterations
needed to complete the calculation. Optionally (this will not affect your grade), you may
also give the time needed to complete it; in that case please include information about your
computer’s configuration.
6. Choose the dense grid and illustrate progress of the algorithm by plotting (one next to another)
density plots with equipotential lines for the electric potential at four different stages of the
algorithm (including the initial and the final stage).
7. ∗ Bonus (up to 8 bonus marks): As an extension of the previous task, animate progress of the
relaxation algorithm (you may us an animated gif image format to do this). If you decide to
do this task, please indicate it in your report and upload the animation to Canvass.
8. Discuss your results.
6See Wolfram’s Mathematica command DensityPlot or ContourPlot to see what a density/contour plot looks like.
7Your visualization should look like the output of the VectorPlot command in Mathematica.
5
3.1 System A
The top a of the box in Figure 2 is at potential Va = 1 V, and the other walls b at V = 0.
3.2 System B
The terminals of the battery or voltage source shown in Figure 2 are reconnected, so that the top
of the box is at Va = 0 and the other walls at Vb = 1 V.
3.3 System C
The box shown in Figure 2 is modified so that the upper surface (which is at Va = 1 V) has a
”finger” that projects down in the interior of the box as shown in Figure 4. The other walls are at
Vb = 0. Discuss your results and compare them to those obtained for System A.
Figure 4: System C: A box with a ”finger”.
4 Deliverables
After completing all tasks you need to write a self-contained report presenting your results. It
should be starting with a short introduction, similar to the opening sections of this document. The
report should be typed (LATEX is recommended, but not compulsory) and submitted in printed form,
including graphs, with all pages stapled together.
You will need to implement the numerical method studied in this project. There are no restrictions imposed on the programming environment (C++, Fortran, Matlab, Mathematica, Octave,
Maxima,… are all allowed). Please keep your code clear, simple, and concise. Your source code
files, the report (pdf format only), and (optional) animation files should be packed into a single ZIP
file and uploaded into the Canvas system by the due date.
In this project, you will be working in groups of four. You need to make sure that workload
is distributed uniformly among all group members. On the title page, please include the following
statement: ”We state that each of us has contributed equally to this project.” followed by signatures
of all group members. The code should be also labelled with the names of group members. Each
group should submit a single report and a single ZIP file.
6
5 Honor Code
The JI Honor Code applies to the project in the same way it does to homework. In particular,
according to the current version of the Honor Code8
Assignments involving collaboration within a group (e.g., lab reports, project reports,
collaborative course work) require that all members of the group whose name appears on
the assignment are jointly and fully responsible for the entirety of the submitted work.
If any section of the submission is found to violate the Honor Code, all group members
whose name appears on the submission are equally and jointly liable for the violation.
8
http://umji.sjtu.edu.cn/academics/academic-integrity/honor-code/, as of 21 October 2016
7
