Relaxation Methods

Modules – Partial Differential Equations

Last edited: 19th of November 2017

Relaxation methods involves solving the system of equations that arises from a finite difference approximation iteratively until a solution is found [1]. We will introduce the Jacobi method, the Gauss-Seidel method and succesive overrelaxation (SOR). The methods can be used to solve an elliptic equation

$$\mathcal L \psi = \rho,$$

where $\mathcal L$ is an elliptical operator and $\rho$ is a source term.

The methods are simple to understand through a concrete example. We will be solving the two-dimensional Poisson equation,

\begin{equation} \nabla^2 \psi(x, t) = \frac{\partial^2 \psi(x, y)}{\partial x^2} + \frac{\partial^2 \psi(x, y)}{\partial y^2} = \rho(x, y), \label{eq:poisson} \end{equation}

for a square domain $\Omega = (0,L)\times (0, L)$, $L\in \mathbb R_+$. Without loss of generality we let $L=1$.

The Poisson Equation in Physics

The Poisson and Laplace $(\rho = 0)$ equations are widely used in physics. In electrostatics, the electric potential is described by $\nabla^2 V\left(\vec r\right) = -\rho \left(\vec r\right)/\epsilon_0,$ where $V$ is the electric potential and $\epsilon$ is the electric permittivity (this is the differential form of Gauss' law for the electric potential. See e.g. [2]). The steady state of the diffusion equation becomes $u_t = \nabla^2 u - \rho = 0$. For a steady incompressible flow in two dimensions, the stream function satisfies $\nabla^2 \psi = -\zeta$, where $\zeta=\nabla\times \vec V$ is the vorticity and the velocity field is given by $\vec V = (u, v) = (\psi_y, -\psi_x)$ (See e.g. [3]).

The solution is unique when the boundary conditions for the entire system is well defined. This is discussed in more detail in the appendix below.

Jacobi Method

Consider the equation

\begin{equation} \frac{\partial \psi(x, y, t)}{\partial t}=\frac{\partial^2 \psi(x, y, t)}{\partial x^2} + \frac{\partial^2 \psi(x, y, t)}{\partial y^2} - \rho(x, y), \label{eq:diffusion} \end{equation}

whose steady state is the Poisson equation \eqref{eq:poisson}. Let $x$, $y$ be discretized to some grid,

$$x_i = x_0 + i\Delta x, \: (\:i=0, 1, ...,N)$$

$$y_i = y_0 + j\Delta y, \: (\:j=0, 1, ..., N)$$ and let, $$t_n = n\,\delta t, \: (\:n = 0, 1, ..., n).$$

Since we are considering the square domain $\Omega$, we have $x_0 = y_0 = 0$ and $x_n = y_m = 1$. To simplify the notation we introduce $u(x_i, y_j, t_n)\equiv u_{i, j}^n$ and $\rho(x_i, y_j) \equiv \rho_{i, j}$.

If we use central differencing in space and and foreward differencing in time (FTCS) on equation \eqref{eq:diffusion}, we get

$$\frac{\psi^{n+1}_{i, j} - \psi^n_{i, j}}{\Delta t} = \frac{\psi^n_{i+1, j}-2\psi^n_{i, j}+\psi^n_{i-1, j}}{(\Delta x)^2} + \frac{\psi^n_{i, j+1}-2\psi^n_{i, j}+\psi^n_{i, j-1}}{(\Delta y)^2} - \rho_{i, j}.$$

For convinience we let $\Delta \equiv \Delta x=\Delta y$. We then obtain

\begin{equation} \psi^{n+1}_{i, j} = \psi^n_{i, j} + \frac{\Delta t}{\Delta^2} \left(\psi_{i+1, j}^n + \psi_{i-1, j}^n + \psi_{i, j+1}^n + \psi_{i, j-1}^n - 4\psi_{i, j}^n\right) - \Delta t \,\rho_{i, j}. \label{eq:poisson_differenced} \end{equation}

It can be shown that the FTCS differencing for the Poisson equation in two dimensions is stable only if $\Delta t/\Delta^2 \leq 1/4$ (see Appendix). If we set $\Delta t/\Delta^2 = 1/4$ we get

\begin{equation} \psi^{n+1}_{i, j} = \frac{1}{4} \left(\psi_{i+1, j}^n + \psi_{i-1, j}^n + \psi_{i, j+1}^n + \psi_{i, j-1}^n\right) - \frac{\Delta ^2}{4} \,\rho_{i, j}. \label{eq:jacobi_method} \end{equation}

The algorithm thus consists of choosing the value at $(x_i,y_j)$ as the average of the neighboring points plus some source term $\rho$.

We can use any initial guess. The iteration is repeated until a convergence criterion is met. We will be using the total relative change between two iterations as a stopping condition. The final solution is the steady state solution of equation \eqref{eq:diffusion}, which is the Poisson equation \eqref{eq:poisson}.


Consider a square domain without any sources. The Poisson equation reduces to the Laplace equation $(\rho = 0)$ inside the domain. The boundary conditions are indicated in figure 1 below. We choose $L=1$, $a=0$ and $b = 1$ (with the necessary dimensions) without loss of generality.

Example - Square setup Figure 1: Setup used in the example. The value of $\psi$ along the solid lines are $a$ and $b$ as indicated. $\psi$ vary evenly between $a$ and $b$ along the dotted lines.

This setup can for example model a conductor (the upper left part, $a$) and a conductor with a uniformly distributed charge (the lower right part, $b$). The electric field is in turn given by $\vec E = \nabla \psi$. The setup can also model the flow of an incompressible and irrotational fluid through a square cavity with uniformal inflow and outflow. The inflow and outflow being the dottet lines, respectivaly the right and bottom. The velocity field is given by $\vec V = (\partial_y \psi, - \partial_x \psi)$. We refer you to our seperate notebooks on stream functions and electromagnetism for more information.

To solve the Laplace equation for this setup we use the Jacobi method. We define the following functions.

In [92]:
import numpy as np

def bc(n=32):
    """ Represents the boundry conditions of psi as in the setupin figure 1
        in a 2D matrix.
        n :     int. Number of discretisations in the x- and y-
    int array, size (n, n). The boundry conditions.
    psi = np.zeros((n, n))
    width = int(n/3)

    # Boundary conditions
    psi[0, 2*width:] = 1  # Lower left part of b
    psi[:width, -1] = 1   # upper right part of b
    psi[0, width:2*width] = np.linspace(0, 1, width)  # Inflow
    psi[width:2*width, -1] = np.linspace(1, 0, width) # Outflow
    # The rest are already zero
    return psi

def jacobi(n=32, maxit=1000, tol=1e-10):
    """ This function uses Jacobi's method to solve the Poisson
    equation for the square domain and boundary conditions as
    indicated in the figure.
        n :     int. Number of discretisations in the x- and y-
        maxit : int. Maximum number of iterations.
        tol :   float. Tolerance in convergence criterion.
                If the total relative change between two
                iterations are less that 'tol', we say that
                the method has converged.
    float array, size(n, n). The estimated solution. 
    psi = bc(n)
    psi_temp = psi.copy()
    itnum = 1             # Iteration counter
    change = tol + 1      # Total relative quadratic change between each iteration
    for itnum in range(maxit):
        # Let each point be the mean of its nearest neighbors
        psi_temp[1:-1, 1:-1] = .25*(psi[2:, 1:-1] + psi[:-2, 1:-1] +
                                    psi[1:-1, 2:] + psi[1:-1, :-2])
        # Update the convergence check for every 50 iteration
        if itnum % 50 == 0:
            if np.sum(np.abs(psi_temp - psi))/np.sum(psi) < tol:
                print("Iterations: ", itnum + 1)
        psi = psi_temp.copy()
    if itnum == maxit - 1:
        print("Maximum number of iterations (%i) reached!"%(maxit))
    return psi

Let's perform one example.

In [93]:
n = 64
tol = 1e-5
psi = jacobi(n, maxit=10000, tol=1e-5)
Iterations:  3551

We then visualize the result for the two physical interpretations mentioned above.

In [94]:
# Needed package
import matplotlib.pyplot as plt
%matplotlib inline

def plot(psi, n):
    # Compute the gradient. v=psi_yy, u=psi_xx
    v, u = np.gradient(psi)

    x = np.linspace(0, 1, n) # x-axis
    vmag = np.sqrt(u**2 + v**2)

    # Create figure
    fig, axarr = plt.subplots(1, 2, figsize=(10, 6), dpi=300)

    # Axes 1: Electric potential and field
    im1 = axarr[0].contourf(x, x, psi, 20)
    axarr[0].streamplot(x, x, -u, -v, color="k")
    axarr[0].set_title("Electric potential and field")
    fig.colorbar(im1, orientation='horizontal', ax=axarr[0],
                  label=r"Electric potential, $V/V_0$")

    # Axes 2: Velocity field and strength
    # Flow of an incompressible and irrotational fluid
    im2 = axarr[1].contourf(x, x, vmag/np.max(vmag), 20)
    axarr[1].streamplot(x, x, v, -u, color="k")
    axarr[1].set_title("Velocity field and strength")
    fig.colorbar(im2, orientation='horizontal', ax=axarr[1],
               label=r"Velocity magnitude, $v/v_0$")
    axarr[1].set_xlim([0, 1])
    axarr[1].set_ylim([0, 1])
In [95]:
plot(psi, n)