## 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<Ny-1 ; i++ ){ // interior nodes

for(j=1 ; j<Nx-1 ; j++ ){ // interior nodes

phinew[i,j]=0.25*(phi[i,j-1]+phi[i,j+1]

+phi[i-1,j]+phi[i+1,j]);

}

}

This would result in NxNy passes through the code and would reduce Python to a crawl. We have to do this

same work in a single line of code.

Python subarrays

Note how python specifies sub arrays:

• To extract a sub array of φ say φ[2,3] to φ[5,8], we use

phi[2:4,5:9]

Note that the last element of our range is not included in the subarray.

• To extract the last ten rows and first eight columns of φ we use

phi[-10:,0:8]

Note that phi[-10:-1,0:8] would have got the last ten rows but would have left out the very last

row.

With this information, we can now convert the for loops in C to vector operations in Python

Vectorizing For Loops

for(i=1 ; i<Ny-1 ; i++ ){ // interior nodes

for(j=1 ; j<Nx-1 ; j++ ){ // interior nodes

phinew[i,j]=0.25*(phi[i,j-1]+phi[i,j+1]

+phi[i-1,j]+phi[i+1,j]);

}

}

The way to “parallelize” this block of C code is to look at the slice of the matrix we work with. We are

updating all the interior elements of phi that is phi[1:-1,1:-1].

Now consider something rather strange: The matrix of left neighbours of this slice, i.e., the matrix created by phi[i,j-1]. Clearly it can be described as the matrix phi[1:Ny-2,0:Nx-3], or equivalently,

phi[1:-1,1:-2]. Note that the last row of the phi is Ny-1 and not Ny.

4

Proof: The (i, j)

th element of phi[1:Ny-1,1:Nx-1] is phi[i+1,j+1]. The (i, j)

th element of

phi[1:Ny-1,0:Nx-2] is phi[i+1,j], which is the left neighbour of phi(i+1,j+1).

Similarly construct the matrix of right neighbours, top neighbours and bottom neighbours. Once this is

done, the C code above can be written as a single line:

phi[1:-1,1:-1]=0.25*(phi[1:-1,0:-2]

+ phi[1:-1,2:]

+ top neighbours

+ bottom neighbours);

This is a matrix equation, which means that the for loops are still there, but they are executed in C within

the Python built in functions. Code in this line in Python.

Boundary Conditions

At boundaries where the electrode is present, just put the value of electrode potential itself. At boundaries

where there is no electrode, the gradient of φ should be tangential. This is implemented by requiring that φ

should not vary in the normal direction.

The code to enforce the boundaries has to go in after the Poisson update. At the electrode, we do not

have to do anything, since the edge row is not changed and the electrode potential is static. But at the other

boundaries, we have to make the final row the same as the adjacent row.

For example, at the left boundary, we need to set the normal derivitive of potential to zero. This means

that ∂ φ/∂ x = 0. This is achieved in C by

for( i=1 ; i<Ny-1 ; i++ )

phi[i,0]=phi[i,1];

In Python, we do the same thing by the following vector command:

phi[1:-1,0]=phi[1:-1,1]

What does this line do? The i

th entry of this vector equation says

phi[i,0]=phi[i,1]

which is what we want. Similarly, for the right side we have

for( i=1 ; i<Ny-1 ; i++ )

phi[i,Nx-1]=phi[i,Nx-2];

Convert it to Python code. Do the same for the top boundary.

The central portion is more difficult. During an iteration, the boundary values of the electrode can

get over-written. So we have to recreate those values. But you have already captured the elements of the

electrode in index vector ii, when setting up the potential. Then all you you have to do is:

phi[ii]=1.0

to restore that boundary condition.

Use Vectorized code!

5

Graph the results

First show how the errors are evolving. They should reduce but very slowly. Use a semi-log plot. You will

find that a log-log plot gives a reasonably straight line upto about 500 iterations, but beyond that, you get

into the exponential regime. This is for a 30 by 30 grid.

To see individual data points, plot every 50th point.

Anytime we have something like an exponential, we should extract the exponent. However, note that

it is an exponential decay only for larger iteration numbers. Let us try to fit the entire vector and only the

portion beyond 500. Remember that if the fit is of the form

y = AeBx

the thing to do is to take the log. Then we have

logy = logA+Bx

That is why it looks like a straight line in a semi-log plot.

Extract the above fit for the entire vector of errors and for those error entries after the 500th iteration.

Plot the fit in both cases on the error plot itself. use legend to label the curves (errors, fit1, fit2).

Stopping Condition

So where should you stop? Or more precisely, what is the error? The maximum error scales as

error = AeBk

where k is the iteration number. for our fit, For one choice of values, I got

errork = 0.0014 exp(−0.00226k)

Summing up the terms, we have

Error =

∞

∑

k=N+1

errork

<

∞

∑

k=N+1

AeBk

≈

Z ∞

N+0.5

AeBkdk

= −

A

B

exp(B(N +0.5))

Note that the inequality comes because the errors at any iteration could be either sign. So some of these

errors could cancel. We look for a worst case result and sum up the absolute values. For our values, this

expression evaluates to

Error < 0.63×0.03 = 0.02

Note that the last per-iteration change was 0.0000474, which shows you that that is extermely misleading.

Take a look at the error plot in detail. The error went from 6 × 10−4

to 1.5 × 10−4

in 500 iterations. As

you know, in exponential decay, what matters is the “half life” or the time constant. For this error vector,

the time constant is 1/B = 442. That explains the difference, since 0.02/442 = 0.0000452. The profile

was changing very little every iteration, but it was continuously changing. So the cumulative error was still

large.

A general comment on the algorithm: This method of solving Laplace’s Equation is known to be one of

the worst available. This is because of the very slow coefficient with which the error reduces.

6

Surface Plot of Potential

Next do a 3-D plot of the potential. For the default values I got the following plot (the boundary conditions

were slightly different with V2 = 0):

To create this plot, you need the following lines in your code:

fig1=figure(4) # open a new figure

ax=p3.Axes3D(fig4) # Axes3D is the means to do a surface plot

title(’The 3-D surface plot of the potential’)

surf = ax.plot_surface(Y, X, phi.T, rstride=1, cstride=1, cmap=cm.jet,linewidth=0, antialiased=False)

plot_surface does the real job of plotting. You can play with the options and see what they do.

The current is primarily going from electrode to electrode, but it spreads in the middle.

Contour Plot of the Potential

Obtain a contour plot of the potential. For this look up the contour function. Again, mark the central

electrode by red dots.

Vector Plot of Currents

Now obtain the currents. This requires computing the gradient. The actual value of σ does not matter to

the shape of the current profile, so we set it to unity. Our equations are

jx = −∂ φ/∂ x

jy = −∂ φ/∂ y

This numerically translates to

Jx,i j =

1

2

φi, j−1 −φi, j+1

Jy,i j =

1

2

φi−1, j −φi+1, j

Create the arrays Jx, Jy. Then call the quiver command:

quiver(y,x,Jy[::-1,:],Jx[::-1,:])

What this does is to create the vector plot in the desired direction.

Plot the current density using quiver, and mark the electrode via red dots. Here is what I got:

0 5 10 15 20 25

0

5

10

15

20

25 The vector plot of the current flow

The current fills the entire cross-section and then flows in the direction of the grounded electrode. Hardly

any current flows out of the top part of the wire. Why not?

Most of the current is in the narrow region at the bottom. That is what will get strongly heated. If you

feel motivated, take the currents calculated above, and now solve

∇·(κ∇T) = q =

1

σ

|J|

2

where the right side represents heat generated from J~· ~E (ohmic loss). The boundary condition is T = 300

at the wire and the ground while ∂T/∂n = 0 at the other three edges. This will now show us where the

heating is taking place (though it is obvious from the above plot), and how hot the plate will become.

8