Solution of 2-D Poisson equation¶

We consider the Poisson equation $$ -(u_{xx} + u_{yy}) = f, \quad\textrm{in}\quad \Omega = [0,1] \times [0,1] $$ with boundary condition $$ u = 0, \quad\textrm{on}\quad \partial\Omega $$ We will take $$ f = 1 $$ Make a uniform grid with $n$ points in each direction, with spacing $$ h = \frac{1}{n-1} $$ The grid points are $$ x_i = ih, \quad 0 \le i \le n-1 $$ $$ y_j = jh, \quad 0 \le j \le n-1 $$ The finite difference approximation at $(i,j)$ is $$ -\frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{\Delta x^2} - \frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{\Delta y^2} = f_{i,j}, \qquad 1 \le i \le n-2, \quad 1 \le j \le n-2 $$

In [9]:

%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import numpy as np
from matplotlib import pyplot as plt

The next function computers the $L^2$ norm of the residual which is used to test for convergence.

In [10]:

def residual(h,u,f):
    n = u.shape[0]
    res = 0.0
    for i in range(1,n-1):
        for j in range(1,n-1):
            uxx = (u[i-1,j]-2*u[i,j]+u[i+1,j])/h**2
            uyy = (u[i,j-1]-2*u[i,j]+u[i,j+1])/h**2
            res += (uxx + uyy + f[i,j])**2
    return np.sqrt(res)

Next, define some problem parameters.

In [11]:

xmin,xmax = 0.0,1.0
n = 21
h = (xmax - xmin)/(n-1)
x = np.linspace(xmin,xmax,n)
X,Y = np.meshgrid(x,x,indexing='ij')
plt.plot(X,Y,'.')
plt.title('Grid points')
plt.axis('equal')
f    = np.ones((n,n))
TOL   = 1.0e-6
itmax = 2000

Gauss-Jacobi¶

In [12]:

uold = np.zeros((n,n))
unew = np.zeros((n,n))
res0 = residual(h,unew,f)

for it in range(itmax):
    for i in range(1,n-1):
        for j in range(1,n-1):
            unew[i,j] = 0.25*(h**2 * f[i,j]
                              + (uold[i-1,j] + uold[i+1,j] + uold[i,j-1] + uold[i,j+1]))
    uold[:,:] = unew
    res = residual(h,unew,f)
    if res < TOL * res0:
        break

print("Number of iterations = %d" % it)
plt.contourf(X,Y,unew,cmap=plt.cm.jet)
plt.axis('equal')
plt.colorbar();

Number of iterations = 1102

Gauss-Seidel¶

In [13]:

u    = np.zeros((n,n))
res0 = residual(h,u,f)

for it in range(itmax):
    for i in range(1,n-1):
        for j in range(1,n-1):
            u[i,j] = 0.25*(h**2 * f[i,j]
                              + (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1]))
    res = residual(h,u,f)
    if res < TOL * res0:
        break

print("Number of iterations = %d" % it)
plt.contourf(X,Y,u,cmap=plt.cm.jet)
plt.axis('equal')
plt.colorbar();

Number of iterations = 552

SOR¶

In [14]:

omega= 2.0/(1 + np.sin(np.pi*h))
u    = np.zeros((n,n))
res0 = residual(h,u,f)

for it in range(itmax):
    for i in range(1,n-1):
        for j in range(1,n-1):
            z = 0.25*(h**2 * f[i,j]
                       + (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1]))
            u[i,j] = (1.0-omega)*u[i,j] + omega*z
    res = residual(h,u,f)
    if res < TOL * res0:
        break

print("Number of iterations = %d" % it)
plt.contourf(X,Y,u,cmap=plt.cm.jet)
plt.axis('equal')
plt.colorbar();

Number of iterations = 58