Finite Difference Schemes

Most of this material comes from Rychtmyer and Morton's monograph and Strikwerda's textbook.

Take domain $(x,t)$ and define a lattice $(x_m, t_n)=(mh, nk)$, where $m$ and $n$ are integers. Typically, $n\ge 0$. We limit the presentation to grids which are uniform in both $x$ and $t$ (see Figure 22).

Figure 22: Space-time Lattice about $x_m,t_n$

here $\left\{\begin{array}{l} h\mbox { is }x\mbox { grid spacing}\\ k\mbox { is }t\mbox { grid spacing} \end{array}\right.$

Notation: Let $U^{n}_m=U(x_m, t_n)=U(mh, nk)$ be the value of $U$ on the lattice. Let $u^n_m$ be an approximation of $U_m^n$ on the same lattice location. We're going to consider mostly grids with constant grid spacing.

As in the ODE case, the most important properties of any numerical scheme for the approximation of a PDE (not just hyperbolic ones) are:

• Convergence
• Consistency
• Stability

Definition: Convergence: for one-step schemes approximating a ANY PDE to be convergent we compare $U(x,t)$ and $u^n_m$: if $U^0_m$ converges to $U_0(x)$ as $mh\rightarrow x$ then $U^n_m$ converges to $U(x,t)$ at $(m_h, nk)$ converges to $(x,t)$ as $h,k\rightarrow 0$. As $h,k\rightarrow 0$ the approximation gets uniformly closer to exact solution on the lattice.

Some properties of a scheme that we should be interested in are:

• Order of Accuracy
• Dissipation, Dispersion
• Speed and Efficiency

Of course, this list of properties is not exhaustive and the properties of importance are different for different types of PDE's.

A fundamental theorem of Finite Difference approximations of PDE's is the Lax-Richtmyer Equivalence Theorem

THE LAX-RICHTMYER EQUIVALENCE THEOREM

A consistent finite difference scheme for a PDE for which the initial value problem is well-posed is convergent if and only if it's stable.

Proof: See Chapter 10 Strikwerda's book.

So while proof of convergence would be a function-analytic exercise, we could instead check for consistency and stability and get convergence as a bonus. This is nice since stability and consistency is usually easier to check than convergence.

What's consistency?

definition:

Given a PDE of the form $PU=f$    and finite difference scheme $P_{k,h}u=f$, we say the FD scheme is consistent with the PDE if, for any smooth $\phi(t,x)$,

$$P\phi-P_{k_1,h}\phi\rightarrow 0$$$$\quad \mbox {as } h,k\rightarrow 0,$$

the convergence being pointwise convergence at each grid point.

Basic idea in Finite Difference Methods: replace derivatives by finite difference approximations. What we obtain is a pointwise approximation on a grid (no information on points not belonging to the lattice)

For the equation

$$U_t+a U_x=0$$

some schemes are

$$\mbox {\lq\lq forward space-forward time''}\quad \displaystyle \frac{u_m^{n+1}-u^n_m}{k}+a\frac{u^n_{m+1}-u^n_m}{h}=0$$

$$\mbox {\lq\lq forward time-centered space''}\quad \displaystyle \frac{u^{n+1}_m-u^n_m}{k}+a\frac{u^n_{m+1}-u^n_{m-1}}{2h}=0$$

$$\mbox {\lq\lq leapfrog''}\quad \displaystyle \frac{u_m^{n+1}-u^{n-1}_{m}}{2k}+a\frac{u^n_{m+1}-u^n_{m-1}}{2h}=0$$

$$\mbox {\lq\lq Lax-Friedrichs''}\quad \displaystyle \frac{u_m^{n+1}-\frac12\left(u^n_{m+1}+u^n_{m-1}\right)}{k}+a\frac {u^n_{m+1}-u^n_{m-1}}{2h}=0$$

and their computational cells are illustrated in the Figure 23.

Figure 23: Computational cells for some elementary schemes for the approximation of hyperbolic schemes

So $$u^{n}_m$$ is an approximation to $U(x,t)$ at $x=mh, t=mk$.

Assume that $U$ is sufficiently regular and continuous:

take $U(x_m\pm h, t)=U(x_m, t_n)\pm h\displaystyle \frac{\partial U}{\partial x}(x_m, t_n)+\displaystyle \frac12h^2\frac{d^2U}{dx^2}(x_m, t_n)+\cdots$

$$\Rightarrow \frac{U(x_m\pm h, t_n)-U(x_m, t_n)}{h}=\frac{\partial U}{dx}(x_m, t_n)+{\cal O}(h)$$

$$\Rightarrow \frac{(U_{m+h}, t_n)-U(x_m-h, t_n)}{2h}=\frac{\partial U}{dx}(x_m, t_n)+)+{\cal O} (h^2)$$

same procedure leads to finite difference approximation scheme to $$\frac{\partial }{\partial t}$$ derivatives.

Example check convergence of the Lax-Friedrichs scheme:

$$P_{k,h}\phi=\displaystyle \frac{\phi^{n+1}_m-\frac12\left(\phi^n_{m+1}+\phi^n_{m-1} \right)}{k}+a\frac{\phi^n_{m+1}-\phi^n_{m+1}}{2h}$$

$$\mbox {let } \phi^n_m\equiv \phi(t_n, x_m) \mbox { then }$$

$$\phi^n_{m\pm 1}=\phi^n_m\pm h\phi_x+\frac12\phi_{xx}\pm \frac{h^3}{6}\phi_{xxx}+{\cal O}(h^4)$$

$$\therefore\quad\frac12(\phi^n_{m+1}+\phi^n_{m-1})=\phi^n_m+\frac12h^2 \phi_{xx}+{\cal O}(h^4)$$

$$\frac{\phi^n_{m+1}-\phi^n_{m-1}}{2h}=\phi_x+\frac{1}{6}h^2\phi_{xxx}+{\cal O}(h^4)$$

substituting

$$P_{k,h}\phi=\phi_t+a\phi_x+\frac12k\phi_{tt}-\frac12\frac{h^2}{k}\phi_{xx}+\frac{1}{6}ah^2\phi_{xxx}+{\cal O}\left(h^4+\frac{h^4}{k}+k^2\right)$$

$$\mbox { so } P_{k,h}\phi\rightarrow 0 \mbox { as }k,h\rightarrow 0, \mbox { i.e. consistent as long as }\frac{h^2}{k}\rightarrow 0$$

so what happens when using finite $h$ and $k$? If $k$ is significantly smaller than $h^2$ then $$\frac{h^2}{k}\ge {\cal O}(1)$$ quantity. Thus, we would effectively be solving the problem $\phi_t+a\phi_x-L\phi_{xx}=0$ (due to finite truncation error terms) which is not what we set out to do in the first place!!

A Fundamental Theorem in FD approximations of Hyperbolic PDE's is the Courant-Friedricks-Lewy Condition (CFL), which will be related to the stability of a scheme.

Stability For the homogeneous problem $U_t+a U_x=f$, i.e. with $f=0:$

definition The IVP for the first order hyperbolic pde $U_t + a U_x =0$ is well-posed if for any time $T\ge t_0\quad \exists\quad C_T$ constant such that any solution $U(t,x)$ satisfies

$$\int^{\infty}_{-\infty}\vert U(t,x)\vert^2dx\le C_T\int^{\infty}_{-\infty}\vert U(t_0, x)\vert^2dx,$$$$\mbox {for }t_0\le t\le T$$

definition A finite difference scheme $P_{k,h}u^n_m=0$ for a first-order equation is stable in a stability region ${\mathbb{D}}$ if there's an integer $J$ such that for any positive time $T>t_0\quad \exists\quad C_T$ constant such that

$$h\sum^{\infty}_{m=-\infty}\vert u^n_m\vert^2\le C_Th\sum^J_{j=0}\sum^{\infty}_ {m=\infty}\vert u^j_m\vert^2$$

   for $$t_0\le t_0+nk\le T$$    with$$\quad (k,h)\in{\mathbb{D}}$$

$J=0$ for 1-step schemes and $J>0$ for multistep schemes, with data at first $J+1$ levels.

Example

Show that

$$\sum^{\infty}_{m=-\infty}\vert v^{n+1}_m\vert^2\le (\vert\alpha\vert+\vert\beta\vert)^2\sum^{\infty}_{m=-\infty}\vert v^n_m\vert^2$$$$\mbox { for}$$

$$\mbox {and thus the scheme } v^{n+1}_m=\alpha v^n_m+\beta v^n_{m+1}$$

$$\mbox {is stable if } (\vert\alpha\vert+\vert\beta\vert)\le 1.$$

$\Box$

The Courant-Friedricks-Lewy Condition (CFL)

definition: Explicit FD scheme can be written as

$$v^{n+1}_m=$$$$\mbox {a finite sum of }v^{n'}_{m'}\quad n'\le n$$

Theorem: For $U_t + a U_x =0$ with explicit scheme of the form $u_m^{n+1}= \alpha u^n_{m-1}+\beta u^n_m+\alpha u^n_{m+1}$ with $k/h=\lambda$ constant, a necessary and sufficient addition for stability is the CFL condition

$$\vert a_i\lambda\vert\le 1$$

For systems of equations, where $u$ is vector and $\alpha, \beta, \alpha$ are matrices, we require that $\vert a_i\lambda\vert\le 1$ for all e'values $a_i$ of the matrix $a$.

Proof: (Heuristic) See Figure 24

Figure 24: Heuristic proof of theorem, based on the finite time of travel of information from cell to cell.

Take Scalar Case: Take $\vert a\lambda\vert>1$ first: Consider $(t,x)=(1,0)$. The solution at $u(1,0)$ depends on value of $u_0(x)$ at either $x=a$ or $x=-a$ (depending on sign of speed ``a''). But from finite difference scheme we have

$$u^n_0$$$$\mbox {depends on } u^0_m \mbox{ only for }m\le n \mbox {( by the form of the scheme)}$$

Since $h=\lambda^{-1}k\Rightarrow mh\le k\lambda^{-1} n=\lambda^{-1}$ since $k=\displaystyle \frac{1}{n}$ (assumes that $t_{\mbox {final}}=1$), $\therefore u^n_0$ depends on $x$ only for $\vert x\vert\le\lambda^{-1}<\vert a\vert$. Thus $u^n_0$ cannot converge to $u(1,0)$ as $h\rightarrow 0$ with $$\frac{h}{k}=1$$. For $\vert a\lambda\vert\le 1$ things work out (convince yourself).

$\Box$

Theorem: (explicit schemes) There's no explicit, unconditionally stable, consistent finite difference schemes for hyperbolic systems of pde's.

Proof: Omitted. $\Box$

Remark: Unconditionally stable means that we can choose any $k$ and $h$ and still remain in the region of stability for the particular scheme. In many instances, a physical problem may require that we time-step an approximation over many many time steps. An explicit scheme is attractive here, because it is very efficient in storage (and usually easy to code). However, we need to consider how long a computation is going to actually take in real clock time: if we are restricted by a very small time step, then it may take a very long time to solve a problem. An alternative is to go to a higher order explicit scheme (but this usually means more communication which is of concern in parallel computing) and this buys us a little longer time steps. However, we might consider a low order implicit scheme which might buy us significantly bigger time steps (but usually a lot more communication). A recent popular alternative are the ``Semi-Lagrangian Methods''.

Example

$$\frac{u^{n+1}_m-u^n_m}{k}+a\frac{u^{n+1}-u^{n+1}_{m-1}}{h}=0$$

Implicit Case computational cell

$$(1+a\lambda)u^{n+1}_{m}=u^n_m+a\lambda u^{n+1}_{m-1}\quad \lambda= \quad \frac{k}{h}$$

square both sides

$$(1+a\lambda)^2\vert u^{n+1}_m\vert^2 \le \vert u^n_m\vert^2+2a\lambda\vert u^n_m\vert u^{n+1}_{m-1}\vert +a^2\lambda^2\vert u^{n+1}_{m-1}\vert^2$$

$$\le (1+a\lambda)\vert u^n_m\vert^2+a\lambda(1+a \lambda)\vert u^{n+1}_{m-1}\vert^2$$

Taking sums over all $m$:

$$(1+a\lambda)^2\sum^{\infty}_{m=-\infty}\vert u_m^{n+1}\vert^2\le ... ...^n_m\vert^2+a\lambda (1+a\lambda)\sum^{\infty}_{m=\infty}\vert u_m^{n+1}\vert^2$$

$$\therefore \sum^{\infty}_{m=-\infty}... ... u^{n+1}_m\vert^2\le \sum^{\infty}_{m=\infty}\vert u^n_m\vert^2\quad \therefore \quad \mbox {stable for all $\lambda$ with $a>0$.}$$

Example

Explicit Case:

$$\frac{u^{n+1}_{m}-u^n_m}{k}+a\frac{u^n_{m+1}-u^n_m}{h}=0$$

$$u^{n+1}_{m}=(1+a\lambda)u^n_m -a\lambda u^n_{m+1},\quad \lambda= \frac{k}{h}$$

square both sides:

$$\vert u_m^{n+1}\vert^2\le (1+\lambda a)^2\vert u^n_m\vert^2 -2a \lambda(1+a\lambda)\vert u^n_m\vert{u^n_{m+1}}\vert$$

$$+ (a\lambda)^2\vert u^n_{m+1}\vert^2$$

Take sum over $m$:

$$\sum^{\infty}_{m=-\infty}\vert u^{n+1}_m\vert^2\le \vert 1+2a\lambda\vert\sum^{\infty}_ {m=-\infty}\vert u^h_m\vert^2$$

require $\vert 1+2a\lambda\vert\le 1\ldots$ only if $a<0$ and for special $\lambda$