Homework 3: Solving Laplace's Equation

Mathematical Background: Laplace's Equation

A wide variety of scientific phenomena can be described by this equation:

This is a second order partial differential equation. (DON'T PANIC!) It is referred to as "Laplace's equation" after its discoverer, the French mathematician Pierre-Simon Laplace (1749-1847). It has uses in the fields of astronomy, atmospheric physics, electromagnetism, and fluid mechanics and thus has wide application in electrical and mechanical engineering. Some engineers owe their jobs to this equation!

Solving Laplace's Equation

Purely mathematical solutions of Laplace's equation are rare. Approximate solutions on a computer (known as numerical solutions) are much more common. In this assignment, we will use the computer to find one such solution. We'll start with a square array of values phi[i][j] which is N by N for a given integer N. This is a grid which will contain the values of the original φ at the grid points (array cells).

Actually, we don't solve for every value of phi[][]. The phi[][] values on the edges of the grid are fixed to given values and not allowed to change. These are the boundary conditions.

Our solution of Laplace's equation finds values for the grid points inside the boundary. Our solution is iterative: We start with a poor solution and repeat a procedure called relaxation a given number of times that makes it a little better every time it's repeated.

As usual, we'll break the program down into several functions.

Part 1: initPhi()

initPhi() sets the initial values of all elements of the phi[][] array. This means the boundary as well as interior grid points. Its prototype is:

void initPhi(double phi[N][N])

N is the number of rows and columns in phi[][]. It is a #defined constant that will have been set before printPhi(). (You can use it everywhere you would use an int variable, except that you can't assign a value to it.)

The boundary conditions initPhi() "enforces" are simple: 1.0 along the top border (including corners) and 0.0 on all other borders and corners. Here's a table with specifics:

This sets values on the borders of the grid. initPhi() will set all interior values to 0.0.

Part 2: printPhi()

printPhi() prints out the contents of the phi[][] array that's passed to it. Its prototype is:

void printPhi(double phi[N][N])

The output of the function is obvious: One line per row i of the grid and each line contains the contents of all grid cells phi[i][j] for all possile j. Use " %4.2f" for the printf() calls.

Here's the result of calling printPhi() on the result of initPhi() for an N of 5:

1.00 1.00 1.00 1.00 1.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00

At this point, it would be a good idea for you to create a test program that does just this to make sure both functions are working.

Part 3: relax()

relax() performs a single relaxation step.

This will make phi[][] converge closer to the solution of Laplace's equation. It's a bit more involved than the first two functions, so here's the pseudocode:

void relax(
    double phi_k[N][N],   /* input -- the current iteration of phi[][] */
    double phi_kp1[N][N]) /* output -- the next iteration of phi[][] */
{
    /*
     * for i from 1 to N-2, inclusive,
     *     for j from 1 to N-2, inclusive,
     *         set phi_kp1[i][j] to the mean of its four nearest
     *          neighbors in phi_k: phi_k[i-1][j], phi_k[i][j-1],
     *          phi_k[i][j+1], and phi_k[i+1][j]
     */
}

Note the use of mnemonics in coverting math to C: phi_k is the C equivalent of the mathematical φ^(k) (the kth iteration of φ) and phi_kp1 is the C equivalent of φ^(k + 1) (the k + 1st iteration of φ). The kp1 stands for "k + 1".

Here's the result of calling printPhi() on the result of a single relax() for the above example:

1.00 1.00 1.00 1.00 1.00
0.00 0.25 0.25 0.25 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00

You can verify for yourself that each new interior cell is the mean of its four neigbors in the old array.

Laplace's equation is solved when a relax() iteration produces no change to any element of the array. We assume here that 10 iterations is sufficient. Is that correct? Is it too much?

Part 4: main()

The pseudocode for main() is:

int main(void)
{
    /*
     * declare two double arrays, phiCurrent[N][N] and phiNew[N][N] 
     * initialize phiCurrent and phiNew with initPhi calls
     * print "initial phi" message
     * use printPhi to print phiCurrent
     * for n (the iteration count) from 1 to 10:
     *     print a newline (for clarity)
     *     call relax to relax phiCurrent into phiNew
     *     print "after n iterations" message
     *     use printPhi to print phiNew
     *     copy all elements of phiNew to phiCurrent
     */
    return 0;
}

This function puts all the others together and shows the results of multiple repeated relaxations of the original phi[][].

Because of the generality of Laplace's equation, this same code could be used to simulate (in 2D) the electrical potential underneath a charged plate or the steady-state distribution of temperature below a heated surface.

It's not required for the homework, but here are a few interesting experiments to perform:

• The value of N is arbitrary. Try setting it to, say, 10.
• Change the boundary conditions by forcing the bottom row to be -1.0. (This might then represent electrical potential in the region between two charged plates: one positive and one negative.)
• We assumed that 10 iterations were enough to stabilize phi[][] (i.e. it wouldn't change after that many iterations). Is that true? If so, is it still true for larger N?
• We assumed that we had a square region. What if it were rectangular? What if it were circular?