Lecture 15
I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the this notebook button on the side of this page.
These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.
We are considering boundary value problems of the following form: Suppose we have a differential operator \mathcal L, a “nice” domain \Omega \subset \mathbb R^d, and data f. Seek u : \Omega \to \mathbb R such that
\begin{align} \mathcal L[u] &= f &&\text{in } \Omega \nonumber\\ u &= 0 &&\text{on } \partial \Omega . \tag{BVP} \end{align}
This is a problem with homogeneous Dirichlet boundary conditions (the values of the solution are prescribed on the boundary (Dirichlet) and are zero (homogeneous)).
We also consider the specific problem (in 1d):
Example 15.1 (Linear Advection-Diffusion-Reaction Equations) Consider
\begin{align} \mathcal L[u] := -\Delta u + p \cdot \nabla u + q u \end{align}
In the context of semiconductor physics, u(x) represents the concentration of charge carriers at a position x \in \Omega. Each of these terms relates to the physics in question:
The same equation also models e.g pollutant transport in a river, drug delivery in a blood vessel, asset price in finance, heat transfer in fluids, ……
By approximating the continuous problem with finite differences, we obtain a discrete problem: find \bm U \in \mathbb R^{n+1} (with index starting at 0 or equivalently U: \Omega_h \to \mathbb R where \Omega_h := \{x_j\}_{j=0}^n \to \mathbb R with U_j := U(x_j)) such that
\begin{align} \begin{split} &\mathcal L_h [\bm U] = \bm f \\ &U_0 = U_n = 0 \end{split}\tag{$\text{BVP}_h$} \end{align}
where \mathcal L_h depends on the choice of finite difference scheme.
Last time, we approximated the ADR equation with the following finite difference scheme on the mesh \{x_j\}_{j=0}^n with x_j = \frac{j}{n} and mesh size h= \frac{1}{n}:
\begin{align} \mathcal L_h &= -L + p D^0 + qI \nonumber\\ % &= -\frac{1}{h^2} \begin{pmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ & 1 & -2 & \ddots \\ & & \ddots & \ddots & 1\\ & & & 1 & -2 \end{pmatrix} % + \frac{p}{2h} \begin{pmatrix} 0 & 1 & & & \\ -1 & 0 & 1 & & \\ & -1 & 0 & \ddots \\ & & \ddots & \ddots & 1\\ & & & -1 & 0 \end{pmatrix} % + q I, \nonumber \end{align}
a (n-1)\times(n-1) matrix (acting on the interior points of the grid).
We saw that this method is:
Thus, in these particular cases, the finite difference approximation converges with order p=2.
We now consider the following periodic problem:
\begin{align} \mathcal Lu &:= -u'' + p u' + q u = f \qquad \text{on } (0,1) \nonumber\\ &u(0) = u(1) \\ &u'(0) = u'(1) \nonumber \end{align}
Exercise 15.1 Write down a consistent finite difference approximation to \mathcal L on the grid \{x_j\}_{j=0}^n with x_j = jh and h = \frac{1}{n} (that also incorporates the boundary conditions).
Hint: the matrices are now n\times n.
We thus consider the problem
\begin{align} \mathcal L_h &= -L_{\rm per} + p D^0_{\rm per} + qI \nonumber\\ % &= -\frac{1}{h^2} \begin{pmatrix} -2 & 1 & & & 1 \\ 1 & -2 & 1 & & \\ & 1 & -2 & \ddots \\ & & \ddots & \ddots & 1\\ 1 & & & 1 & -2 \end{pmatrix} % \nonumber\\ &\qquad + \frac{p}{2h} \begin{pmatrix} 0 & 1 & & & -1\\ -1 & 0 & 1 & & \\ & -1 & 0 & \ddots \\ & & \ddots & \ddots & 1\\ 1 & & & -1 & 0 \end{pmatrix} % + q I. \nonumber \end{align}
We would like to prove stability (and thus convergence) of this numerical method. We do this in the 2-norm: we want to show that
\begin{align} \| x \|_2 \leq C \| Ax\|_2 \end{align}
where A := -L_{\rm per} + p D^0_{\rm per} + q I.
We have seen that, if (\lambda_k, v_k) are normalised eigenpairs of A and
\begin{align} A = \sum_k \lambda_k v_k v_k^\star \end{align}
with \lambda_k \not=0, then A is invertible with eigenpairs (\lambda_k^{-1}, v_k). Therefore, if we able to compute the eigenvalues of A := -L_{\rm per} + p D^0_{\rm per} + q I and they are nonzero we can invert A and write
\begin{align} \| x\|_2 &= \|A^{-1} Ax \|_2 \nonumber\\ % &\leq \|A^{-1}\|_2 \| Ax \|_2 \nonumber\\ % &= \rho(A^{-1}) \| Ax \|_2 \nonumber\\ % &= \max_{\lambda \in \sigma(A)} \frac{1}{|\lambda|} \| Ax \|_2. \end{align}
We will use the discrete Fourier transform to diagonalise \mathcal L_h.
Suppose we have some U_j for j = 0,1,...,n-1 (e.g. representing the values of some function on a grid \{x_j\}_{j=0}^{n} (and U_n = U_0)). The idea behind the discrete Fourier transform (DFT) is to write U_j with respect to the Fourier basis: for k = 0,\dots,n-1,
\begin{align} (\bm e_k)_j &:= \frac{1}{\sqrt{n}} e^{\frac{2\pi ijk}{n}}. \end{align}
Exercise 15.2 Show that \bm e_k satisfy the periodic boundary conditions.
We define the inner product (u,v) := \sum_{j=0}^{n-1} \overline{u_j} v_j on \mathbb C^{n} (notice this extends the dot product to complex valued vectors).
Some people put the complex conjugate in the other component \sum_{j=0}^n u_j \overline{v_j}.
We say that u, v are orthogonal if (u,v) = 0. A set \{u^{(j)}\}_j is orthonormal if (u^{(j)}, u^{(j')}) = \delta_{jj'} (i.e. =1 if j=j' and =0 if j\not=j').
Exercise 15.3 Show that \{ \bm e_k \}_{k=0}^{n-1} forms an orthonomal basis of \mathbb C^n.
We can write U with respect to this basis: there exists constants \hat{U}_k such that
\begin{align} U_j = \sum_{k=0}^{n-1} \hat{U}_k (\bm e_k)_j. \end{align}
The constants \hat{U}_k are the Fourier coefficients of U and we can write U = \sum_{k=0}^{n-1} \hat{U}_k \bm e_k.
Exercise 15.4 Use the orthonormality of \bm e_k to write down a formula for \hat{U}_k.
Extra. Show Parseval’s identity: \|U\|_2 = \|\hat{U}\|_2.
Exercise 15.5 Write down (i) D^0_{\rm per} \bm e_k and (ii) -L_{\rm per} \bm e_k. Hence write down the eigenvalues \lambda_k of \mathcal L_h.
Recall, we are approximating the solution to -u'' + p u' + q u = f with PBCs u(0) = u(1) and u'(0) = u'(1) with the finite difference \mathcal L_h = -L_{\rm per} + p D^0_{\rm per} + qI.
Exercise 15.6 Take p = 0 and q>0. Show that the solution U_h of \mathcal L_h U_h = F converges to the solution of the continuous problem as h \to 0.
Exercise 15.7 Let F_j = \sum_{k=0}^{n-1} \hat{F}_k (\bm e_k)_j and U_j = \sum_{k=0}^{n-1} \hat{U}_k (\bm e_k)_j. Rewrite the problem \mathcal L_h U = F in terms of the Fourier coefficients (\hat{U}_k, \hat{F}_k). Explain how you can use DFT to solve the linear system of equations -\mathcal L_h U = F without inverting the matrix \mathcal L_h directly.
Here, we solve (-L_{\rm per} + q I) U = F doing this:
using Plots
q = 1
f(x) = x
n = 2^10
h = 1/n
# this is an inefficient way to do this (see A5)
DFT = [ exp( -2π*1im*j*k/n )/sqrt(n) for k∈0:n-1, j∈0:n-1 ]
Uhat = DFT*f.(0:h:1-h) ./ ( (2/h^2)*[1-cos(2π*k/n) for k∈0:n-1] .+ q )
U = DFT' * Uhat
U = [U..., U[1]]
plot( 0:h:1, real.(U), m=1 )
Previous reading: