Section 2.2 Finite Differences Method: Hyperbolic PDEs
Before analyzing the algorithms to solve numerically the advection PDE, we illustrate a concrete physical problem from which this PDE arises. Let \(R\subset\Bbb R^3\) be a region bounded by a closed surface \(\partial R\) containing a mass of fluid \(M_R\) with density \(u\text{.}\) In particular, this means that the quantity of fluid contained in a subregion \(P\subset R\) is given by
\begin{equation*}
M_P = \iiint_P u(\mathbf x) dV.
\end{equation*}
Now, suppose that the fluid in some neighborhood of \(R\) moves with speed \({\mathbf v}(t,{\mathbf x})\text{.}\) Then we have the following continuity equation:
\begin{equation*}
\frac{dM(t)}{dt} = \Sigma(t) - \iint_{\partial R} u(t,{\mathbf x})\,{\mathbf v}(t,{\mathbf x})\cdot {\mathbf dS}({\mathbf x}),
\end{equation*}
where \(\Sigma\) is the contribution given by sources and sinks inside \(R\) and the last term represents the net amount of fluid leaving (or entering) the region \(R\text{:}\) for every point \({\mathbf x}\in\partial R\text{,}\) the amount, per time unit, that is entering/leaving \(R\) through \({\mathbf x}\) is equal to
\begin{equation*}
u(t,{\mathbf x})\,{\mathbf v}(t,{\mathbf x})\cdot {\mathbf dS}
\end{equation*}
(recall that \({\mathbf dS}\) at \({\mathbf x}\) is perpendicular to \(\partial R\) at that point). 
\begin{equation*}
\iint_{\partial V} {\mathbf F}\cdot {\mathbf{dS}} = \iiint_V {\mathbf\nabla}\cdot{\mathbf F} dV
\end{equation*}
to convert the integral continuity equation into its differential form:
\begin{equation*}
\frac{\partial u}{\partial t} + {\mathbf\nabla}\cdot\left(u {\mathbf v}\right) = 0.
\end{equation*}
Finally, in case of an incompressible fluid (i.e. \({\mathbf\nabla}\cdot{\mathbf v}=0\)) this equation reduces to the advection equation
\begin{equation*}
\frac{\partial u}{\partial t} + {\mathbf v} \cdot {\mathbf\nabla} u = 0.
\end{equation*}
1D Advection equation. We consider in particular the simple case of a 1-dimensional fluid with constant velocity field \(c\text{,}\) in which case the advection equation reduces to
\begin{equation*}
u_t+cu_x=0.
\end{equation*}
The most natural initial condition for such equation is specifying the shape of the solution at \(t=0\text{,}\) namely \(u(0,{\mathbf x})=f({\mathbf x})\text{.}\) 
\begin{equation*}
(\partial_t+c\partial_x) u(t,x) = 0,\; u(0,x)=f(x).
\end{equation*}
Hence \(u\) is constant on the integral lines of the vector field \(\partial_t+c\partial_x\text{.}\) These are the solutions of the ODE
\begin{equation*}
\begin{cases}
\dot t(s) = 1\cr
\dot x(s) = c\cr
\end{cases}
\end{equation*}
namely \(x(t)=x_0+ct\text{.}\) This means that \(u(x,t)\) is equal to the value of \(f\) at the point \((0,x_0)\) such that \(x_0+c t = x\text{,}\) namely \(x_0=x-ct\text{.}\) In other words, the solution is
\begin{equation*}
u(x,t) = f(x-ct),
\end{equation*}
namely the profile of the solution at time \(t\) is the same as the profile at the initial time shifted by \(ct\) to the right. Below we solve numerically the 1D advection equation with several different instances of the Finite Differences Method (Section 1.17) to test the weak and strong points of this method for hyperbolic PDEs. Recall that the basic idea of this method is replacing the "phase space" \([0,T]\times[0,L]\) with the uniform lattice
\begin{equation*}
t_n = n\Delta t,\;0\leq n\leq N;\;\;x_m = m\Delta x,\;0\leq m\leq M
\end{equation*}
and \(u_t\) and \(u_x\) with one of their discrete approximations. 
Definition 2.2.2.
The quantity
\begin{equation*}
c_0=\frac{\Delta x}{\Delta t}
\end{equation*}
has the dimension of a velocity. We call \(c_0\) the critical velocity of the lattice.Subsection 2.2.1 FTBS (Downwind Method)
. In this first method, we use the forward Euler approximation for time derivatives and the backward one for space derivatives (Forward in Time, Backward in Space):
\begin{equation*}
u_t(t_n,x_m)\to\frac{u_{n+1,m}-u_{n,m}}{\Delta t}
\end{equation*}
\begin{equation*}
u_x(t_n,x_m)\to\frac{u_{n,m}-u_{n,m-1}}{\Delta x}.
\end{equation*}
Then the PDE
\begin{equation*}
u_t+cu_x=0
\end{equation*}
becomes the linear system
\begin{equation*}
u_{n+1,m} = u_{n,m} - \nu\left(u_{n,m}-u_{n,m-1}\right),
\end{equation*}
with
\begin{equation*}
0\leq n\leq N-1,\;1\leq m\leq M,\;\nu=c/c_0.
\end{equation*}
FTBS Diagram. This scheme can be summarized by the following diagram: 
\begin{equation*}
u(0,x)=f(x)
\end{equation*}
becomes here the vector
\begin{equation*}
u_{_{0,0}}=f(0),\dots,u_{_{0,m}}=f(x_m),\dots,u_{_{0,M}}=f(L).
\end{equation*}
Notice that the solution of the FTBS system cannot give the values of the variables \(u_{n,0}\text{,}\) \(n\geq1\text{,}\) namely the values of the solution at the left endpoint at all times, since they would involve the coefficients \(u_{n-1,-1}\text{!}\) This is due to the fact that the PDE is defined for all \(x\in\Bbb R\) but numerical methods can only solve the problem within a finite box. The impossibility to get the \(u_{n,0}\text{,}\) \(n\geq1\text{,}\) corresponds to the fact that the characteristics of the PDE hit the boundary at all points \((t,0)\text{.}\) We need therefore to provide an extra boundary condition on the segment \((t,x)\text{,}\) \(0\leq x\leq L\text{,}\) corresponding to the values reached, on that segment, by the solution of the PDE.
Remark 1: this extra boundary condition is not really "extra". The point is that some of the characteristics meet first the time axis than the space axis (see Figure 2.2.1) and so the initial condition must be assiged also on that side. This is not a restriction in any way: the values on the x axis are anyway independent on each other (as long as the function is continuous).
Remark 2: different schemes will entail different kinds of extra boundary conditions!
FTBS local order. Let us set \(h=\Delta t\) and \(k=\Delta x\text{.}\) From the Taylor expansion formula we know that
\begin{equation*}
u(t+h,x) - u(t,x) = u_t(t,x)h + \frac{1}{2}u_{tt}(\tau,x)h^2,
\end{equation*}
\begin{equation*}
u(t,x+k) - u(t,x) = u_x(t,x)k + \frac{1}{2}u_{xx}(t,\xi)k^2.
\end{equation*}
Using these formulae in the method expression
\begin{equation*}
u_{n+1,m} - u_{n,m} - \nu\left(u_{n,m}-u_{n,m-1}\right) = 0
\end{equation*}
we get that the left expression is equal to
\begin{equation*}
u_t(t,x)h + \frac{1}{2}u_{tt}(\tau,x)h^2 - \nu \left(u_x(t,x)k + \frac{1}{2}u_{xx}(t,\xi)k^2\right).
\end{equation*}
Since, by hypohesis, \(u_t(t,x)h - \nu u_x(t,x)k=0\text{,}\) we see that the local error order is 2. Notice also that the second order error
\begin{equation*}
\frac{1}{2}u_{tt}(\tau,x)h^2 -\nu\frac{1}{2}u_{xx}(t,\xi)k^2
\end{equation*}
is equal to zero for \(\nu=1\text{,}\) so we get the best possible precision with this method precisely when \(\nu=1\text{.}\)
FTBS consistency: the CFL condition. Courant, Friedrichs and Lewy formulated the following simple necessary condition for the convergence of a finite difference scheme.
Definition 2.2.3.
For any variable \(u_{n,m}\) we define numerical and exact domains of dependence as follows: The numerical domain of dependence of \(u_{n,m}\) is the convex envelope of all coefficients \(u_{i,j}\) that \(u_{n,m}\) depends upon, explicitly or implicitly, in the linear system determined by the method. Its exact domain of dependence is the set of all points that lie on characteristic curves passing through \(u_{n,m}\text{.}\)
Theorem 2.2.4 (Courant, Friedrichs and Lewy).
A necessary condition for the convergence of a PDE Finite-Differences method is that the the exact domain of dependence of each variable be a subset of its numerical domain of dependence.Proof.
\begin{equation*}
x(t) = x_0-ct.
\end{equation*}
Hence, in case of the FTBS scheme, there cannot be convergence for \(c<0\) while, for \(c\geq0\text{,}\) we must also have \(c/c_0\leq1\text{,}\) namely \(0\leq\nu\leq1\text{.}\)
Remark: the CFL condition is necessary but not sufficient!
FTBS convergence: the von Neumann condition. Consider a solution of the form
\begin{equation*}
u_{n,m} = \lambda^n e^{ikm\Delta x}
\end{equation*}
(this is called a Fourier mode of the solution). Then, substituting this expression in the algorithm
\begin{equation*}
u_{n+1,m} = (1-\nu)u_{n,m} +\nu u_{n,m-1},
\end{equation*}
we obtain that
\begin{equation*}
\;\;\;\lambda^{n+1} e^{ikm\Delta x} = (1-\nu)\lambda^n e^{ikm\Delta x} + \nu\lambda^n e^{ik(m-1)\Delta x}\text{,}
\end{equation*}
namely
\begin{equation*}
\lambda = 1 - \nu + \nu e^{-ik\Delta x}.
\end{equation*}
The FTBS algorithm is stable when
\begin{equation*}
|\lambda|^2 = 1 - 4\nu(1-\nu)\sin^2\frac{k\Delta x}{2}\leq1.
\end{equation*}
This happens, independently on \(\Delta x\text{,}\) precisely when \(0\leq\nu\leq1\text{,}\) since
\begin{equation*}
\max_{0\leq\nu\leq1}4\nu(1-\nu)=1.
\end{equation*}
This is precisely the same range in which the method satisfies the CFL condition.
Implementation in Python. The code below uses FTBS to solve numerically the advection PDE
\begin{equation*}
u_t(t,x)+cu_x(t,x)=0,\;u(0,x)=u_0(x),\;(t,x)\in[0,1]\times[0,1],
\end{equation*}
whose exact solution is
\begin{equation*}
u(t,x) = u_0(x-ct).
\end{equation*}
The following python code produces an animation showing the time evolution of the solution. The output does not show in the web page so to see the animation you need to cut&paste the code and run it on your machine.
c=0.4, tSteps=600, xSteps=600. In this case \(\nu=0.4\text{.}\)
c=1, tSteps=600, xSteps=600. In this case, \(\nu=1\text{.}\)Subsection 2.2.2 FTFS (Upwind)
If we use forward differences also for space, we get the algorithm
\begin{equation*}
u_{n+1,m} = u_{n,m} - \nu\left(u_{n,m+1}-u_{n,m}\right),
\end{equation*}
Below we show the numerical domain of dependence. 
Ultimately, we showed that every Advection PDE can be solved numerically with either the FTBS algorithm or with the FTFS algorithm as long as we choose \(\Delta t\) and \(\Delta x\) so that \(|c|\leq c_0\text{.}\)
Subsection 2.2.3 FTCS (Centered)
Based on what seen above, a good candidate for an algorithm working for both \(c\geq0\) and \(c\geq0\) is the one where the centered differences approximation is used for the space derivative:
\begin{equation*}
u_x(t_n,x_m)\simeq\frac{u_{n,m+1}-u_{n,m-1}}{2\Delta x}.
\end{equation*}
The advection PDE \(u_t=cu_x\) becomes this way the linear system
\begin{equation*}
u_{n+1,m} = u_{n,m} - \frac{\nu}{2}\left(u_{n,m+1}-u_{n,m-1}\right)\,,
\end{equation*}
with \(0\leq n\leq N\text{,}\) \(1\leq j\leq M-1\text{,}\) \(\nu=c/c_0\text{.}\) Notice that for FTCS we need to set boundary conditions at both \(x=0\) and \(x=L\text{.}\)
FTCS Diagram. This scheme can be summarized by the following diagram: 
\begin{equation*}
|\lambda|^2=1+\nu^2\sin^2(k\Delta x)\geq1!
\end{equation*}
We say that the FTCS algorithm is unconditionally unstable, namely unstable for any value of \(\Delta t\) and \(\Delta x\text{.}\)Subsection 2.2.4 The Lax Method
It turns out that there is a simple way to fix FTCS's problem, namely replacing the \(u_{n,m}\) term in the FTCS system with the average over its two neighbours:
\begin{equation*}
u_{n+1,m} = \frac{u_{n,m+1}+u_{n,m-1}}{2} - \nu\frac{u_{n,m+1}-u_{n,m-1}}{2}\,.
\end{equation*}
The Lax Diagram This scheme can be summarized by the following diagram: 
\begin{equation*}
|\lambda|^2 = 1-(1-\nu^2)\sin^2(k\Delta x)
\end{equation*}
so the method is stable independently on \(\Delta x\) as long as \(\nu^2\leq1\) (conditional stability).
Implementation in Python. The code below uses the Lax algorithm to solve numerically the advection PDE
\begin{equation*}
u_t(t,x)+cu_x(t,x)=0,\;u(0,x)=u_0(x),\;(t,x)\in[0,1]\times[0,1],
\end{equation*}
whose exact solution is
\begin{equation*}
u(t,x) = u_0(x-ct).
\end{equation*}
The dispersion problem. The Lax method works nice from the point of view of the "phase", in the sense that the velocity of the profiles is close to the constant velocity of the exact solution. What is much less satisfying is the evident damping problem: the profiles tend to become lower and larger, namely the numerical algorithm has a diffusion effect, a qualitative behavior that is absent in the exact solution. A way to explain this phenomenon is to re-write the Lax algorithm as
\begin{equation*}
u_{n+1,m}-u_{n-1,m} + \nu(u_{n,m+1}-u_{n,m-1}) = (u_{n,m+1}-2u_{n,m}+u_{n-1,m}) - (u_{n+1,m}-2u_{n,m}+u_{n-1,m}).
\end{equation*}
After diving both sides by \(\Delta t\text{,}\) in the continuous limit this linear system becomes
\begin{equation*}
u_t + c u_x = \frac{(\Delta x)^2}{2\Delta t}u_{xx} - \frac{(\Delta t)^2}{2}u_{tt}
\end{equation*}
While the last term becomes negligible for \(\Delta t\to0\text{,}\) the combination \((\Delta x)^2/\Delta t\) does not necessarily do that. Hence, the Lax method is also a solver for a diffusion PDE (see the section about the heat PDE) and so diffusion is inherently inside the algorithm.
Subsection 2.2.5 The Lax-Wendroff Method
The methods we illustrated so far use all the same idea: given the exact solution \(\hat u(t,x)\) of the advection PDE \(u_t+cu_x=0\text{,}\) one expands \(\hat u\) in Taylor series in \(t\) about \(t_0\) and gets
\begin{equation*}
\hat u(t,x) = \hat u(t_0,x) + \Delta t\, \hat u_t(t_0,x) + O((\Delta t)^2),
\end{equation*}
where \(\Delta t=t-t_0\text{.}\) Then, after replacing \(\hat u_t\) with \(-c\hat u_x\) and using either forward or central differencing for the derivative with respect to \(x\text{,}\) one gets, respectively, the FTBS and FTCS algorithms (and the Lax one, which is a deformation of FTCS).
A natural way to improve the precision of the solution is using one more term of the Taylor expansion:
\begin{equation*}
\hat u(t,x) = \hat u(t_0,x) + \Delta t\, \hat u_t(t_0,x) + \frac{(\Delta t)^2}{2} \hat u_{tt}(t_0,x) + O((\Delta t)^3) =
\end{equation*}
\begin{equation*}
= \hat u(t_0,x) - c\Delta t\, \hat u_x(t_0,x) + c^2\frac{(\Delta t)^2}{2} \, \hat u_{xx}(t_0,x) + O((\Delta t)^3).
\end{equation*}
The Lax-Wendroff algorithms consists in dropping all terms beyond the second order and replacing the derivatives with the central differencing and backward differencing for the second derivative:
\begin{equation*}
u_{n+1,m} = u_{n,m} - \nu\frac{u_{n,m+1}-u_{n,m-1}}{2} + \nu^2\frac{u_{n,m+1}-2u_{n,m}+u_{n,m-1}}{2}
\end{equation*}
\begin{equation*}
= \alpha u_{n,m+1} + \beta u_{n,m} + \gamma u_{n,m-1},
\end{equation*}
where
\begin{equation*}
\alpha = -\frac{\nu(1-\nu)}{2},\; \beta = 1-\nu^2,\; \gamma = \frac{\nu(1+\nu)}{2}.
\end{equation*}
Lax-Wendroff CFL and Von Neumann conditions. The CFL condition for this case is the same as for the Lax method, so a necessary condition for its consistence is that \(|c|\leq c_0\text{.}\) Moreover, in this case
\begin{equation*}
|\lambda|^2 = 1-\nu^2(1-\nu^2)(1-\cos(k\Delta x))^2
\end{equation*}
so, just like the Lax method, it is stable independently on \(\Delta x\) as long as \(\nu^2\leq1\text{.}\) Unlike the Lax method, though, the diffusion in the Lax-Wendroff algorithm is much weaker.
Implementation in Python. The code below uses the Lax algorithm to solve numerically the advection PDE
\begin{equation*}
u_t(t,x)+cu_x(t,x)=0,\;u(0,x)=u_0(x),\;(t,x)\in[0,1]\times[0,1],
\end{equation*}
whose exact solution is
\begin{equation*}
u(t,x) = u_0(x-ct).
\end{equation*}