Partial Differential Equations 1: Boundary Value Problems

In this lecture, we will explore numerically solving partial differential equations (PDEs) that have the form of boundary value problems using various techniques of relaxation.

Learning objectives: After completing this lecture, you should be able to:

  1. Implement the Jacobi and Gauss-Seidel methods of solving PDEs
  2. Speed up calculation by implementing successive over-relaxation (SOR)
  3. Solve for the electrostatic potential $\phi$ for 2-dimensional (2D) problems using SOR including fixed charges
  4. Calculate the surface charge density on the surfaces of a simulation of the 2D poission equation in physical units (C/m$^2$)
In [ ]:
# Notebook code
import numpy as np
import matplotlib.pyplot as plt
from time import time
import sys
sys.path.append('resource/asnlib/public/')
sys.path.append('../answer_checking/') 
from validate_answers import *
In [ ]:
%matplotlib notebook

Poisson's equation in 2D

We will look specifically in this course at the solutions to a specific, linear partial differential equation known as Poisson's equation:

$$ \nabla^2 \phi = \frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} + \frac{\partial^2\phi}{\partial z^2} = \rho(x,y,z) $$

This equation is hopefully well known to you from your course in electrostatics: is is the equation that determines the electrostatic potential $\phi(\vec{r})$ for a given distribution of charges $\rho(\vec{r})$.

In this lecture, we will discuss how to solve a two-dimensional version of the equation. By working in 2D, we will keep the calculation times relatively short and it will be easy to visualize our results. The technique, however, can easily be extended into 3D in an obvious way.

What does this mean physically? Of course, Maxwell's equations apply to problems in 3-dimensions: we can't just cut one of the dimensions off! However, if our system is translationally invariant in one dimension (say, the $z$ direction for example), then we know that $\phi$ cannot depend on $z$ since if we move in the $z$ direction, everything looks the same. In this case, the derivative in the $z$ direction will be zero:

$$ \frac{\partial^2\phi}{\partial z^2} = 0 $$

leading us to a two-dimensional version of the Poisson equation:

$$ \frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} = \rho(x,y) $$

What would something like this look like physically? A good example is the two wires that make up a simple transmission line:

If we approximate the wires as infinitely long, or are interested only in places near the wires and also far from their ends, the 2D Poisson equation above will give us a good approximation of the electrostatic potential in the plane perpendicular to the wire direction.

We will start with the simple case where $\rho(x,y) = 0$, leading to Laplace's equation:

$$ \frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} = 0 $$

We will use this to explore various techniques for solving it using finite differences and different relaxation techniques, then move on to the solving physical problems in electrostatics.

Finite Difference method

The first step in solving the problem numerically is to map the continuous $x,y$ plane into discrete points. Here, we will consider discrete points with a spacing $a$.

The next step is to approximate the partial derivatives with finite differences. Doing so will lead to:

$$ \frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} \approx \frac{ \phi(x+a,y) + \phi(x-a,y) + \phi(x,y+a) + \phi(x,y-a) - 4 \phi(x,y)} {a^2} $$

All we need to do now is to equate this to zero and find the solutions of the resulting coupled (linear) equations. How many equations will we have? If we choose a grid of 100x100 points in our $x,y$ plane, we will have 10,000 equations! That's a lot! Coding this into a matrix is quite daunting...

Jacobi method

Fortunately, there is a simpler method for solving these 10,000 equations, which is known as "relaxation".

The problem we are looking at is a boundary-value problem: to be able to solve for the potential on our grid, we need to know the value of the potential on the boundaries of our simulation. This is our "known" starting point, and we need to then find the solutions in the middle that solve the Laplace (or the Poisson) equation.

The way relaxation works is the following: we start with our fixed boundary condition on the outside of our simulation "box". We then start with an initial guess at the potential everywhere inside the box. What do we choose for this? It turns out that it doesn't matter what we choose: the techniques we will look at will always work no matter what you chose (which is handy!).

(Note that here, we will consider only "metallic" fixed-potential boundary conditions, also know as the Dirichlet condition. One alternative is also the Neumann boundary condition, specifying no electric field perpendicular to the boundary, or periodic boundary conditions. In addition, there are more exotic boundary conditions that allow, among other things, the simulation of an open boundary to infinity.)

Once we have our boundary set and our initial guess, we then iterate through all the points in our simulation and replace the value of the potential at iteration $N+1$ with the value that would solve the equation we are solving, given the value of it's neighbors at that step of the iteration $N$. With $\rho(x,y) = 0$, this becomes:

$$ \phi_{N+1}(x,y) = \frac{\phi_N(x+a,y) + \phi_N(x-a,y) + \phi_N(x,y+a) + \phi_N(x,y-a)}{4} $$

Note that for $\rho = 0$, the lattice spacing $a$ falls out: we don't even need to know the spacing to do the calculation, just the aspect ratio.

Miraculously (?), if we just keep doing this for a while, the potential will magically converge to the correct solution!

How do we know when to stop? One way is to keep track of the largest change in $\phi$ of any of the pixels in our simulation. If this change in potential goes below a target accuracy that we configure, then we can stop iterating our simulation.

Exercise 1(a) Consider a square box with the left, right, and bottom sides grounded, and the top side set to a voltage of 1V. In our finite element technique, we will divide the two-dimensional $x-y$ space into an $M \times M$ grid, illustrated here for $M = 11$:

The outer boundary points are fixed by the boundary condition, and are not updated when we are solving for $\phi$: instead they remain at a fixed potential. We need to keep track of this so that we do not accidentally update any points that should correspond to a fixed external potential.

Note that our update equation above does not depend on the spacing between our grid points: the Laplace equation has the same solution regardless if the size of the box in 1 meter of 1 nm! For the Laplace equation, we do not have to keep track what the spacing between our points is.

Your task is to fill in the code below to use Jacobi's method to calculate $\phi(x,y)$.

During your calculation, keep track of the absolute value of the largest change delta_max in $\phi$ during that iteration in a list so that we can plot later the convergence as a function of the number of iterations. Use a threshold of 0.1 mV for $\phi$.

To benchmark the computational speed of our simulation, we will keep track of the running time and the code speed.

In [ ]:
target_accuracy = 1e-4

# The size of the simulation grid
M = 101
phi = np.zeros([M,M])
phi[-1,1:-1] = 1

# A second matrix for storing the values calculated for the next
# iteration
phinew = phi.copy()

# A matrix to keep track of which points in our simulation 
# are fixed voltages. In this case, it is only the boundaries
# of the simulation
fixed = np.zeros([M,M])
fixed[0,:] = 1
fixed[-1,:] = 1
fixed[:,0] = 1
fixed[:,-1] = 1

# Keep track of the biggest update. 
# When this is smaller than our target accuracy, 
# we will stop iterating
delta_max = 1
delta_max_list = []

t1 = time()
while delta_max > target_accuracy:
    for i in range(M):
        for j in range(M):
            if fixed[i,j]:
                #phinew[i,j] = ...
                ###
### Your code here
###

            else:
                #phinew[i,j] = ....
                ###
### Your code here
###

    #delta_max =....
    ###
### Your code here
###

    delta_max_list.append(delta_max)
    print("N_iter %d delta_max %e\r" % (len(delta_max_list), delta_max), end='')
    # Now that we're done, phi becomes phinew. A sneaky trick: if we swap the two, 
    # we don't need to copy the whole array! 
    phi,phinew = phinew,phi
t2 = time()

print("\nTotal running time: %.2f min" % ((t2-t1)/60))
print("Code speed: %.1f iterations per second" %(len(delta_max_list)/(t2-t1)))

answer_13_1a_1 = phi.copy()
answer_13_1a_2 = delta_max_list.copy()
In [ ]:
question = "answer_13_1a"
num = 2

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

To explore your data, you can use this function to make a handy interactive "data explorer" that I wrote (see code in resource/asnlib/public/explore_image.py if you are interested) that can show you both color scale plots and also line cuts:

In [ ]:
from explore_image import *

Now let's use this function to explore the $\phi$ that you calculated:

For the interactive plots, I will have selected the matplotlib notebook driver because it is much faster. However, you must not forget to turn off the "power button" in the top right corner before moving on to the next code cells.

In [ ]:
explore_image(answer_13_1a_1)

Exercise 1(b): Make a plot of delta_max as a function of the iteration number. The y-axis should be a log scale. As always, your plot should have appropriate labels with units.

In [ ]:
###
### Your code here
###

Gauss-Seidel method

A speedup of the Jacobi method can be obtained by using the already updated calculated values of $\phi$ instead of waiting for the next round.

For doing this, we have to have a different way of keeping track of delta_max than we used above, and so the code will look a bit different. A handy (and computationally efficient) way of doing this in python is to keep track of all the delta in a matrix and then us np.max(np.abs(delta)) to find the largest one afterwards.

(You might think: I could add an if statement in my loop to keep track of delta_max: why do I need a full matrix? This is what I initially thought too but it turns out that this slows your code down a lot and you're much faster just keeping the full matrix and using np.max() at the end. Why? Python is an interpreted language, which means adding anything to a hand-coded loop is just a really, really bad idea if you don't have to. Because np.max() does the calculation with "vectorized": it runs with pre-compiled, and not interpreted, loops. We saw this in the first lecture with the 700-fold speed improvement of np.average compared to writing our own for loop to calculate the average of a vector.)

Exercise 2(a) Perform the calculation of exercise 1(a) with the Gauss-Seidel method.

In [ ]:
target_accuracy = 1e-4

M = 101
phi = np.zeros([M,M])
delta = phi.copy()
phi[-1,1:-1] = 1

# A matrix to keep track of which points in our simulation 
# are fixed voltages. In this case, it is only the boundaries
# of the simulation
fixed = np.zeros([M,M])
fixed[0,:] = 1
fixed[-1,:] = 1
fixed[:,0] = 1
fixed[:,-1] = 1

# Keep track of the biggest update. 
# When this is smaller than our target accuracy, 
# we will stop iterating
delta_max = 1
delta_max_list = []

t1 = time()
while delta_max > target_accuracy:
    for i in range(M):
        for j in range(M):
            if not fixed[i,j]:
                #delta[i,j] = ...
                #...
                ###
### Your code here
###

    delta_max = np.max(np.abs(delta))
    delta_max_list.append(delta_max)
    print("N_iter %d delta_max %e\r" % (len(delta_max_list), delta_max), end='')   
t2 = time()
print("\nTotal running time: %.2f min" % ((t2-t1)/60))
print("Code speed: %.1f iterations per second" %(len(delta_max_list)/(t2-t1)))

answer_13_2a_1 = phi.copy()
answer_13_2a_2 = delta_max_list.copy()
In [ ]:
question = "answer_13_2a"
num = 2

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

You can check that it looks the same:

In [ ]:
# Notebook code
explore_image(answer_13_2a_1)

Although the code is not that much faster in total time (because of the slightly more complicated stuff going on in the for loops), it converges in fewer iterations.

Exercise 2(b): Plot delta_max vs iteration number for both the Jacobi method and the Gauss-Seidel method. Use a logscale for the y-axis.

In [ ]:
###
### Your code here
###

Successive over-relaxation

Although we did not gain much speed with the Gauss-Seidel technique, there is a small trick we can play with it to make it MASSIVELY faster.

The trick is as follows: instead of correcting phi to the value that would solve the equation based on it's neighboring values, we go a bit further:

$$ \phi_{new}(x,y) = \phi_{old}(x,y) + \omega \delta $$

where $\delta$ is the correction that you would have calculated in the Gauss-Seidel method and $\omega$ is a number between 0 and 2. With $\omega > 1$, this technique is known as successive over-relaxation (SOR). It turns out, if you pick the right $\omega$, you simulation can run MUCH MUCH faster.

Picking the optimal value of $\omega$ is not easy, and depends in detail on the size and even configuration of your simulation. The Gauss-Seidel method corresponds to $\omega = 1$. For $\omega > 2$, the method does not converge (it is unstable). And, for nonlinear Poisson equations (in which $\rho$ depends on $\phi$ in a nonlinear way, which I studied in my PhD), under-relaxing with $\omega < 1$ can help unstable simulations converge. In general, for regular linear PDEs and large simulations, a value close to 2 is a good guess, and guessing a bit high is better than guessing to low. But even if you guess a bit wrong, it still can speed up your simulation very significantly.

(In particular, it can change the scaling with simulation size from $N^2$ to $N^{1.5}$, which is a huge deal! You can read chapter 5 of my PhD thesis if you would like to learn some more about the speed and applications of relaxation technique.)

As you will see, even for these small simulation sizes, it can significantly speed up the calculation.

Exercise 3(a): Modify your code from 2(a) to implement SOR. Use an SOR paramter of 1.95 (which I found worked well by trail-and-error...).

In [ ]:
SOR_parameter = 1.95
target_accuracy = 1e-4

###
### Your code here
###


print("\nTotal running time: %.2f min" % ((t2-t1)/60))
print("Code speed: %.1f iterations per second" %(len(delta_max_list)/(t2-t1)))

answer_13_3a_1 = phi.copy()
answer_13_3a_2 = delta_max_list.copy()
In [ ]:
question = "answer_13_3a"
num = 2

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

By now, I guess you trust that the technique works, but you can also take a look at $\phi$ just to check if you want:

In [ ]:
# Notebook code 
explore_image(answer_13_3a_1)

Exercise 3(b): Make a plot of delta_max vs iteration number for the three techniques.

In [ ]:
###
### Your code here
###

Note that SOR initially is slower than the other two techniques (and sometimes the convergence parameter can even INCREASE at the start with SOR), but it quickly overtakes th other two and converges very quickly to the correct answer.

Including fixed charges

In some cases, one might want to simulate the case where there are external fixed charges in our simulation.

To do this, we need to add an additional term to our update equations to account for the external charge. In the Jacobi method, we would have:

$$ \phi_{N+1}(x,y) = \frac{\phi_N(x+a,y) + \phi_N(x-a,y) + \phi_N(x,y+a) + \phi_N(x,y-a)}{4} + \frac{a^2\rho(x,y)}{4 \epsilon_0} $$

Note that now that externally imposed charge $\rho$ is not zero, one needs to now account for the lattice spacing correctly to calculate the correct potential to achieve the correct values of the source charges.

For implementing SOR, we would then follow the Gauss-Seidel approach for the updating of $\phi$ and also over-relax.

In our simulation, we will also now need an additional matrix rho to specify the positions of fixed charge densities in the simulation.

Exercise 4(a): Peform a simulation of a 2d "point source" consisting of a single pixel at position i,j = 50,50 with rho = 1. In the simulation, take $a = 1$ mm. (Note that our "point source" 2D simulation actually corresponds to the voltage around a thin, infinitely long wire.) You can use SOR if you want, and a reasonble value for $\omega$ is 1.9

In [ ]:
target_accuracy = 1e-4
M = 101
a = 1e-3
eps0 = 8.85e-12

# Our fixed charge input matrix
rho_fixed = np.zeros([M,M])
rho_fixed[M//2,M//2] = 1

### BEGIN SOLUTION
phi = np.zeros([M,M])
delta = phi.copy()

# A matrix to keep track of which points in our simulation 
# are fixed voltages. In this case, it is only the boundaries
# of the simulation
fixed = np.zeros([M,M])
fixed[[0,-1],:] = 1
fixed[:,[0,-1]] = 1

delta_max = 1
delta_max_list = []

t1 = time()
SOR_parameter = 1.9

while delta_max > target_accuracy:
    for i in range(M):
        for j in range(M):
            if not fixed[i,j]:
                delta[i,j] = (phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1])/4 \
                    + a**2/4/eps0*rho_fixed[i,j] - phi[i,j]
                phi[i,j] += delta[i,j]*1.9
    delta_max = np.max(np.abs(delta))
    delta_max_list.append(delta_max)
    print("N_iter %d delta_max %e\r" % (len(delta_max_list), delta_max), end='')

t2 = time()
print("\nTotal running time: %.2f min" % ((t2-t1)/60))
print("Code speed: %.1f iterations per second" %(len(delta_max_list)/(t2-t1)))
### END SOLUTION

answer_13_4a_1 = phi.copy()
answer_13_4a_2 = rho_fixed.copy()
In [ ]:
question = "answer_13_4a"
num = 2

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback
In [ ]:
# Notebook code
explore_image(answer_13_4a_1)

Calculating charge distributions

In all of the above simulations, even ones with external charge $\rho=0$, there was charge present in the resulting calculated physical problem. How do I know this? If there is no charge, there will be no change in potential. OK, but where is this charge?

With the exception of the last simulation, the charge in the simulation was induced charge that formed due to the voltage on the applied to the metals boundaries.

Even in the last simulation of exercise 4, where there is very clearly charge in the middle of the simulation where we placed the fixed charge, there is also charge induced on the boundary walls: the image charge induced by the electrical fields from the point charge we added in the middle.

How do we calculate this charge? Well, if we know $\phi(x,y)$, we can just use the Poisson equation:

$$ \rho(x,y) = \epsilon_0 \nabla^2 \phi(x,y) $$

where numerically, we use the following to calculate $\nabla^2 \phi$:

$$ \frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} \approx \frac{ \phi(x+a,y) + \phi(x-a,y) + \phi(x,y+a) + \phi(x,y-a) - 4 \phi(x,y)} {a^2} $$

A subtle point arise at the boundary: for the outer boundary points, don't have neighbors. For example, a point on the left of the simulation has no point at $x-a$.

What do I do? One could add special code for calculating the charge for these boundary conditions, but this would be a lot of cumbersome work.

A very simple and pragmatic way around this for the "metallic" boundary condition is just to add an extra grid point: make the boundaries two grid points thick. Since the induced charge in a (non-quantum-mechanical) metal is a perfect surface layer, we know that these "extra" boundary points must have zero charge. We can then just set them all to zero, and then use the Poisson equation above to calculate the charge density everywhere, just skipping the outermost pixels.

Let's implement this for the problem of exercise 4.

Exercise 5(a): Implement exercise 4(a) with an extra boundary grid point on all sides.

In [ ]:
target_accuracy = 1e-4
M = 103
a = 1e-3
eps0 = 8.85e-12

phi = np.zeros([M,M])
delta = phi.copy()

# A matrix to keep track of which points in our simulation 
# are fixed voltages. In this case, it is only the boundaries
# of the simulation
fixed = np.zeros([M,M])
# fixed[[___,___,___,___],:] = 1
# fixed[:,[___,___,___,___]] = 1
###
### Your code here
###



# The external fixed charge
rho_fixed = np.zeros([M,M])
# rho_fixed[___,___] = 1
###
### Your code here
###


# ...
###
### Your code here
###



answer_13_5a_1 = rho_fixed.copy()
answer_13_5a_2 = phi.copy()
In [ ]:
question = "answer_13_5a"
num = 2

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

Now that we have the calculate potential with the "extra boundary" point, we can calculate the charge density in our simulation, including the charge density induced on the boundary walls.

As we are working with 2D simulations of infinitely long objects, it is actually more convenient to work instead with the 2D surface charge density $\sigma$ (in units of C/m$^2$) instead of the 3D charge density $\rho$ (in units of C/m$^3$). For this, we can use the following equation:

$$ \sigma(x,y) = - \frac{\epsilon_0}{a} [ \phi(x+a,y) + \phi(x-a,y) + \phi(x,y+a) + \phi(x,y-a) - 4 \phi(x,y) ] $$

Exercise 5(b): Write a routine calculate_sigma() that calculate the surface charge density in a 2D Poission simulation. To keep things simple, you can set up your function to take no arguments and instead use variables defined in the global scope. It should return the calculated surface charge density. Use it to calculate the surface charge density for the simulation in question 4(a).

In [ ]:
def calculate_sigma():
    sigma_calc = np.zeros([M,M])
    ###
### Your code here
###

    return sigma_calc

answer_13_5b_1 = calculate_sigma()
In [ ]:
question = "answer_13_5b"
num = 1

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

You can use this code to make a plot of the calculated surface charge density. I've chosen a range of the colorscale such that it is not overwhelmed by the large charge density of the point source.

In [ ]:
plt.imshow(answer_13_5b_1, cmap="RdBu")
plt.colorbar()
plt.clim(vmax=-np.min(answer_13_5b_1)-2e-6)

You can also "explore" the image here and look at some line cuts. I would suggest a "gamma" setting of about 0.1 to make the charge on the edges visible, and turning "autoscale" on.

In [ ]:
explore_image(answer_13_5b_1,zname='$\sigma"')

What do you expect the total charge to be in your simulation? Use np.sum() to check your answer.

In [ ]:
np.sum(answer_13_5b_1)