Discretization

The continuous domain \([a,b]\times [0,\infty)\) can be replaced into finite number of points \((x_i, t_n)\) as follows:

Space discretization: spatial domain \([a, b]\) into \(M+2\) equidistant points

\[x_{i}=a+i \Delta x, \ i=0,...,M+1, \ \Delta x=\frac{b-a}{M-1}\]
  • There are \(M\) interior points \((x_1, x_2, ..., x_M)\).
  • Two boundary points exist \((x_0=a, x_{M+1}=b)\).
  • \(\Delta x\) present mesh size.

Time discretization: temporal domain \([0, T_{out}]\) into a set of time points,

\[t_n = n\Delta t, \ n=0,...,N_{out},...\]
  • \(t_{0}=0\) is initial time.
  • \(T_{out}=\Delta tN_{out}\), where \(\Delta t\) is time step.

In the relation \(c = \lambda\frac{\Delta t}{\Delta x}\), if \(c\) and \(\lambda\) are constants, it implies a fixed relationship between \(\Delta t\) and \(\Delta x\) as follows:

\[\Delta x = \Delta t \times K, \ K>0 \ (\text{constant})\]

Using these points, the continuous function \(q(x,t)\) can be replaced by a finite number of discrete values \(q(x_i, t_n)\). Subsequently, with these discrete values, the partial differential equation (PDE) can be solved approximately. This method is known as finite difference approximations, where the differential operator is replaced by a numerical operator.

Approximation of Derivatives

For a smooth function \(g(z)\), all the derivatives \(g^{(k)}(z_0)\) at \(z_0\) can be obtained using Taylor’s theorem as follows:

\[g(z_0+\Delta z)=g(z_0)+\Delta zg'(z_0)+\frac{1}{2}\Delta z^2g^{(2)}(z_0)+O(\Delta z^3)\] \[g(z_0+\Delta z)-g(z_0)=\Delta zg'(z_0)+\frac{1}{2}\Delta z^2g^{(2)}(z_0)+O(\Delta z^3)\] \[g'(z_0)=\frac{g(z_0+\Delta z) - g(z_0)}{\Delta z}+O(\Delta z)\] \[g'(z_0)\approx \frac{g(z_0+\Delta z) - g(z_0)}{\Delta z}\]

This is a first-order forward approximation. For first-oder backward approximation,

\[g(z_0-\Delta z)=g(z_0)-\Delta zg'(z_0)+\frac{1}{2}\Delta z^2g^{(2)}(z_0)+O(\Delta z^3)\] \[g'(z_0)=\frac{g(z_0) - g(z_0-\Delta z)}{\Delta z}+O(\Delta z).\]

For second-order central approximation,

\[g'(z_0)=\frac{g(z_0+\Delta z) - g(z_0-\Delta z)}{2\Delta z}+O(\Delta z^2).\]

Finite Difference Method (FDM)

Finite difference method approximates the exact value \(q(x_i, t_n)\) at the mesh point \((x_i, t_n)\) with an approximated value \(q_i^n\) as follows:

\[q_i^n\approx q(x_i, t_n)\]

Using forward approximation, the temporal and spatial derivatives in the PDE can be obtained by

\[\frac{\partial q(x_i, t_n)}{\partial t}=\frac{q(x_i, t_{n+1})-q(x_i, t_n)}{\Delta t}+O(\Delta t),\] \[\frac{\partial q(x_i, t_n)}{\partial x}=\frac{q(x_{i+1}, t_{n}) - q(x_i, t_n)}{\Delta x}+O(\Delta x).\]

Various Methods

Godunov’s Method

For Godunov’s method, time derivatives is given by

\[\frac{\partial q(x_i, t_n)}{\partial t}=\frac{q(x_i, t_{n+1})-q(x_i, t_n)}{\Delta t}+O(\Delta t).\]

However, space derivatives depend on the characteristic speed as follows:

\[\frac{\partial q(x_i, t_n)}{\partial x}=\frac{q(x_{i}, t_{n}) - q(x_{i-1}, t_n)}{\Delta x}+O(\Delta x),\quad \text{if}\quad \lambda>0.\] \[\frac{\partial q(x_i, t_n)}{\partial x}=\frac{q(x_{i+1}, t_{n}) - q(x_i, t_n)}{\Delta x}+O(\Delta x),\quad \text{if}\quad \lambda<0.\]

These derivatives (e.g., \(\lambda > 0\)) are substituted into the differential operator in PDE as follows:

\[L_e(q)\equiv\frac{\partial q(x,t)}{\partial t}+\lambda\frac{\partial q(x,t)}{\partial x}=0\] \[L_e(q(x_i,t_n)) = \frac{\partial q(x_i,t_n)}{\partial t}+\lambda\frac{\partial q(x_i,t_n)}{\partial x}\] \[\begin{aligned} L_e\left(q\left(x_i, t_n\right)\right) & =\partial_t q\left(x_i, t_n\right)+\lambda \partial_x q\left(x_i, t_n\right) \\ & =\frac{q\left(x_i, t_{n+1}\right)-q\left(x_i, t_n\right)}{\Delta t}+O(\Delta t) \\ & =\lambda \frac{q\left(x_i, t_n\right)-q\left(x_{i-1}, t_n\right)}{\Delta x}+O(\Delta x) \\ & =0 \end{aligned}\]

Here, the exact solution, \(q(x_i, t_n)\) is replaced by the approximated value, \(q_i^n\), suppressing terms of \(O(\Delta t) + O(\Delta x)\), as follows:

\[\frac{q_i^{n+1}-q_i^n}{\Delta t}+\lambda\left(\frac{q_i^n-q_{i-1}^n}{\Delta x}\right)=0\]

This consists of a numerical operator instead of a differential operator and can be rewritten by

\[q_i^{n+1}=q_i^{n}-c\left(q_i^n-q_{i-1}^n\right)\]

where \(c=\frac{\lambda \Delta t}{\Delta x}\) means CFL number. The equation states that the Godunov scheme is always upwind.

FTCS method

FTCS method means Forward-in-Time Centred-in-Space.

\[\frac{\partial q(x_i, t_n)}{\partial t}=\frac{q(x_i, t_{n+1})-q(x_i, t_n)}{\Delta t}+O(\Delta t).\] \[\frac{\partial q(x_i, t_n)}{\partial x}=\frac{q(x_{i+1}, t_{n})-q(x_{i-1}, t_n)}{2\Delta t}+O(\Delta x^2).\]

The exact solution, \(q(x_i, t_n)\) is replaced by the approximated value, \(q_i^n\), suppressing terms of \(O(\Delta t) + O(\Delta x^2)\), as follows:

\[\frac{q_i^{n+1}-q_i^n}{\Delta t}+\lambda\left(\frac{q_{i+1}^n-q_{i-1}^n}{2\Delta x}\right)=0\] \[q_i^{n+1}=q_i^{n}-\frac{1}{2}c\left(q_{i+1}^n-q_{i-1}^n\right)\]

FTCS seems better than Godunov’s method due to its second-order accuracy in space; however, unfortunately, since FTCS is unconditionally unstable, it is useless. To rescue this, there are serveral approaches as follows:

Explicit Lax-Friedrichs (LF) Scheme

LF method replaces \(q_i^n\) with half of the neighboring values \(q_{i-1}^n\) and \(q_{i+1}^n\). Thus, FTCS method can be written as

\[\frac{q_i^{n+1}-\frac{1}{2}(q_{i+1}^n+q_{i-1}^n)}{\Delta t}+\lambda\left(\frac{q_{i+1}^n-q_{i-1}^n}{2\Delta x}\right)=0\] \[q_i^{n+1}=\frac{1}{2}\left(1+c\right)q_{i-1}^n+\frac{1}{2}\left(1-c\right)q_{i+1}^n\]

For example, if the CFL number is positive (i.e., the characteristic speed is positive), the first term on the right side of the equation have a greater impact on the solution than the second term. Conversely, if the CFL number is negative, the second term has a more pronounced effect on the solution.

Lax-Wendroff (LW) Scheme

LW method employs Cauchy-Kowalewskaya procedure, which replaces time derivatives with space derivatives as follow:

\[\frac{\partial q}{\partial t}=-\lambda\frac{\partial q}{\partial x},\qquad \frac{\partial^2 q}{\partial t^2}=\lambda^2\frac{\partial^2 q}{\partial x^2},\qquad \frac{\partial^{(k)} q}{\partial t^{(k)}}=\left(-\lambda^{k}\right)\frac{\partial^{(k)} q}{\partial x^{(k)}}\]

Using these, a Taylor series in time can be written as

\[q(x_i,t_{n+1})=q(x_i,t_{n})-\Delta t\lambda\partial_x q(x_i,t_{n})+\frac{1}{2}\Delta^2\lambda^2\partial^2_xq(x_i,t_{n})+O(\Delta t^3)\] \[\partial_x q(x_i,t_{n})=\frac{q(x_{i+1},t_n)-q(x_{i-1},t_n)}{2\Delta x}+O(\Delta x^2)\] \[\partial_x^2 q(x_i,t_{n})=\frac{q(x_{i+1},t_n)-2q(x_{i},t_n)+q(x_{i-1},t_n)}{\Delta x^2}+O(\Delta x^2)\] \[q^{n+1}_i=\frac{1}{2}c(1+c)q^{n}_{i-1}+(1-c^2)q^{n}_{i}-\frac{1}{2}c(1-c)q^{n}_{i+1}\]

This equation is second-order accuracy and stable in \(\left\vert c \right\vert\le1\)

Implicit Time Discretization

Implicit form can be obtained following manner:

\[\partial_xq(x_i, t_n)=\frac{q(x_{i+1}, t_{n+1})-q(x_{i-1}, t_{n+1})}{2\Delta x}+O(\Delta x^2)\] \[\frac{q^{n+1}_i-q^{n}_i}{\Delta t}+\lambda\left(\frac{q^{n+1}_{i+1}-q^{n+1}_{i-1}}{2\Delta x}\right)=0\] \[q^{n+1}_{i}=q^n_i-\frac{1}{2}c\left(q^{n+1}_{i+1}-q^{n+1}_{i-1}\right)\] \[-\frac{1}{2}cq^{n+1}_{i-1}+q^{n+1}_{i}+\frac{1}{2}cq^{n+1}_{i+1}=q^n_i\]

From this, the linear algebraic system can be obtained as

\[\bf{A}\bf{X}=\bf{B}\]

where

\[\bf{A}=\left[\begin{array}{lll} -\frac{1}{2} c & 1 & \frac{1}{2} c \end{array}\right],\qquad \bf{X}=\left[\begin{array}{lll} q^{n+1}_{i-1} & q^{n+1}_{i} & q^{n+1}_{i+1} \end{array}\right]^{\text{T}},\qquad \bf{B}=q^{n}_{i}\]

Conservation Form

The general scalar conservation law (=balance law) is given by

\[\frac{\partial q}{\partial t}+\frac{\partial f(q)}{\partial x}=0\]

where \(f(q)\) is the physical flux. This can be rewritten by

\[q^{n+1}_i=q^n_i-\frac{\Delta t}{\Delta x}\left(f_{i+\frac{1}{2}}-f_{i-\frac{1}{2}}\right)\]

where \(f_{i+\frac{1}{2}}=f_{i+\frac{1}{2}}\left(q^n_{i-l},...,q^n_{i},q^n_{i+1},...,q^n_{i+r}\right)\) represents the numerical flux, requiring to satisfy

  • Continuity: \(f_{i+\frac{1}{2}} \left(\mathbb{R}^{l+r}\rightarrow\mathbb{R}\right)\) is a continuous real-value function
  • Consistency: when all arguments of numerical flux are identical, numerical flux become a physical flux, \(f_{i+\frac{1}{2}}(\hat{q},...,\hat{q})=f(\hat{q})\).

Lax-Wendroff (LW) Scheme

\[q^\text{LW}_{i+\frac{1}{2}}=\frac{1}{2}(q^n_i+q^n_{i+1})-\frac{1}{2}\frac{\Delta t}{\Delta x}\left[f(q_{i+1}^n)-f(q_{i}^n)\right]\] \[f^\text{LW}_{i+\frac{1}{2}}=f(q^\text{LW}_{i+\frac{1}{2}})\]

Godunov Centred (GC) Scheme

\[q^\text{GC}_{i+\frac{1}{2}}=\frac{1}{2}(q^n_i+q^n_{i+1})-\frac{\Delta t}{\Delta x}\left[f(q_{i+1}^n)-f(q_{i}^n)\right]\] \[f^\text{GC}_{i+\frac{1}{2}}=f(q^\text{GC}_{i+\frac{1}{2}})\]

Lax-Friedrichs (LF) Scheme

\[f^\text{LF}_{i+\frac{1}{2}}=\frac{1}{2}\left[f(q^n_i)+f(q^n_{i+1})\right]-\frac{1}{2}\frac{\Delta x}{\Delta t}(q_{i+1}^n-q_{i}^n)\]

FORCE Scheme

\[f^\text{FORCE}_{i+\frac{1}{2}}=\frac{1}{2}\left[f^\text{LF}_{i+\frac{1}{2}}+f^\text{LW}_{i+\frac{1}{2}}\right]\]

Time Step, \(\Delta t\)

The choice of \(\Delta t\) for explicit methods is constraint by a stability condition (i.e., CFL condition).

From definition of CFL number,

\[C_{cfl}=\lambda \frac{\Delta t}{\Delta x}\] \[\Delta t=C_\text{cfl}\frac{\Delta x}{\lambda}\]

where \(\lambda\) have to be the largest wave speed present at time level \(n\) \(\left(\lambda = S^n_\text{max}\right)\) and CFL coefficient \((C_\text{cfl})\) is limited as follows:

\[0<C_\text{cfl}\le C_\text{lim}.\]

Here, \(C_\text{lim}\) is linearized stability limit of the scheme. Practically, CFL coefficient is given by

\[C_\text{cfl}=0.9\times C_\text{lim}.\]