Franz J. Vesely > CompPhys Tutorial > Partial Differential Equations  
 
 





 
< >
Part II: Ch. 5



5.1 Initial Value Problems I: Conservative-hyperbolic DE



$ \frac{\textstyle \partial u}{\textstyle \partial t} = - \frac{\textstyle \partial j}{\textstyle \partial x} $



Best (i.e. most stable, exact, etc.): Lax-Wendroff technique

Approach to Lax-Wendroff via:   FTCS $\Longrightarrow$ Lax $\Longrightarrow$ Leapfrog


Subsections


5.1.1 FTCS Scheme; Stability Analysis

Writing $u_{j}^{n} \equiv u(x_{j},t_{n})$ and using DNGF for the time derivative (FT, "forward-time"), and DST for the space derivative (CS, for "centered-space"), we write $\partial u/\partial t = - \partial j/\partial x$ as

$ \begin{eqnarray} \frac{\textstyle 1}{\textstyle \Delta t} \left[ u_{j}^{n+1}- u_{j}^{n} \right] & \approx & - \frac{\textstyle 1}{\textstyle 2 \Delta x} \left[ j_{j+1}^{n}- j_{j-1}^{n} \right] \end{eqnarray} $


$ \begin{eqnarray} && u_{j}^{n+1}= u_{j}^{n}-\frac{\textstyle \Delta t}{\textstyle 2 \Delta x} \left[ j_{j+1}^{n}- j_{j-1}^{n} \right] \end{eqnarray} $


This may be symbolized as follows:


Stability analysis (J. v. Neumann):

At time $t_{n}$, expand $u(x,t)$:
$ \begin{eqnarray} u_{j}^{n}&=& \sum_{k} U_{k}^{n} e^{\textstyle ikx_{j}} \end{eqnarray} $

where $ k=2\pi l/L \;\; (l=0,1,\dots)$. Insert this in $u_{j}^{n+1}=T[u_{j'}^{n}]$ to find each Fourier component's propagation law, $U_{k}^{n+1}=g(k) U_{k}^{n}$.

$\Longrightarrow$ Stable if $| g(k) | \leq 1 \;\;{\rm for \; all \; k}$.


Application to FTCS + advective equation with $j=cu$:

$ \begin{eqnarray} g(k) \, U_{k}^{n} e^{\textstyle ikj \Delta x} &=& U_{k}^{n} \, e^{\textstyle ikj \Delta x} - \frac{\textstyle c \Delta t}{\textstyle 2 \Delta x} U_{k}^{n} [e^{\textstyle ik(j+1)\Delta x} - e^{\textstyle ik(j-1)\Delta x}] \end{eqnarray} $

or
$ \begin{eqnarray} g(k) & = & 1- \frac{\textstyle ic \Delta t}{\textstyle \Delta x} \sin k \Delta x \end{eqnarray} $

Obviously, $\vert g(k)\vert > 1$ for any $k$; the FTCS method is inherently unstable.



5.1.2 Lax Scheme

Replacing in the FTCS formula the term $u_{j}^{n}$ by its spatial average $[u_{j+1}^{n}+ u_{j-1}^{n}]/2$, we approximate $\partial u/\partial t = - \partial j/\partial x$ by

$ \begin{eqnarray} && u_{j}^{n+1} = \frac{1}{2} \left[ u_{j+1}^{n}+ u_{j-1}^{n} \right] -\frac{\textstyle \Delta t}{\textstyle 2 \Delta x} \left[ j_{j+1}^{n}- j_{j-1}^{n} \right] \end{eqnarray} $



Stability / Friedrichs-Löwy condition:

Insert Fourier expanded $u(x)$ in Lax formula to find
$ \begin{eqnarray} g(k)= \cos k \Delta x -i \frac{\textstyle c \Delta t}{\textstyle \Delta x} \sin k \Delta x \end{eqnarray} $

The condition $\vert g(k)\vert \leq 1$ is tantamount to

$ \begin{eqnarray} && \frac{\textstyle |c| \Delta t}{\textstyle \Delta x} \leq 1 \end{eqnarray} $


Region below the dashed line: physically relevant for $u_{j}^{n+1}$, according to $x(t_{n+1})=x(t_{n})\pm |c| \Delta t$

Close scrutiny shows that LAX solves not the original PDE but

$ \begin{eqnarray} \frac{\textstyle \partial u}{\textstyle \partial t} & = & -c \frac{\textstyle \partial u}{\textstyle \partial x} + \frac{\textstyle (\Delta x)^{2}}{\textstyle 2 \Delta t} \frac{\textstyle \partial^{2} u}{\textstyle \partial x^{2}} \end{eqnarray} $

The additional diffusive term makes the method stable. However, it is an artefact and should be small:

$ \begin{eqnarray} |c| \Delta t & >> & \frac{\textstyle \Delta x}{\textstyle 2} \frac{\textstyle |\delta^{2} u|}{\textstyle |\delta u|} \end{eqnarray} $



5.1.3 Leapfrog Scheme (LF)

Use DST for $\partial/\partial t$: $ \partial u /\partial t \approx (u^{n+1}-u^{n-1})/2\Delta t$ to find the leapfrog expression

$ \begin{eqnarray} && u_{j}^{n+1}- u_{j}^{n-1}= -\frac{\textstyle \Delta t}{\textstyle \Delta x} \left[ j_{j+1}^{n}- j_{j-1}^{n} \right] \end{eqnarray} $


Stability requires once more that $c\Delta t/\Delta x \leq 1$ (CFL condition)



5.1.4 Lax-Wendroff Scheme (LW)


  • Lax method with half-step: $\Delta x/2$, $\Delta t/2$: $ \begin{eqnarray} u_{j+1/2}^{n+1/2} & = & \frac{1}{2} \left[ u_{j+1}^{n}+ u_{j}^{n} \right] -\frac{\textstyle \Delta t}{\textstyle 2\Delta x} \left[ j_{j+1}^{n}- j_{j}^{n} \right] \end{eqnarray} $ and analogously for $ u_{j-1/2}^{n+1/2}\;$.

  • Evaluation, e.g. for the advective case $j= C \cdot u$: $ \begin{eqnarray} u_{j+1/2}^{n+1/2} &\Longrightarrow& j_{j+1/2}^{n+1/2} \end{eqnarray} $

  • Leapfrog with half-step: $ \begin{eqnarray} u_{j}^{n+1}&=& u_{j}^{n} -\frac{\textstyle \Delta t}{\textstyle \Delta x} \left[ j_{j+1/2}^{n+1/2}- j_{j-1/2}^{n+1/2} \right] \end{eqnarray} $


Stability:

Once more assuming $j=cu$ and using the ansatz $U_{k}^{n+1}= g(k) U_{k}^{n}$ we find

$ \begin{eqnarray} g(k)&=&1-ia \sin k \Delta x -a^{2}(1-\cos k \Delta x), \end{eqnarray} $

with $a=c \Delta t/ \Delta x$. The requirement $\vert g\vert^{2}\leq 1$ leads once again to the CFL condition, $c\Delta t/\Delta x \leq 1$.




5.1.5 Lax and Lax-Wendroff in Two Dimensions


$ \frac{\textstyle \partial u}{\textstyle \partial t} = -\frac{\textstyle \partial j_{x}}{\textstyle \partial x} -\frac{\textstyle \partial j_{y}}{\textstyle \partial y} $

(advective case: $j_{x}=c_{x}u$ and $j_{y}=c_{y}u$)

Lax scheme:

$ u_{i,j}^{n+1} = \frac{1}{4} \left[ u_{i+1,\,j}^{n}+u_{i,j+1}^{n}+u_{i-1,j}^{n}+u_{i,j-1}^{n}\right] -\frac{\textstyle \Delta t}{\textstyle 2 \Delta x} \left[ j_{x,i+1,j}^{n}-j_{x,i-1,j}^{n}\right] - \frac{\textstyle \Delta t}{\textstyle 2 \Delta y} \left[j_{y,i,j+1}^{n}-j_{y,i,j-1}^{n}\right] $


Figure: Lax method in two dimensions


Lax-Wendroff:
For the second stage (half-step leapfrog) we need $j_{x,i+1/2,j-1/2}^{n+1/2}$ etc., which requires $u_{i+1/2,j-1/2}^{n+1/2}$, which must be determined from $u_{i,j-1/2}^{n},$ $u_{i+1,j-1/2}^{n}$ etc.

But: quantities with half-step spatial indices ($\scriptstyle i+1/2$, $\scriptstyle j-1/2$ etc.) are given at half-step times ($t_{n+1/2}$) only.

Modifying the LW scheme to allow for this, we have Lax-Wendroff in 2 dimensions:
  • Lax method to determine the $u$-values at half-step time $t_{n+1/2}$:

    $ \begin{eqnarray} u_{i+1,j}^{n+1/2}&=&\frac{1}{4} \left[ u_{i+2,j}^{n}+u_{i+1,j+1}^{n}+u_{i,j}^{n}+u_{i+1,j-1}^{n}\right] \\ && \\ && \;\;\;\; -\frac{\textstyle \Delta t}{\textstyle 2\Delta x} \left[ j_{x,i+2,j}^{n}-j_{x,i,j}^{n}\right] -\frac{\textstyle \Delta t}{\textstyle 2\Delta y} \left[ j_{y,i+1,j+1}^{n}-j_{y,i+1,j-1}^{n}\right] \end{eqnarray} $
    etc.

  • Evaluation at half-step time:

    $u_{i+1,j}^{n+1/2} , \dots \; \Longrightarrow \; j_{x,i+1,j}^{n+1/2} , \dots $

  • Leapfrog with half-step:

    $ \begin{eqnarray} u_{i,j}^{n+1}=u_{i,j}^{n} &-&\frac{\textstyle \Delta t}{\textstyle 2\Delta x} \left[ j_{yx,i+1,j}^{n+1/2}-j_{x,i-1,j}^{n+1/2} \right] - \frac{\textstyle \Delta t}{\textstyle 2\Delta y} \left[ j_{y,i,j+1}^{n+1/2}-j_{y,i,j-1}^{n+1/2} \right] \end{eqnarray} $

Lax-Wendroff in two dimensions

Figure: First stage (= Lax) in the 2-dimensional LW method: $ \textstyle \circ $ ... $t_{n},t_{n+1}$, ... $t_{n+1/2}$

For $u_{i,j}^{n+1}$ only the points $\textstyle \circ $ (at $t_{n}$) are used; for $u_{i+1,j}^{n+1}$ we use the points .

Problem: Drift between subgrids $\textstyle \circ $ and .

Solution: If the given PDE contains a diffusive term, this guarantees coupling. Otherwise, artificially add a small diffusive term.

Stability analysis: Fourier modes are now 2-dimensional:
$u(x,y)=\sum_{k} \sum_{l} U_{k,l} e^{\textstyle ikx+ily}$

Assuming $\Delta x = \Delta y$ we find the CFL condition
$ \Delta t \leq \frac{\textstyle \Delta x}{\textstyle \sqrt{2} \sqrt{c_{x}^{2}+c_{y}^{2}}} $




5.1.6 Resumé: Conservative-hyperbolic DE




$ \begin{eqnarray} \frac{\textstyle \partial u}{\textstyle \partial t} & = & - \frac{\textstyle \partial j}{\textstyle \partial x} \end{eqnarray} $


- Use Lax-Wendroff!

- If not, use at least Lax, but see that in addition to CFL

$ \begin{eqnarray} |c| \Delta t &>>& \frac{\textstyle \Delta x}{\textstyle 2} \frac{\textstyle |\delta^{2} u|}{\textstyle |\delta u|} \end{eqnarray} $

- Forget FTCS and Leapfrog!


To test the various methods, let us apply them to the 1D wave equation. When altering the propagation velocity $c$, the time step $\Delta t$, or the grid width $\Delta x$, keep in mind the operation regions of the different algorithms:

FTCS - unstable
LAX, LEAPFROG, LAX-WENDROFF - CFL condition $|c| \Delta t / \Delta x \leq 1$

Applet Hpde1: Start




vesely 2005-10-10

< >