Simple example of relaxation method used to solve Laplace's equation

From Computational Physics by Mark Newman, Section 9.1

Program laplace.py downloaded directly from Computational Physics - Programs and data - University of Michigan

This is a 2-D problem with $V = 1$ on one side, and $V=0$ on other three sides. Griffiths, Section 3.3.1 solves a very similar problem analytically using separation of variables and Fourier series; in Griffiths, the direction perpendicular to the $V=!$ side extends to infinity. A simple extension of Newman's code that enlarges the grid in one dimension leads to results that agree with Griffiths.

M.L.'s notebook laplace.ipynb does the same calculation as this notebook, but without Newman's explicit looping.

Instead of looping, the array phi is shifted up one grid step, shifted down one grid step, shifted left one grid step, and shifted down one grid step, and the results are summed and then divided by 4.

`for i in range(20000): phi_new = (np.roll(phi_old, 1, axis=0) + np.roll(phi_old, -1, axis=0) \

           + np.roll(phi_old, 1, axis=1) + np.roll(phi_old, -1, axis=1))/4
phi_new[:, 0] = phi_new[:,m] = 0  # Reset boundaries
phi_new[n,:] = 0                  # Reset boundary
phi_new[0,:] = 1                  # Reset boundary`

This is MUCH faster than Newman's looping procedure!

from numpy import empty,zeros,max
from pylab import imshow,gray,show

# Constants
M = 100         # Grid squares on a side
V = 1.0         # Voltage at top wall
target = 1e-6   # Target accuracy

# Create arrays to hold potential values
phi = zeros([M+1,M+1],float)
phi[0,:] = V
phiprime = empty([M+1,M+1],float)

# Main loop
delta = 1.0
count = 0
while delta>target:
    count += 1
    # Calculate new values of the potential
    for i in range(M+1):
        for j in range(M+1):
            if i==0 or i==M or j==0 or j==M:
                phiprime[i,j] = phi[i,j]
            else:
                phiprime[i,j] = (phi[i+1,j] + phi[i-1,j] \
                                 + phi[i,j+1] + phi[i,j-1])/4

    # Calculate maximum difference from old values
    delta = max(abs(phi-phiprime))

    # Swap the two arrays around
    phi,phiprime = phiprime,phi

# Make a plot
imshow(phi)
gray()
show()