Description
1 Introduction
We wish to solve for the currents in a resistor. The currents depend on the shape of the resistor and we also
want to know which part of the resistor is likely to get hottest.
2 The Resistor Problem
A wire is soldered to the middle of a copper plate and its voltage is held at 1 Volt. One side of the plate is
grounded, while the remaining are floating. The plate is 1 cm by 1 cm in size.
∂ φ
∂n = 0
∂ φ
∂n = 0
1 Volt
∂ φ
∂n = 0
As a result, current flows. The current at each
point can be described by a “current density” ~j. This
current density is related to the local Electric Field
by the conductivity:
~j = σ~E
Now the Electric field is the gradient of the potential,
~E = −∇φ
and continuity of charge yields
∇·~j = −
∂ ρ
∂t
Combining these equations we obtain
∇·(−σ∇φ) = −
∂ ρ
∂t
Assuming that our resistor contains a material of constant conductivity, the equation becomes
∇
2
φ =
1
σ
∂ ρ
∂t
For DC currents, the right side is zero, and we obtain
∇
2
φ = 0
3 Numerical Solution in 2-Dimensions
Laplace’s equation is easily transformed into a difference equation. The equation can be written out in
Cartesian coordinates
∂
2
Assuming φ is available at points
xi
, y j
we can write
∂ φ
∂ x
(xi
,y j)
=
φ(xi+1/2
, y j)−φ(xi−1/2
, y j)
∆x
and
∂
2φ
∂ x
2
(xi
,y j)
=
φ(xi+1, y j)−2φ(xi
, y j) +φ(xi−1, y j)
(∆x)
2
Combining this with the corresponding equation for the y derivitives, we obtain
φi, j =
φi+1, j +φi−1, j +φi, j+1 +φi, j−1
4
(1)
Thus, if the solution holds, the potential at any point should be the average of its neighbours. This is a very
general result and the above calculation is just a special case of it.
So the solution process is obvious. Guess anything you like for the solution. At each point, replace
the potential by the average of its neighbours. Keep iterating till the solution converges (i.e., the maximum
change in elements of φ is less than some tolerance).
But what do we do at the boundaries? At boundaries where the electrode is present, just put the value of
potential itself. At boundaries where there is no electrode, the current should be tangential because charge
can’t leap out of the material into air. Since current is proportional to the Electric Field, what this means
is the gradient of φ should be tangential. This is implemented by requiring that φ should not vary in the
normal direction.
4 Code
Import the libraries
First import the libraries needed by python.
2 hsetup 2i≡
from pylab import *
import mpl_toolkits.mplot3d.axes3d as p3
Define the Parameters
These are the parameters that control the run. They have the following defaults, which should be corrected
by the user via commandline arguments (see sys.argv)
Nx=25; // size along x
Ny=25; // size along y
radius=8;// radius of central lead
Niter=1500; // number of iterations to perform
Allocate the potential array and initialize it
First initialize φ to zero. Note that the array should have Ny rows and Nx columns. (the x and y are flipped
to make the solution fit the diagram). When you are done, φ should look like so:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
Now use “where” to locate the points in the region with a radius of 8 of the centre. To do this, you need
two arrays, one with x-coordinates and the other with y-coordinates with the same shape as the array. This
is done with meshgrid:
Y,X=meshgrid(y,x)
where x and y are vectors containing the x and y positions. Note that you should set x = 0 and y = 0 in the
middle of the region. Your condition to find the region that is 1 volt is then
X ∗X +Y ∗Y ≤ 0.35 ∗ 0.35
Use where and capture the elements that satisfy this condition in an index array “ii”. Note that ii will have
two columns, one giving the x-coordinate and the second the y-coordinate. Set the potential at these points
to 1 by
phi[ii]=1.0
Plot a contour plot of the potential in Figure 1. Mark the V = 1 region by marking those nodes red.
Perform the iteration
During the iteration, we will use a second potential array, called oldphi to hold the values of the previous
iteration. This is to keep track of how much the array changed during the current iteration.
You cannot create it by assignment:
oldphi=phi
The problem is that a python array is a pointer. So the assignment above merely creates a new pointer to
the existing data. Changing oldphi will change phi as well! That is not what we want. We want a new
memory area. The way to do this with an assignment operation is to use the copy member function:
oldphi=phi.copy()
3
To keep track of errors, allocate an error vector, errors with Niter elements.
Iterate Niter times and record the errors. The structure of your for loop should be as follows:
for k in range(Niter)
save copy of phi
update phi array
assert boundaries
errors[k]=(abs(phi-oldphi)).max();
#end
Updating the Potential
We first have to update the potential using Eq. 1. This is the central part of the algorithm. How can we do
it? In C, it is just a nested pair of for loops:
for(i=1 ; i
