Welcome to the fourth notebook of Relax and hold steady: elliptic problems, the fifth module of "Practical Numerical Methods with Python".
In the previous notebooks, we examined traditional relaxation methods (iterative methods). Modern engineering problems involve millions of grid points and as the number of grid points grows, the convergence rate of iterative methods starts to suffer. Enter multigrid.
Consider the discretized 1D Poisson equation:
\begin{equation} \frac{p_{i+1}-2p_i+p_{i-1}}{\Delta x^2}=b_i ,\ \ \ \ \ 0 \le i \le N_x-1 \end{equation}where $N_x$ is the number of grid points. Since the indices of the grid points start from $0$, the index of the last point is $N_x-1$. We start with Dirichlet boundary conditions: $$p_0=p_{N_x-1}= 0$$
If $p{^n_i}$ and $p{^{exact}_i}$ denote the approximated solution at the n-th iteration in a relaxation method and the exact solution, respectively, we can define the error at this iteration as
\begin{equation} e{^n_i} = p{^{exact}_i} - p{^n_i} ,\ \ \ \ \ 0 \le i \le N_x-1 \end{equation}That is, the error is the difference between the exact and the estimated solution.
When we specify an initial guess, $p{^0_i}$, the initial error will be given by
$$e{^0_i} = p{^{exact}_i} - p{^0_i}$$When we solve an elliptic PDE by any iterative method we try to eliminate this error.
If we were to apply a Fourier transform to the errors $e^n_i$, we would find that the errors are actually composed of many different single-wavenumber waves. This yields an interesting question:
When we eliminate errors using a relaxation method, do these wavenumbers make an impact on the convergence rate?
We can answer this with some simple tests.
For simplicity, consider the 1D Laplace equation with Dirichlet boundaries
$$p_0=p_{N_x-1}= 0$$The exact solution is $p{^{exact}_i}=0$ and Equation $()$ now simplifies to
\begin{equation} e{^n_i}=-p{^n_i} \end{equation}We want to examine error convergence as a function of wavenumber so we need a function that is dependent on wavenumber. We can define $p^0_i$ as follows:
\begin{equation} p{^0_i} = \sin{\left(\frac{ik\pi}{N_x-1} \right)} ,\ \ \ \ \ 0 \le i \le N_x-1 \end{equation}where $k$ is the wavenumber. Let's take a look at our initial condition for a few different values of $k$.
import numpy
from matplotlib import pyplot, rcParams
%matplotlib inline
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
k_values = [1, 16, 31, 48, 63]
nx = 65
i = numpy.arange(nx-1)
pyplot.figure(figsize=(10,6))
for j,k in enumerate(k_values):
pyplot.subplot(2,3,j+1)
p = numpy.sin(i * k * numpy.pi / (nx-1))
pyplot.plot(numpy.linspace(0,1,nx-1), p)
pyplot.title("$k = $"+str(k));
pyplot.tight_layout();
Those are certainly different! Now we want to investigate how the error decreases for each of these different initial conditions!
We have a few functions to help streamline the error analysis.
def laplace1d_IC(nx, k):
'''
Generates initial guess for Laplace 1D eq. under a given number
of grid points (nx) and wavenumber k within domain [0,1]x[0,1]
Parameters:
----------
nx: int, number of grid points in x direction
k: float, wavenumber
Returns:
-------
p: 1D array of float, initial guess of the unknowns
b: 1D array of float, 0th-order derivative term in Poisson eq.
x: 1D array of float, linspace coordinates in x
dx: float, grid spacing in x
'''
dx = 1.0/(nx-1)
x = numpy.linspace(0,1,nx)
##initial conditions
i = numpy.arange(0, nx)
p = numpy.sin(i * k * numpy.pi / (nx-1))
b = numpy.zeros(nx)
return p, b, x, dx
To solve our example problem we'll once again use Gauss-Seidel with weave.blitz
. You will notice that in contrast to the previous notebook, the Poisson and Laplace solving functions perform only a single iteration per function call. The reason for this will become clear as we delve deeper into multigrid methods.
from scipy import weave
def poisson1d_GS_SingleItr(nx, dx, p, b):
'''
Gauss-Seidel method for 1D Poisson eq. with Dirichlet BCs at both
ends. Only a single iteration is executed. **blitz** is used.
Parameters:
----------
nx: int, number of grid points in x direction
dx: float, grid spacing in x
p: 1D array of float, approximated soln. in last iteration
b: 1D array of float, 0th-order derivative term in Poisson eq.
Returns:
-------
p: 1D array of float, approximated soln. in current iteration
'''
expr = "p[1:-1] = 0.5 * (p[:-2] + p[2:] - dx * dx * b[1:-1])"
weave.blitz(expr, check_size=0)
return p
We used a relative $L_2$-norm to calculate the error in the previous notebooks, but that doesn't work very well when the exact solution is a line of zeros. Instead of $L_2$, we can use the root-mean-square of the difference of two solutions to compare error between different grid sizes in a meaningful way.
def RMS(p):
'''
Return the root mean square of p.
Parameters:
----------
p: array
Returns:
-------
Root mean square of p
'''
return numpy.sqrt(numpy.sum(p**2) / p.size)
Time to get started! We will examine the error associated with the wavenumbers
$$k = [1, 16, 31, 48, 63 ]$$for two different grid sizes,
$$nx = [65, 33]$$Note: We are storing the error for each wavenumber-grid combination in a dictionary. A dictionary is an unordered collection with keys and values. Each key in our dictionary is a tuple comprised of the grid-size and the wavenumber. The corresponding value for a given key (in the code below) is a NumPy array containing the RMS error for each iteration.
# define different Np
nx = [65, 33]
# define wavenumbers
k = [1, 16, 31, 48, 63]
# number of iterations we will run
Nitr = 100
# initialize a space to store root mean square errors
err = {}
# start iteration
for nxi in nx:
for ki in k:
key = (nxi, ki)
err[key] = numpy.empty(Nitr+1)
p, b, x, dx = laplace1d_IC(nxi, ki)
err[key][0] = RMS(p)
for itr in range(1, Nitr+1):
p = poisson1d_GS_SingleItr(nxi, dx, p, b)
err[key][itr] = RMS(p)
# plot errors in each case and at each iteration
pyplot.figure(figsize=(9, 4.1))
for n, nxi in enumerate(nx):
pyplot.subplot(1, 2, n+1)
for ki in k:
pyplot.semilogy(numpy.array(range(Nitr+1)), err[(nxi, ki)],
lw=2, label='k='+str(ki))
pyplot.legend(loc=3, fontsize='small', ncol=3, mode='expand')
pyplot.xlabel('Iteration')
pyplot.ylabel('RMS Error')
pyplot.title('Nx = '+str(nxi))
pyplot.ylim((1e-8, 1))
pyplot.tight_layout();
Why are there only three lines in the $N_x=33$ plot? In fact, there are five lines, but two of them are hiding. What's happening?
Any grid has a limited resolution and if we try to resolve a rapidly oscillating function on a coarser grid, we encounter a phenomenon called aliasing. The high wavenumber error can't be fully resolved in the domain, so it wraps around.
In the example above, the error for wavenumbers $k = [48, 63]$ is exactly the same as for wavenumbers $k=[16, 1]$, respectively. That's a visual observation, but we also check the numbers using the error dictionary.
numpy.allclose(err[(33,1)],err[(33,63)])
numpy.allclose(err[(33,16)],err[(33,48)])
The implication here is that no error can have a wavenumber greater than $N_x-1$.
Another important observation in the plots above is the rate of convergence.
Compare the rate of convergence for wavenumbers $k=[1,16,31]$ in both plots above. What do you see? The error converges more quickly for these wavenumbers in the $N_x=33$ plot.
But what about the wavenumbers $k = [48, 63]$? They clearly converge faster in the $N_x=65$ domain.
These observations are at the heart of multigrid.
If we assume that a solution exists to the Poisson equation
\begin{equation} \nabla^2 p = b \end{equation}then we know that at each point in the domain, there will be some error between the calculated value, $p^n_i$ and the exact solution $p^{exact}_i$. We may not know what the exact solution is, but we know it's out there. Earlier, we defined the error at any point $i$ as
\begin{equation} e^n_i = p^{exact}_i - p^n_i \end{equation}Note: We are talking about error at a specific point, not a measure of error across the entire domain.
What if we recast the Poisson equation in terms of a not-perfectly-relaxed $p^n$?
\begin{equation} \nabla^2 p^n \approx b \end{equation}Equation $(7)$ is an approximation because $p^n \neq p$. To "fix" the equation, we need to add an extra term to account for the difference in the Poisson equation -- that extra term is called the residual. We can write out the modified Poisson equation like this:
\begin{equation} r^n + \nabla^2 p^n = b \end{equation}Rearranging and discretizing with 2nd-order central difference yields
\begin{equation} r^n_i = b_i - \frac{p{^n_{i+1}}-2p{^n_i}+p{^n_{i-1}}}{\Delta x^2} \end{equation}If we combine Equation (6) with Equation (9), it gives us the following result:
\begin{equation} r^n_i = \frac{e{^n_{i+1}}-2e{^n_i}+e{^n_{i-1}}}{\Delta x^2} \end{equation}Still with us? What sort of system does Equation $(10)$ look like? The residual can be used as a source term to relax the error itself. Then, according to Equation $(6)$, we can add this relaxed error back to the corresponding $p_i$ term to increase the accuracy of the solution.
def residual(dx, pn, b, r):
'''
Calculate the residual for the 1D Poisson equation.
Parameters:
----------
pn: 1D array, approximated solution at a certain iteration n
b: 1D array, the b(x) in the Poisson eq.
Return:
----------
The residual r
'''
# r[0] = 0
r[1:-1] = b[1:-1] - (pn[:-2] - 2 * pn[1:-1] + pn[2:]) / dx**2
# r[-1] = 0
return r
One of the important bookkeeping operations in multigrid is the transfer of data from one grid to another.
Transfering data from a fine grid to a coarse grid is called restriction. There are many ways to implement restriction. Here we introduce two restriction operators. The first is called injection, the other is full weighting.
Injection
Injection is the simplest method of restriction. Where two grids have matching coordinates, it copies over the corresponding data. In the locations where one grid has points and the other doesn't, it simply ignores that data.
If we have a coarse grid with half as many points as a fine grid, we can write out the injection operator as follows:
\begin{equation} v_{i, c} = v_{2i, f} \end{equation}With Dirichlet boundary conditions, the boundaries are mapped one-to-one. Neumann boundary conditions require a little extra consideration as we'll see later.
Full weighting
In contrast to injection, full weighting considers the contributions of all points when mapping from the fine grid to the coarse grid.
The Dirichlet boundary conditions are again mapped one-to-one.
Different restriction operators can impact the rate of convergence. We implement full weighting in this example but you should also try the same problem using injection and compare the performance of the two methods.
The following code block shows a function for full weighting. Note that we use stepped slicing here.
The syntax for stepped slicing is
array[start:end:stepsize]
For example, x[1:10:2]
would return the values of x[1]
, x[3]
, x[5]
, x[7]
, and x[9]
.
def full_weighting_1d(vF, vC):
'''
Transfer a vector on a fine grid to a coarse grid with full weighting
. The number of elements (not points) of the coarse grid is
half of that of the fine grid.
Parameters:
----------
vF: 1D numpy array, the vector on the fine grid
vC: 1D numpy array, the vector on the coarse grid,
size(vC) = (size(vF) + 1) / 2
Output: vC
'''
vC[0] = vF[0]
vC[1:-1] = 0.25 * (vF[1:-3:2] + 2. * vF[2:-2:2] + vF[3:-1:2])
vC[-1] = vF[-1]
return vC
After we relax the residual equation on the coarse grid, we need an operation that can transfer the relaxed errors back to the fine grid.
The simplest way to move from the coarse grid to the fine grid is linear interpolation. The values of the points common to both grids are set equal while the values at the remaining points on the fine grid are calculated via linear interpolation, as shown below.
\begin{equation} \begin{array}{c} v_{2i,f}=v_{i, c} ,\ \ \ \ \ 0 \le i \le N_{x,c}-1 \\ v_{2i+1,f}=\frac{1}{2}\left(v_{i,c}+v_{i+1,c}\right) ,\ \ \ \ \ 0 \le i \le N_{x,c}-2 \end{array} \end{equation}def interpolation_1d(vC, vF):
'''
Transfer a vector on a coarse grid to a fine grid by linear
interpolation. The number of elements (not points) of the coarse
grid is a half of that of the fine grid.
Parameters:
----------
vC: 1D numpy array, the vector on the coarse grid,
vF: 1D numpy array, the vector on the fine grid
size(vF) = size(vC) * 2 - 1
Output: vF
'''
vF[::2] = vC[:]
vF[1:-1:2] = 0.5 * (vC[:-1] + vC[1:])
return vF
Aliasing http://nbviewer.ipython.org/github/ketch/teaching-numerics-with-notebooks/blob/master/Aliasing.ipynb
William L. Griggs, Van Emden Henson, and Steve F. McCormick, A Multigrid Tutorial, SIAM, Philadephia, 2000
U. Trottenberg, C. W. Oosterlee, Anton Schüller, Multigrid, Academic press, 2000
Painless Conjugate Gradient
from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())