Demo - 1D Poisson’s equation

Authors

Mikael Mortensen (mikaem at math.uio.no)

Date

Nov 14, 2019

Summary. This is a demonstration of how the Python module shenfun can be used to solve Poisson’s equation with Dirichlet boundary conditions in one dimension. Spectral convergence, as shown in Figure Convergence of 1D Poisson solvers for both Legendre and Chebyshev modified basis function, is demonstrated. The demo is implemented in a single Python file dirichlet_poisson1D.py, and the numerical method is is described in more detail by J. Shen [She94] and [She95].

https://rawgit.com/spectralDNS/spectralutilities/master/figures/poisson1D_errornorm.png

Convergence of 1D Poisson solvers for both Legendre and Chebyshev modified basis function

Model problem

Poisson’s equation

Poisson’s equation is given as

(1)\[ \nabla^2 u(x) = f(x) \quad \text{for }\, x \in [-1, 1],\]
(2)\[ u(-1)=a, u(1)=b, \notag\]

where \(u(x)\) is the solution, \(f(x)\) is a function and \(a, b\) are two possibly non-zero constants.

To solve Eq. (1) with the Galerkin method we need smooth continuously differentiable basis functions, \(v_k\), that satisfy the given boundary conditions. And then we look for solutions like

(3)\[ u(x) = \sum_{k=0}^{N-1} \hat{u}_k v_k(x),\]

where \(N\) is the size of the discretized problem, \(\hat{\mathbf{u}} = \{\hat{u}_k\}_{k=0}^{N-1}\) are the unknown expansion coefficients, and the basis is \(\text{span}\{v_k\}_{k=0}^{N-1}\).

The basis functions can, for example, be constructed from Chebyshev, \(T_k(x)\), or Legendre, \(L_k(x)\), polynomials and we use the common notation \(\phi_k(x)\) to represent either one of them. It turns out that it is easiest to use basis functions with homogeneous Dirichlet boundary conditions

(4)\[ v_k(x) = \phi_k(x) - \phi_{k+2}(x),\]

for \(k=0, 1, \ldots N-3\). This gives the basis \(V^N_0 = \text{span}\{v_k(x)\}_{k=0}^{N-3}\). We can then add two more linear basis functions (that belong to the kernel of Poisson’s equation)

(5)\[ v_{N-2} = \frac{1}{2}(\phi_0 + \phi_1),\]
(6)\[ v_{N-1} = \frac{1}{2}(\phi_0 - \phi_1).\]

which gives the inhomogeneous basis \(V^N = \text{span}\{v_k\}_{k=0}^{N-1}\). With the two linear basis functions it is easy to see that the last two degrees of freedom, \(\hat{u}_{N-2}\) and \(\hat{u}_{N-1}\), now are given as

(7)\[ u(-1) = \sum_{k=0}^{N-1} \hat{u}_k v_k(-1) = \hat{u}_{N-1} = a,\]
(8)\[ u(+1) = \sum_{k=0}^{N-1} \hat{u}_k v_k(+1) = \hat{u}_{N-2} = b,\]

and, as such, we only have to solve for \(\{\hat{u}_k\}_{k=0}^{N-3}\), just like for a problem with homogeneous boundary conditions (for homogeneous boundary condition we simply have \(\hat{u}_{N-2} = \hat{u}_{N-1} = 0\)). We now formulate a variational problem using the Galerkin method: Find \(u \in V^N\) such that

(9)\[ \int_{-1}^1 \nabla^2 u \, v \, w\, dx = \int_{-1}^1 f \, v\, w\, dx \quad \forall v \, \in \, V^N_0.\]

Note that since we only have \(N-3\) unknowns we are only using the homogeneous test functions from \(V^N_0\).

The weighted integrals, weighted by \(w(x)\), are called inner products, and a common notation is

(10)\[ \int_{-1}^1 u \, v \, w\, dx = \left( u, v\right)_w.\]

The integral can either be computed exactly, or with quadrature. The advantage of the latter is that it is generally faster, and that non-linear terms may be computed just as quickly as linear. For a linear problem, it does not make much of a difference, if any at all. Approximating the integral with quadrature, we obtain

\[\begin{split}\begin{align*} \int_{-1}^1 u \, v \, w\, dx &\approx \left( u, v \right)_w^N, \\ &\approx \sum_{j=0}^{N-1} u(x_j) v(x_j) w(x_j), \end{align*}\end{split}\]

where \(\{w(x_j)\}_{j=0}^{N-1}\) are quadrature weights. The quadrature points \(\{x_j\}_{j=0}^{N-1}\) are specific to the chosen basis, and even within basis there are two different choices based on which quadrature rule is selected, either Gauss or Gauss-Lobatto.

Inserting for test and trialfunctions, we get the following bilinear form and matrix \(A\in\mathbb{R}^{N-3\times N-3}\) for the Laplacian (using the summation convention in step 2)

\[\begin{split}\begin{align*} \left( \nabla^2u, v \right)_w^N &= \left( \nabla^2\sum_{k=0}^{N-3}\hat{u}_k v_{k}, v_j \right)_w^N, \quad j=0,1,\ldots, N-3\\ &= \left(\nabla^2 v_{k}, v_j \right)_w^N \hat{u}_k, \\ &= a_{jk} \hat{u}_k. \end{align*}\end{split}\]

Note that the sum in \(a_{jk} \hat{u}_{k}\) runs over \(k=0, 1, \ldots, N-3\) since the second derivatives of \(v_{N-1}\) and \(v_{N}\) are zero. The right hand side linear form and vector is computed as \(\tilde{f}_j = (f, v_j)_w^N\), for \(j=0,1,\ldots, N-3\), where a tilde is used because this is not a complete transform of the function \(f\), but only an inner product.

The linear system of equations to solve for the expansion coefficients of \(u(x)\) is given as

(11)\[ A \hat{\mathbf{u}} = \tilde{\mathbf{f}}.\]

Now, when the expansion coefficients \(\hat{\mathbf{u}}\) are found by solving this linear system, they may be transformed to real space \(u(x)\) using (3), and here the contributions from \(\hat{u}_{N-2}\) and \(\hat{u}_{N-1}\) must be accounted for. Note that the matrix \(A\) (different for Legendre or Chebyshev) has a very special structure that allows for a solution to be found very efficiently in order of \(\mathcal{O}(N)\) operations, see [She94] and [She95]. These solvers are implemented in shenfun for both bases.

Method of manufactured solutions

In this demo we will use the method of manufactured solutions to demonstrate spectral accuracy of the shenfun Dirichlet bases. To this end we choose an analytical function that satisfies the given boundary conditions:

(12)\[ u_e(x) = \sin(k\pi x)(1-x^2) + a(1+x)/2 + b(1-x)/2,\]

where \(k\) is an integer and \(a\) and \(b\) are constants. Now, feeding \(u_e\) through the Laplace operator, we see that the last two linear terms disappear, whereas the first term results in

(13)\[ \nabla^2 u_e(x) = \frac{d^2 u_e}{dx^2},\]
(14)\[ = -4k \pi x \cos(k\pi x) - 2\sin(k\pi x) - k^2 \pi^2 (1 - x^2) \sin(k \pi x).\]

Now, setting \(f_e(x) = \nabla^2 u_e(x)\) and solving for \(\nabla^2 u(x) = f_e(x)\), we can compare the numerical solution \(u(x)\) with the analytical solution \(u_e(x)\) and compute error norms.

Implementation

Preamble

We will solve Poisson’s equation using the shenfun Python module. The first thing needed is then to import some of this module’s functionality plus some other helper modules, like Numpy and Sympy:

from shenfun import inner, div, grad, TestFunction, TrialFunction, Function, \
    project, Dx, Array, Basis
import numpy as np
from sympy import symbols, cos, sin, exp, lambdify

We use Sympy for the manufactured solution and Numpy for testing.

Manufactured solution

The exact solution \(u_e(x)\) and the right hand side \(f_e(x)\) are created using Sympy as follows

a = -1
b = 1
k = 4
x = symbols("x")
ue = sin(k*np.pi*x)*(1-x**2) + a*(1 + x)/2. + b*(1 - x)/2.
fe = ue.diff(x, 2)

# Lambdify for faster evaluation
ul = lambdify(x, ue, 'numpy')
fl = lambdify(x, fe, 'numpy')

These solutions are now valid for a continuous domain. The next step is thus to discretize, using a discrete mesh \(\{x_j\}_{j=0}^{N-1}\) and a finite number of basis functions.

Note that it is not mandatory to use Sympy for the manufactured solution. Since the solution is known (14), we could just as well simply use Numpy to compute \(f_e\) at \(\{x_j\}_{j=0}^{N-1}\). However, with Sympy it is much easier to experiment and quickly change the solution.

Discretization

We create a basis with a given number of basis functions, and extract the computational mesh from the basis itself

N = 32
SD = Basis(N, 'Chebyshev', bc=(a, b))
#SD = Basis(N, 'Legendre', bc=(a, b))
X = SD.mesh(N)

Note that we can either choose a Legendre or a Chebyshev basis.

Variational formulation

The variational problem (12) can be assembled using shenfun’s TrialFunction, TestFunction and inner() functions.

u = TrialFunction(SD)
v = TestFunction(SD)
# Assemble left hand side matrix
A = inner(v, div(grad(u)))
# Assemble right hand side
fj = Array(SD, buffer=fl(X))
f_hat = Function(SD)
f_hat = inner(v, fj, output_array=f_hat)

Note that fl(X) returns a Numpy array of the correct shape and type of the left hand side of (3), evaluated on all quadrature points X. We wrap this Numpy array in an Array class (fj = Array(SD, buffer=fl(X))), because an Array is required as input to the inner() function.

Solve linear equations

Finally, solve linear equation system and transform solution from spectral \(\{\hat{u}_k\}_{k=0}^{N-1}\) vector to the real space \(\{u(x_j)\}_{j=0}^{N-1}\) and then check how the solution corresponds with the exact solution \(u_e\).

u_hat = A.solve(f_hat)
uj = SD.backward(u_hat)
ue = ul(X)
print("Error=%2.16e" %(np.linalg.norm(uj-ue)))
assert np.allclose(uj, ue)

Convergence test

A complete solver is given in Sec. Complete solver. This solver is created such that it takes in two commandline arguments and prints out the \(l_2\)-errornorm of the solution in the end. We can use this to write a short script that performs a convergence test. The solver is run like

>>> python dirichlet_poisson1D.py 32 legendre
Error=6.5955040031498912e-10

for a discretization of size \(N=32\) and for the Legendre basis. Alternatively, change legendre to chebyshev for the Chebyshev basis.

We set up the solver to run for a list of \(N=[12, 16, \ldots, 48]\), and collect the errornorms in arrays to be plotted. Such a script can be easily created with the subprocess module

import subprocess

N = range(12, 50, 4)
error = {}
for basis in ('legendre', 'chebyshev'):
    error[basis] = []
    for i in range(len(N)):
        output = subprocess.check_output("python dirichlet_poisson1D.py {} {}".format(N[i], basis), shell=True)
        exec(output) # Error is printed as "Error=%2.16e"%(np.linalg.norm(uj-ua))
        error[basis].append(Error)

The error can be plotted using matplotlib, and the generated figure is shown in the summary’s Fig. Convergence of 1D Poisson solvers for both Legendre and Chebyshev modified basis function. The spectral convergence is evident and we can see that after \(N=40\) roundoff errors dominate as the errornorm trails off around \(10^{-14}\).

import matplotlib.pyplot as plt
plt.figure(figsize=(6, 4))
for basis, col in zip(('legendre', 'chebyshev'), ('r', 'b')):
    plt.semilogy(N, error[basis], col, linewidth=2)
plt.title('Convergence of Poisson solvers 1D')
plt.xlabel('N')
plt.ylabel('Error norm')
plt.savefig('poisson1D_errornorm.png')
plt.legend(('Legendre', 'Chebyshev'))
plt.show()

Complete solver

A complete solver, that can use either Legendre or Chebyshev bases, chosen as a command-line argument, can be found here.

PSS18

Ambrish Pandey, Janet D. Scheel, and Jörg Schumacher. Turbulent superstructures in rayleigh-bénard convection. Nature Communications, 9(1):2118, 2018. doi:10.1038/s41467-018-04478-0.