EE2703 - Week 5 assignment on Python -  Laplace Equation - Solved

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.

                                                                                            by the conductivity:

                                                                                 Now the Electric field is the gradient of the potential,

1   Volt          
and continuity of charge yields

∇  j = −

∂t

Combining these equations we obtain 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 ∇·(−σ∇φ) = − 

∂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

 

Assuming φ is available at points xi,  we can write

                                                                                                      φ(xi+1/2,yj)−φ(xi−1/2,yj)

∂x x

and

 ∂x 

Combining this with the corresponding equation for the y derivitives, we obtain

                                                                                       φi,j                                                   (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()

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.

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 ith 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!

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.0014exp(−0.00226k)

Summing up the terms, we have

∞

Error = ∑ errork

k=N+1

∞

< ∑ AeBk

k=N+1

Z ∞

                                                                                                           ≈               AeBkdk

N+0.5

A

                                                                                                             = −      exp(B(N +0.5))

B

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.

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,li

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

 

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:

 

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

 

where the right side represents heat generated from  (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.