Relaxation Methods¶
Modules – Partial Differential Equations¶
Last edited: February 7th 2018
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(\mathbf r\right) = -\rho \left(\mathbf 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 \mathbf V$ is the vorticity and the velocity field is given by $\mathbf 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}.
Example¶
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.
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 $\mathbf 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 $\mathbf 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.
import numpy as np
def bc(n=32):
""" Represents the boundry conditions of psi as in the setupin figure 1
in a 2D matrix.
Parameters:
n : int. Number of discretisations in the x- and y-
direction
Returns:
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.
Parameters:
n : int. Number of discretisations in the x- and y-
direction
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.
Returns:
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)
break
psi = psi_temp.copy()
if itnum == maxit - 1:
print("Maximum number of iterations (%i) reached!"%(maxit))
return psi
Let's perform one example.
n = 64
tol = 1e-5
psi = jacobi(n, maxit=10000, tol=1e-5)
We then visualize the result for the two physical interpretations mentioned above.
# Needed package
import matplotlib.pyplot as plt
plt.style.use('bmh')
%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$")
axarr[0].set_xlabel("$x/L$")
axarr[0].set_ylabel("$y/L$")
# 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_xlabel("$x/L$")
axarr[1].set_ylabel("$y/L$")
axarr[1].set_xlim([0, 1])
axarr[1].set_ylim([0, 1])
plot(psi, n)