Numerical PDEs - Casey Chu

Numerical PDEs
I took these notes in spring 2017 for CME 306, taught by Professor Lexing Ying.
Notation
If $f$ is $|\alpha|$ -differentiable, then $D^\alpha f \in C^{|\alpha|-1}$ is the usual derivative.
If $f$ is weakly differentiable, then $D^\alpha f : \R^n \to \R$ satisfies $\langle D^\alpha f, \phi \rangle = (-1)^{|\alpha|} \langle f, D^\alpha\phi \rangle \qquad\text{for all } \phi \in C^\infty_0.$
If $f$ is differentiable in the distributional sense, then $D^\alpha f : C^\infty_0 \to \R$ satisfies $(D^\alpha f) \{ \phi \} = (-1)^{|\alpha|} \langle f, D^\alpha\phi \rangle \qquad\text{for all } \phi \in C^\infty_0.$
$C^k(\Omega)$ is the space of $k$ -differentiable functions on $\Omega$ , equipped with the norm
$||f||_{C^k(\Omega)} = \max_{|\alpha| \le k} ||D^\alpha f||_{C(\Omega)}$ and the seminorm $|f|_{C^k(\Omega)} = \max_{|\alpha| = k} ||D^\alpha f||_{C(\Omega)}$
$L^2(\Omega)$ is the space of square-integrable functions on $\Omega$ , equipped with the inner product $\langle f, g \rangle = \int_\Omega fg$ and induced norm $||f|| = \int_\Omega f^2.$
$H^k(\Omega)$ is the Sobolev space $H^k(\Omega) = \{ f \in L^2(\Omega) \,|\, D^\alpha f \in L^2(\Omega) \text{ for all } |\alpha| \le k \}$ equipped with the inner product $\langle f, g \rangle_k = \sum_{|\alpha| \le k} \langle D^\alpha f, D^\alpha g \rangle_{L^2(\Omega)}.$
- Poincare inequality: $||f|| \le C ||D f||$ .
- Note $|f|_1^2 \le ||f||_1^2 \le (C^2 + 1)|f|_1^2$ for $H_0^1$
Finite-difference methods
In finite-difference methods, derivatives are approximated with forward/backwards/central difference operators: $\partial U_j = \frac{U_{j+1} - U_j}{h}, \qquad \bar\partial U_j = \frac{U_{j} - U_{j-1}}{h}, \qquad \hat\partial U_j = \frac{U_{j+1} - U_{j-1}}{2h}.$
They may be combined to yield, for example, approximations to the second derivative:
$\partial\bar\partial U_j = \frac{U_{j+1} - 2U_j + U_{j-1}}{h^2}, \qquad \Delta_h U_j = \sum_\ell \partial_\ell\bar\partial_\ell U_j.$
Error estimates (truncation error) are obtained using Taylor’s theorem: $\begin{aligned} |\partial u_j - u'(x_j)| &\le Ch |u|_{C^2(x_{j-1}, x_j)} \\ |\hat\partial u_j - u'(x_j)| &\le Ch^2 |u|_{C^2(x_{j-1}, x_{j+1})} \\ |\partial\bar\partial u_j - u''(x_j)| &\le Ch^2 |u|_{C^4(x_{j-1}, x_{j+1})} \\ |\Delta_h u_j - \Delta u(x_j)| &\le Ch^2 |u|_{C^4}. \end{aligned}$
Poisson’s equation
To solve Poisson’s equation with finite differences, replace $\begin{aligned} -\Delta u &= f \quad\text{ in }\Omega &\quad\Longrightarrow\quad -\Delta_h U_j &= f_j \quad \text{ for } x_j \in \Omega \\ u &= 0 \quad\text{ on }\partial\Omega & U_j &= 0 \quad \text{ for } x_j \in \partial\Omega, \end{aligned}$ assuming that all neighbors of every point $x_j$ are contained in $\bar\Omega$ . This yields a system of linear equations in the $U_j$ , which we may solve.
Consistency is given by the error estimate $|U-u|_\Omega \le Ch^2 |u|_{C^4}$ by the following argument.
- Discrete maximum principle: If $-\Delta_h U \le 0$ for all $x_j \in \Omega$ , then $U$ attains its maximum on $\partial \Omega$ . (Proof: assume it attains an interior maximum; then its neighbors are also maxima.)
- Bound on growth: $|U|_{\bar\Omega} \le |U|_\Gamma + C |\Delta_h U|_\Omega$ . (Proof: Recalling for a paraboloid $W$ , $\Delta_h W = \Delta W = \text{constant}$ , construct a paraboloid for which $\Delta_h W_j = -|\Delta_h U|_\Omega$ and $0 \le W_j \le C$ . Then apply the discrete maximum principle to the mesh function $V^\pm_j = \pm U_j - W_j$ .)
- This implies uniqueness (and hence existence), applying the bound to $U_1 - U_2$ .
- This implies the error estimate, applying the bound to $U- u$ and using the truncation error on $\Delta_h$ : $|U- u|_\Omega \le C |\Delta_h U - \Delta_h u| = |f - \Delta_h u| = |\Delta u - \Delta_h u| \le Ch^2 |u|_{C^4}.$
Initial-value problems
Suppose we are solving an initial value problem of the form $u_t = p(\tfrac{\partial}{\partial x})u,$ where $p$ is some polynomial. With finite differences, we construct an evolution operator $E_k$ that describes the approximate solution at the next time step $t+k$ : $U_j^{n+1} = E_k U_j^n$ and compute iterates upon an initial condition $U^0_j$ . A necessary condition for convergence to the correct solution is both consistency and stability.
For analysis, view $E_k$ as a convolution operator acting on a function $v(x)$ : $E_kv = \sum_p a_p v(x - hp)\, dx,$ so that the Fourier transform of both sides is $(E_k v)\hat{} = \hat{E}(h\xi) \hat{v}(\xi),$ with $\hat{E}(\xi) = \sum_p a_p e^{-ip\xi}.$
For an operator $E_k$ to be consistent, we require $E_k v \sim E(k) v$ as $k \to 0$ and $h \to 0$ , where $E(t)$ is the true time evolution operator (that evolves by time $t$ ). Generally, a scheme created via finite differences will be consistent. However, if we need to constrct an operator by hand as a linear combination of lattice points, there is a nice characterization in the Fourier domain. By Taylor expansion in $t$ and using the PDE, the true evolution operator is $E(k)v = \exp(k\tfrac{\partial}{\partial t})v = \exp(kp(\tfrac{\partial}{\partial x}))v.$ Under the Fourier transform, $\frac{\partial}{\partial x} \mapsto i\xi$ , so this equation becomes $(E(k)v)\hat{} = \exp(kp(i\xi))\hat{v}.$ Therefore, $E_k$ is accurate to order $r$ , i.e. $E_k v = E(k)v + O(h^{r+1} v^{(r+1)}),$ if and only if $\hat{E}_k(h\xi) = \exp(kp(i\xi)) + O((h\xi)^{r+1}),$ or, rescaling: $\hat{E}_k(\xi) = \exp(kp(\tfrac{i\xi}{h})) + O(\xi^{r+1}).$
For an operator to be stable w.r.t. some norm, we require $||E_k^n u|| \le ||u||$ for all $n$ , where now we think of $u$ as the error. If an operator is unstable, the error’s magnitude may increase exponentially. (Note that the PDE itself must be stable, which is true for the heat equation and transport equation.) Von Neumann stability analysis uses the $L_2$ -norm, which especially nice because the Plancherel theorem says that $||v||^2_2 = (2\pi)^{-1}||\hat{v}||^2_2$ . Thus, the stability condition can be written $||\hat{E}_k^n \hat{u}||_{2} \le ||\hat{u}||_{2}$ , which holds if and only if $||\hat{E}_k||_\infty \le 1$ . A necessary condition for stability is the CFL condition (Courant-Friedrichs-Lewy), which says that a necessary condition for stability is that the domain of $E_k$ contains the domain of $E(k)$ .
If $E_k$ is both consistent and stable, we may show convergence via the following error estimate. Similar to above, $\begin{aligned} ||U^n - u_n||_2^2 &= ||\big(\hat{E}(h\xi)^n - \exp(nkp(i\xi))\big) \hat{v}||_2^2 \\ &= ||\big(\hat{E}(h\xi) - \exp(kp(i\xi))\big) \sum_{j=0}^{n-1} \hat{E}(h\xi)^{n-1-j}\exp(jkp(i\xi)) \hat{v}||_2^2 \\ &\le Cn||\big(\hat{E}(h\xi) - \exp(kp(i\xi))\big)\hat{v}||_2^2 \\ &\le Cn||(h\xi)^{r+1}\hat{v}||_2^2 \\ &\le C||nh^{r+1}v^{(r+1)}||_2^2. \end{aligned}$
Heat equation
For the heat equation, the condition for $r$ th-order accuracy becomes $\hat{E}_k(\xi) = e^{-\lambda \xi^2} + O(\xi^{r+2}),$ where $\lambda = \frac{k}{h^2}$ . (We may increase the order by 1 in error term because the $p$ only has even powers.)
The $\theta$ -method uses the evolution operator given by $\bar\partial_t U_j^{n+1} = \theta\partial_x \bar\partial_x U_j^{n+1} + (1-\theta)\partial_x \bar\partial_x U_j^n.$ The explicit forward Euler scheme uses $\theta = 0$ , is a second-order scheme, and is $L_2$ -stable for $\lambda \le \frac{1}{2}$ .
The implicit forward Euler scheme uses $\theta = 1$ and is so-named because we need to solve a linear system to determine $U^{n+1}$ . It is stable for all $\lambda$ . The mesh sizes $h$ and $k$ need not be related anymore, and the max-norm error estimate can be written $||U^n - u^n||_\infty \le C t(h^2 + k) \max_{\tau \le t} |u(\cdot, t)|_{C^4}.$
The Crank-Nicolson scheme uses $\theta = \frac{1}{2}$ . It is stable for all $\lambda$ . The max-norm error estimate can be written $||U^n - u^n||_\infty \le C t(h^2 + k^2) \max_{\tau \le t} |u(\cdot, t)|_{C^6}.$
Transport equation
For the scalar transport equation $u_t = au_x$ , the condition for $r$ th-order accuracy becomes $\hat{E}_k(\xi) = e^{i\lambda a\xi} + O(\xi^{r+1}),$ where $\lambda = \frac{k}{h}$ . The von Neumann stability condition for the following three schemes is $|a\lambda| \le 1$ .
The forward Euler scheme uses the evolution operator given by $\partial_t U_j^n = \begin{cases} a \partial_x U_j^n & \text{ if } a > 0, \\ a \bar\partial_x U_j^n & \text{ if } a < 0, \end{cases}$ where the dependence on the sign of $a$ is called upwinding (necessary to satisfy the CFL condition).
The Friedrichs scheme is a first-order method based on the central difference $\partial_t U^n = a\hat\partial_x U^n + \frac{1}{2}\frac{h}{\lambda}\partial_x \bar\partial_x U^n,$ where the artificial diffusion stabilizes it.
The Lax-Wendroff scheme is a method constructed explicitly to be second-order accurate: $U_j^{n+1} = aU^n_{j-1} + bU^n_{j} + cU^n_{j+1},$ where $a = \tfrac{1}{2}(a^2\lambda^2 + a\lambda), \qquad b = 1 - a^2\lambda^2, \qquad c = \tfrac{1}{2}(a^2\lambda^2 - a\lambda).$
The Wendroff box scheme is a second-order method suitable for mixed initial-boundary value problems: $\partial_t U_{j+1/2}^n + a_{j+1/2}^{n+1/2} \partial_x U_j^{n+1/2} + b_{j+1/2}^{n+1/2} U_{j+1/2}^{n+1/2} = f_{j+1/2}^{n+1/2},$ where averages of surrounding grid points (up to four) are used if the grid point doesn’t exist.
Finite-element method
The Galerkin method is the following. We relax the strong formulation of the PDE $Au = f$ into the weak formulation $( Au, v ) = ( f, v ) \qquad \forall v \in H$ for some bilinear form $( \cdot, \cdot )$ and space of test functions $H$ . Define a bilinear form $\phi$ as $\phi(u, v) = ( Au, v ),$ often symmetric.
If we relax $H$ to some subspace $W$ with basis $v_1, \ldots, v_n$ , we may express $u = \sum_i u_i v_i,$ so that this becomes $\sum_i u_i \phi( v_i, v_j ) = ( f, v_j ) \qquad j = 1, \ldots, n.$ In matrix form, this is the linear system $Au = b$ , with the matrices $A_{ij} = \phi(v_i, v_j), \qquad b_j = ( f, v_j ).$ The matrix $A$ is called the stiffness matrix and $b$ is called the load vector.
Let’s consider Poisson’s equation $-\nabla \cdot (a\nabla u) = f$ with Dirichlet boundary conditions, with $a \ge a_0 > 0$ (a coercivity condition). Here, the space of test functions is the Sobolev space $H_0^1$ with the bilinear forms $(v,w) = \int_\Omega vw, \qquad \phi(v,w) = -a\int_\Omega \nabla v \cdot \nabla w.$ For the finite-element method, the subspace is the space of piecewise linear functions on meshpoints $\{ P_j \}$ , with basis of the “pyramid functions” $\{ \Phi_i \}$ interpolated from $\Phi_i(P_j) = \delta_{ij}.$
Error analysis is most natural in the energy norm $||v||_a = \sqrt{\phi(v,v)}$ . Since $\phi(U, \chi) = (f, \chi)$ and $\phi(u, \chi) = (f, \chi)$ for all $\chi \in W$ by construction, we have $\phi(U-u, \chi) = 0,$ so choosing $\chi = \chi_* -U$ , $\begin{aligned} ||U-u||_a^2 &= \phi(U-u, U-u) + \phi(U-u, \chi_* -U) \\ &= \phi(U-u, \chi_* - u) \\ &\le ||U-u||_a ||\chi_* - u||_a, \end{aligned}$ so $||U-u||_a = \min_{\chi \in W} ||\chi_* - u||_a.$ We assumed that $\phi$ is symmetric to apply Cauchy-Schwarz.
We may be interested in bounding with the $L_2$ norm. For this, first convert norms from above using the relation $c||v||_{L_2} \le ||v||_a \le C||v||_{L_2},$ where the first follows is a coercivity condition (which follows from $a \ge a_0$ and Poincaré’s inequality) and the second is a boundedness condition (???). Then set $\chi = I_h u$ to be the linear interpolation of $u$ and use the interpolation bound $||I_h v - v||_{L_2} \le Ch^2|v|_{H^2}$ to arrive at $||U - u||_{L_2} \le Ch^2|u|_{H^2} \le Ch^2 |\Delta u|_{L_2},$ where the last inequality is true when $\Omega$ is convex.
We may be also interested in bounding with the $H^1$ seminorm $|v|_{H^1} = || \nabla v||_{L_2}$ . For this, first convert norms from above using the relation $c|v|_{H^1} \le ||v||_a \le C|v|_{H^1},$ which the first follows from our $a \ge a_0$ condition and the second is Poincaré’s inequality. Then set $\chi = I_h u$ to be the linear interpolation of $u$ and use the interpolation bound $|I_h v - v|_{H^1} \le Ch|v|_{H^2}$ to arrive at $|U - u|_{H^1} \le Ch|u|_{H^2} \le Ch ||\Delta u||_{L_2},$ where the last inequality is true when $\Omega$ is convex.
Finite-volume method
Conservative schemes discretize the spatial domain into cells $C_j = [x_j, x_{j+1}]$ and locally solve the integral equation $\frac{d}{dt}\int_{x_j}^{x_{j+1}} u(x,t) \,dx = f(u(x_j,t)) - f(u(x_{j+1},t)).$ Approximate $U_j \approx \frac{1}{h}\int_{x_j}^{x_{j+1}}u(x,t) \,dx$ to arrive at the conservative scheme $h\partial_t U_j^n = \hat{f}(U_{j-1}^n, U_{j}^n) - \hat{f}(U_{j}^n, U_{j+1}^n)$ with numerical flux $\hat{f}(U_L, U_R)$ .
A conservative scheme is consistent if $\hat{f}$ locally satisfies a Lipschitz condition $|\hat{f}(U_L, U_R) - f(u)| \le C \max(|U_L- u|, |U_R- u|).$ If a consistent conservative scheme converges, then it converges to a weak solution (Lax and Wendroff).
If additionally the scheme is monotone, i.e. the scheme is written $U_j^{n+1} = G(U_{j-1}^n, U_{j}^n, U_{j+1}^n)$ and $G$ is monotonically increasing in its three arguments, then the scheme converges to an entropy solution (Crandell and Majda). This condition may be expanded to read $\partial_1\hat{f} \ge 0, \qquad \partial_2\hat{f} \le 0, \qquad \frac{h}{k} \ge \partial_1 \hat{f}(U_j, U_{j+1}) - \partial_2\hat{f}(U_{j-1}, U_j).$ Monotone schemes are at best first-order.
Define $\alpha = \max_{|u| \le A} |f'(u)|$ and $A \ge |u(x,0)|$ , where we can think of $\alpha$ as the maximum velocity of flow. The following schemes are monotone with $\frac{h}{k} \ge \alpha,$ a CFL condition.
The Godunov scheme uses the numerical flux $\hat{f}(U_L, U_R) = \begin{cases} \min_{U_L \le u \le U_R} f(u) & U_L \le U_R \\ \max_{U_R \le u \le U_L} f(u) & U_R \le U_L, \end{cases}$ derived from solving the Riemann problem at each mesh point (requires casework).
The Lax-Friedrichs scheme uses the numerical flux $\hat{f}(U_L, U_R) = \tfrac{1}{2}[f(U_L) + f(U_R) - \alpha(U_R - U_L) ].$