Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license © 2014 L.A. Barba, C.D. Cooper, G.F. Forsyth. Based on CFD Python, © 2013 L.A. Barba, also under CC-BY license.¶

Relax and hold steady¶

Welcome to the second notebook of "Relax and hold steady: elliptic problems", Module 5 of the course "Practical Numerical Methods with Python". Are you relaxed yet?

In the previous notebook, you learned to use Jacobi iterations to solve Laplace's equation. The iterations relax the solution from an initial guess to the final, steady-state solution. You also saw again that the way we treat boundary conditions can influence our solution. Using a first-order approximation of the Neumann boundary messed up our spatial convergence in the whole domain! (We expected second-order spatial convergence from the central difference scheme, but we got closer to first order.) This was easily fixed by using a second-order scheme for the Neumann boundary. It's always good to check that you get the expected order of convergence.

A word of warning: in this course module, we will introduce a different use of the word "convergence". Before, we used it to refer to the decay of the truncation errors (in space and time) with a decrease in the grid spacing ($\Delta x$ and $\Delta t$). Now, we also have a relaxation scheme, and we use the word convergence to mean that the iterative solution approaches the exact solution of the linear system. Sometimes, this is called algebraic convergence. We'll concern ourselves with this in the next lesson. But first, let's play with Poisson.

Poisson equation¶

The Poisson equation has a forcing function that drives the solution to its steady state. Unlike the Laplace equation, Poisson's equation involves imposed values inside the field (a.k.a., sources):

$$\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b$$

In discretized form, this looks almost the same as the Laplace Equation, except for the source term on the right-hand side:

$$\frac{p_{i+1,j}^{k}-2p_{i,j}^{k}+p_{i-1,j}^{k}}{\Delta x^2}+\frac{p_{i,j+1}^{k}-2 p_{i,j}^{k}+p_{i,j-1}^{k}}{\Delta y^2}=b_{i,j}^{k}$$

As before, we rearrange this to obtain an equation for $p$ at point $i,j$, based on its neighbors:

$$p_{i,j}^{k+1}=\frac{(p_{i+1,j}^{k}+p_{i-1,j}^{k})\Delta y^2+(p_{i,j+1}^{k}+p_{i,j-1}^{k})\Delta x^2-b_{i,j}^{k}\Delta x^2\Delta y^2}{2(\Delta x^2+\Delta y^2)}$$

It's slightly more complicated than the Laplace equation, but nothing we can't handle.

An example problem¶

Let's consider the following Poisson equation:

$$\begin{equation} \nabla^2 p = -2\left(\frac{\pi}{2}\right)^2\sin\left( \frac{\pi x}{L_x} \right) \cos\left(\frac{\pi y}{L_y}\right) \end{equation}$$

in the domain

$$\left\lbrace \begin{align*} 0 &\leq x\leq 1 \\ -0.5 &\leq y \leq 0.5 \end{align*} \right.$$

where $L_x = L_y = 1$ and with Dirichlet boundary conditions

$$p=0 \text{ at } \left\lbrace \begin{align*} x&=0\\ y&=0\\ y&=-0.5\\ y&=0.5 \end{align*} \right.$$

To solve this equation, we assume an initial state of $p=0$ everywhere, apply the boundary conditions and then iteratively relax the system until we converge on a solution.

To start, let's import the libraries and set up our spatial mesh.

The Jacobi iterations need an exit condition, based on some norm of the difference between two consecutive iterations. We can use the same relative L2-norm that we wrote for the Laplace exit condition, so we saved the function into a helper Python file (helper.py) for easy importing.

Now, what value to choose for the exit condition? We saw in the previous notebook that with an exit tolerance of $10^{-8}$, we could converge well for the different grids we tried, and observe second-order spatial convergence (with the second-order Neumann BC). We speculated in the end that we might be able to use a less stringent exit tolerance, since the spatial error was a lot larger (around $0.0002$ for the finer grid). Here, we'll try with $2\times 10^{-7}$. Go ahead and try with different values and see what you get!

It's time to write the function to solve the Poisson equation. Notice that all of the boundaries in this problem are Dirichlet boundaries, so no BC updates required!

There's also one extra piece we're adding in here. To later examine the convergence of the iterative process, we will save the L2-norm of the difference between successive solutions. A plot of this quantity with respect to the iteration number will be an indication of how fast the relaxation scheme is converging.

We can use the plot_3d function we wrote in the previous notebook to explore the field $p$, before and after the relaxation. We saved this plotting function into the helper Python file, so we can re-use it here.

Now we initialize all of the problem variables and plot!

That looks suitably boring. Zeros everywhere and boundaries held at zero. If this were a Laplace problem we would already be done!

But the Poisson problem has a source term that will evolve this zero initial guess to something different. Let's run our relaxation scheme and see what effect the forcing function has on p.

It took 3,125 iterations to converge to the exit criterion (that's quite a lot, don't you think? Let's now take a look at a plot of the final field:

Something has definitely happened. That looks good, but what about the error? This problem has the following analytical solution:

$$\begin{equation} p(x,y) = \sin{\left(\frac{x\pi}{L_x} \right)}\cos{\left(\frac{y\pi}{L_y} \right)} \end{equation}$$

Time to compare the calculated solution to the analytical one. Let's do that.

That seems small enough. Of course, each application problem can have different accuracy requirements.

Algebraic convergence¶

Remember that we saved the L2-norm of the difference between two consecutive iterations. The purpose of that was to look at how the relaxation scheme converges, in algebraic sense: with consecutive solutions getting closer and closer to each other. Let's use a line plot for this.

It looks like in the beginning, iterations started converging pretty fast, but they quickly adopted a slower rate. As we saw before, it took more than 3,000 iterations to get to our target difference between two consecutive solutions (in L2-norm). That is a lot of iterations, and we would really like to relax faster! No worries, we'll learn to do that in the next notebook.

Spatial convergence¶

For a sanity check, let's make sure the solution is achieving the expected second-order convergence in space.

That looks pretty much second order! Remember that the boundary conditions can adversely affect convergence, but Dirichlet boundaries are "exact" and will never impact your convergence.

Final word¶

We have used the difference between two consecutive solutions in the iterative process as a way to indicate convergence. However, this is not in general the best idea. For some problems and some iterative methods, you could experience iterates stagnating but the solution not converging.

Convergence of an iterative solution of a system $A \mathbf{x} = \mathbf{b}$ means that:

$$\begin{equation} \lim_{k \rightarrow \infty} \mathbf{x}^k = \mathbf{x} \end{equation}$$

The error in the solution is actually $\mathbf{x}-\mathbf{x}^k$, but we're looking at $\mathbf{x}^{k+1}-\mathbf{x}^k$ for our exit criterion. They are not the same thing and the second could tend to zero (or machine precision) without the first being comparably small.

A discussion of better ways to apply stopping criteria for iterative methods is a more advanced topic than we want cover in this course module. Just keep this in mind as you continue your exploration of numerical methods in the future!

---

The cell below loads the style of the notebook¶