MATH566 Lesson 27: Numerical PDE - Introduction

Overview

• First-order PDEs

  advection equation

  convection equation

  characteristic solution

• Second-order PDE classification, canonical forms

  hyperbolic, wave equation

  parabolic, heat equation

  elliptic, Poisson equation

• Reformulating second-order PDEs as first-order PDE system, eigenproblems

• Overview of numerical method development: finite differences, finite volume, finite element, spectral methods

• Finite difference example: leapfrog discretization of wave equation

First-order PDEs

• Many (most) phenomena depend on multiple independent variables

• Natural phenomena are governed by conservation laws (mass, momentum, energy, charge): change in quantity $q(t,x)$ in an infinitesimal volume at time $t$ and position $x$ equals difference of what is going out/in and what was produced

  $=-+\sigma (t,x,q)$

  $f(t,x,q)$ is the flux of quantity $q$, $\sigma (t,x,q)$ is the source of quantity $q$.

• Advection: transport of quantity $q$ in space $x$ and time $t$ by velocity field $u$

  • Constant velocity advection IBVP, $u=$constant, flux $f=uq$

    ${\displaystyle q_t+u{q}_x=0,q(x,t=0)=f\left(x\right),q(x=0,t)=g\left(t\right),q(x=1,t)=h\left(t\right).}$

  • Variable velocity advection, $u(x,t)$, same equations as above

  Examples: transport of a pollutant in a river, drug in the blood stream

• Convection: transport of quantity $q$ in space-time $(x,t)$ by a velocity field that depends on $q$, e.g., Burgers' equation for $q(t,x)$, $q_t\equiv \partial q/\partial t$, $q_x\equiv \partial q/\partial x$

  ${\displaystyle q_t+q{q}_x=0(\mathrm{nonlinear},\mathrm{similar}\mathrm{IBVP}\mathrm{conditions}\mathrm{as}\mathrm{above})}$

Analytical geometry refresher - quadratic classification

• Theory of conic sections highlights quadratic forms in $(x,y)$ coordinates

  ${\displaystyle A{x}^2+2Bxy+C{y}^2+Dx+Ey-G=0\Rightarrow }$ ${\displaystyle \left[\begin{array}{ll}x& y\end{array}\right]\begin{array}{l}\left[\begin{array}{ll}A& B\\ B& C\end{array}\right]\left[\begin{array}{l}x\\ y\end{array}\right]+\end{array}\left[\begin{array}{ll}D& E\end{array}\right]\begin{array}{l}\left[\begin{array}{l}x\\ y\end{array}\right]-G=0\end{array},𝑴=\left[\begin{array}{ll}A& B\\ B& C\end{array}\right],𝒒=\left[\begin{array}{l}x\\ y\end{array}\right]}$

  $M$ symmetric $\Rightarrow$ orthogonal diagonalizable, real eigenvalues, $𝑴=𝑼𝚲{𝑼}^T$.

  ${\displaystyle {𝒒}^T(𝑼𝚲{𝑼}^T)𝒒+{𝒄}^T𝒒=-F,𝒛={𝑼}^T𝒒\Rightarrow {𝒛}^T𝚲𝒛+{𝒎}^T𝒛=-F}$

• Denote $𝚲=\mathrm{diag}(a,b)$, consider $G=0$ (homogeneous), ${𝒎}^T={𝒄}^T𝑼=\left[\begin{array}{ll}2s& 2t\end{array}\right]$. Quadratic form becomes ${𝒛}^T𝚲𝒛+{𝒎}^T𝒛=0$ under change of coordinates

  ${\displaystyle 𝒛=\left[\begin{array}{l}u\\ v\end{array}\right]={𝑼}^T𝒒={𝑼}^T\left[\begin{array}{l}x\\ y\end{array}\right]}$

  • $ab>0\Rightarrow {\xi }^2+{\eta }^2=0$ an ellipse, $\xi =\sqrt{a}u+s/\sqrt{a}$, $\eta =\sqrt{b}v+t/\sqrt{b}$

  • $ab<0\Rightarrow {\xi }^2-{\eta }^2=0$, a hyperbola

  • $ab=0\Rightarrow {\xi }^2=\eta$, a parabola

Second order PDE classification

• Mathematical physics highlights certain ubiquitous PDEs of form

  ${\displaystyle A{u}_{xx}+2B{u}_{xy}+C{u}_{yy}+D{u}_x+E{u}_y+Fu=G}$

• Simplest case: $A,\dots ,G$ constant, linear PDE, classified similar to quadratics

  ${\displaystyle (\left[\begin{array}{ll}{\partial }_x& {\partial }_y\end{array}\right]\begin{array}{l}\left[\begin{array}{ll}A& B\\ B& C\end{array}\right]\left[\begin{array}{l}{\partial }_x\\ {\partial }_y\end{array}\right]+\end{array}\left[\begin{array}{ll}D& E\end{array}\right]\left[\begin{array}{l}{\partial }_x\\ {\partial }_y\end{array}\right]+F)u=G}$

• As in case of quadratics, changes of variables lead to canonical forms

  • Poisson equation $u_{\xi \xi }+u_{\eta \eta }=f$, an elliptical PDE

  • Wave equation $u_{\xi \xi }-u_{\eta \eta }=f$, a hyperbolic PDE

  • Heat equation $u_{\xi \xi }-u_{\eta }=f$, a parabolic PDE

• The above classification is of special relevance to numerical analysis since different numerical methods are applicable for each type of equation

Reformulating second-order PDEs as first-order systems

• Homogeneous wave equation $u_{tt}-u_{xx}=0$, $v\equiv u_t,w=u_x,u_{tx}=u_{xt}\Rightarrow$

  ${\displaystyle \begin{array}{l}{v}_t=w_x\\ w_t=v_x\end{array},𝒒=\left[\begin{array}{l}v\\ w\end{array}\right],{𝒒}_t+𝑨{𝒒}_x=𝟎,𝑨=\left[\begin{array}{ll}0& -1\\ -1& 0\end{array}\right]}$ ${\displaystyle 𝑨=𝑼𝚲{𝑼}^T=\frac{1}{\sqrt{2}}\left[\begin{array}{ll}1& 1\\ -1& 1\end{array}\right]\left[\begin{array}{ll}-1& 0\\ 0& 1\end{array}\right]\frac{1}{\sqrt{2}}\left[\begin{array}{ll}1& -1\\ 1& 1\end{array}\right]}$

  Second-order wave equation leads to a first order system with real eigenvalues

• Homogeneous Poisson equation $u_{xx}+u_{yy}=0$, $v\equiv u_x,w=u_y,u_{xy}=u_{yx}\Rightarrow$

  ${\displaystyle \begin{array}{l}{v}_y=w_x\\ w_y=-v_x\end{array},𝒒=\left[\begin{array}{l}v\\ w\end{array}\right],{𝒒}_y+𝑨{𝒒}_x=𝟎,𝑨=\left[\begin{array}{ll}0& -1\\ 1& 0\end{array}\right]}$ ${\displaystyle 𝑨=𝑼𝚲{𝑼}^T=\frac{1}{\sqrt{2}}\left[\begin{array}{ll}i& 1\\ -i& 1\end{array}\right]\left[\begin{array}{ll}i& 0\\ 0& -i\end{array}\right]\frac{1}{\sqrt{2}}\left[\begin{array}{ll}i& -i\\ 1& 1\end{array}\right]}$

  Second-order elliptic equation leads to a first order system with imaginary eigenvalues

PDE numerical method development

• Basic ideas: discretize both operators (${\partial }_t,{\partial }_x$) or discretize only one operator (typically ${\partial }_x$) and reduce to an ODE system (typically in $t$)

• Approaches:

  • finite difference discretization of differentiation operators, $u_i^n=u(nk,ih)$

    ${\displaystyle u_{tt}-u_{xx}=0\Rightarrow \frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{h^2}-\frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{k^2}=0}$

  • finite difference discretization of ${\partial }_{xx}$ operator only, $u_i\left(t\right)=u(t,ih)$

    ${\displaystyle \frac{d^2{u}_i}{d{t}^2}=\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}}$

  • introduce a piecewise approximation in space for $u(t,x)$

    ${\displaystyle u(t,x)\cong U_i\left(t\right)+[U_{i+1}\left(t\right)-U_i\left(t\right)](x-x_i)/(x_{i+1}-x_i)}$

    a finite element method.

• Different approximations of $u(t,x)$ lead to finite volume, spectral methods.

Numerical PDE example: leap-frog in time, centered in space

• Consider the wave equation $u_{tt}-u_{xx}=0$ with initial, boundary conditions

  ${\displaystyle u(0,x)=\sin x,u_t(0,x)=0,u(t,0)=0,u(t,\pi )=0.}$

  This is known as the plucked string problem, and models a guitar string plucked at midpoint.

• Construct a numerical method by introducing a centered derivative approximation in space, $x_j=jh$, $h=\pi /m$

  ${\displaystyle u_{xx}(t,jh)\cong \frac{u_{i+1}\left(t\right)-2u_i\left(t\right)+u_{i-1}\left(t\right)}{h^2},j=1,\dots ,m-1}$

• Replace above approximation in wave equation at $x=x_j$

  ${\displaystyle \frac{d^2{u}_j\left(t\right)}{d{t}^2}=\frac{u_{i+1}\left(t\right)-2u_i\left(t\right)+u_{i-1}\left(t\right)}{h^2},j=1,\dots ,m-1}$

• Transform second-order ODE into two first-order ODEs

  ${\displaystyle \frac{d{v}_j}{dt}=\frac{u_{i+1}\left(t\right)-2u_i\left(t\right)+u_{i-1}\left(t\right)}{h^2},\frac{d{u}_j}{dt}=v_j}$

Numerical PDE example: leap-frog in time, centered in space

• Obtain a system of ODEs

  ${\displaystyle \frac{d𝒒}{dt}=𝑴𝒒,𝒒=\left[\begin{array}{l}𝒖\\ 𝒗\end{array}\right],𝒖=\left[\begin{array}{l}{u}_1\\ ⋮\\ u_{m-1}\end{array}\right],𝒗=\left[\begin{array}{l}{v}_1\\ ⋮\\ v_{m-1}\end{array}\right],𝑴=\left[\begin{array}{ll}𝟎& 𝑰\\ 𝑫& 𝟎\end{array}\right]}$

• Note the block matrix structure with $𝑰$ the identity matrix, and

  ${\displaystyle 𝑫=\frac{1}{h^2}\left[\begin{array}{lllll}-2& 1& & & \\ 1& -2& 1& & \\ & \ddots & \ddots & \ddots & \\ & & 1& -2& 1\\ & & & 1& -2\end{array}\right],𝑰,𝑫\in {ℝ}^{(m-1)\times (m-1)}}$

• Apply leap-frog to the ODE system with time step $k$, ${𝒒}^n=𝒒(nk)$

  ${\displaystyle \frac{{𝒒}^{n+1}-{𝒒}^{n-1}}{2k}=𝑴{𝒒}^n}$

Stability considerations

• Leap-frog has a stability region $z=\lambda k$ on the slit from $z=-i$ to $z=+i$

• Eigenvalues of $𝑴$ are required. These can be determined analytically

  ${\displaystyle 𝑴𝒒=\lambda 𝒒\Rightarrow \left[\begin{array}{ll}𝟎& 𝑰\\ 𝑫& 𝟎\end{array}\right]\left[\begin{array}{l}𝒖\\ 𝒗\end{array}\right]=\lambda \left[\begin{array}{l}𝒖\\ 𝒗\end{array}\right]\Rightarrow \{\begin{array}{l}𝒗=\lambda 𝒖\\ 𝑫𝒖=\lambda 𝒗\end{array}\Rightarrow 𝑫𝒖={\lambda }^2𝒖.=\mu 𝒖}$

• The eigenvalues of $𝑫$ are

  ${\displaystyle {\mu }_l=-\frac{4}{h^2}{\sin }^2\left(\frac{lh}{2}\right),l=1,2,\dots ,m-1,h=\pi /m}$

• The eigenvalues of $𝑴$ are therefore

  ${\displaystyle {\lambda }_l=\sqrt{{\mu }_l}=\pm \frac{2i}{h}\sin \left(\frac{lh}{2}\right),l=1,2,\dots ,m-1}$

  and are purely imaginary, and therefore leap-frog numerical solutions can be made stable by an appropriate time step restriction.