Basic Notion

General Form

In the previous post, I studied various numerical methods to solve the partial difference equation. These finite difference methods can be written in the general from as follows:

\[q_i^{n+1}=H\left(q^n_{i-l},q^n_{i-l+1},...,q^n_{i},...,q^n_{i+r}\right),\]

where \(H\) is a real-valued function of \(l+r+1\) real variables and depend on the specific numerical method. For example, Godunov’s method is expressed in the general form as follows:

For \(\lambda>0\),

\[H=cq^n_{n-1}+\left(1-c\right)q^n_{i},\]

For \(\lambda<0\),

\[H=\left(1+c\right)q^n_{n}-cq^n_{i+1}.\]

Monotone Scheme

The scheme is monotone (i.e., oscillation-free), if satisfying

\[\frac{\partial}{\partial q^n_k}H\left(q^n_{i-l},q^n_{i-l+1},...,q^n_{i},...,q^n_{i+r}\right)\ge 0,\]

which can apply to non-linear as well as linear scheme. This means the derivatives with respect to all arguments \(\left(i-l\le k \le i+r\right)\) have to be non-negative.

Linear Scheme

Linear schemes are defined as the special class of general form for the linear advection equation as follows:

\[q_i^{n+1}=\sum^{k=r}_{k=-l}b_kq^n_{i+k}\]

where the coefficients \(b_k\) are constant (i.e., independent on the solution \(q\)).

It is noted that linear scheme is monotone if and only if all its coefficients are non-negative.

Godunov’s Theorem

There are no monotone, linear schemes for the linear advection equation with constant \(\lambda\), of accuracy two or higher. In other words, linear monotone schemes are at most first-order accurate, thus necessary condition for a numerical scheme to be oscillation-free and of high-order of accuracy is to be non-linear \((b_k=b_k(q))\). In other words, the scheme must be non-linear, even when applied to linear equations.

\(L^\infty\)-Stability

Scheme is \(L^\infty\)-stability, if for all given mesh point \(i\) satisfying

\[C_m \le q^n_i \le C_M\]

the solution for all mesh point (\(i\)) satisfies

\[C_m \le q^{n+1}_{i} \le C_M\]

This means that if the initial data is bounded below and above by constants \(C_m\) and \(C_M\) respectively, then the solution at the nex time level is also bounded below and above by the same constants \(C_m\) and \(C_M\).

Local Truncation Error (LTE)

The general form is associated numerical operator as follows:

\[L_a\left(q^n_i\right)=q_i^{n+1}-H\left(q^n_{i-l},q^n_{i-l+1},...,q^n_{i},...,q^n_{i+r}\right)=0\]

where \(a\) means the ‘approximate’ and \(q^n_i\) is the approximated value. This numerical operator is different from the differential operator given by

\[L_e\left(q\right)=\partial_t q+ \lambda\partial_x q=0\]

where \(e\) denotes the ‘exact’ and \(q\) is the exact value.

To evaluate error of the numerical operator, the local truncation error (LTE) is defined as

\[\tau^n_i=\frac{1}{\Delta t}L_a\left(q\left(x_i,t_n\right)\right)\]

where \(\Delta t\) has a role of scaling factor. It is noted that because the numerical operator \(L_a\) is applied to the exact solution \(q\) evaluated at the mesh point \((x_i,t_n)\), the result is no longer zero.

To calculate the local truncation error, following steps are performed.

  • First of all, it assume that the solution \(q(x,t)\) of the PDE is sufficiently smooth (for use of Taylor expansions).
  • Taylor expansions in space and time about the point \((x_i, t_n)\) give
\[\begin{aligned} \tau^n_i&=\partial_tq\left(i\Delta x,n\Delta t\right)+\lambda \partial_x q\left(i\Delta x,n\Delta t\right) \\ &+O(\Delta t^k) + O(\Delta x^m) \end{aligned}\]

where \(k\) and \(m\) are positive integers. The first line is identically same as differential operator, therefore, remaining terms are given by

\[\tau^n_i=O(\Delta t^k) + O(\Delta x^m)\]

where the right terms denotes the leading term of \(\tau^n_i\).

Order of Accuracy

From this expression, the numerical scheme has order of accuracy \(k\) in time and \(m\) in space. When the fixed relationship between \(\Delta x\) and \(\Delta t\) is considered, the order of accuracy is defined as

\[p=\min\{k,m\}\]

Consistency

If \(\tau_i^n \rightarrow0\) as \(\Delta x \rightarrow0\) (or \(\Delta t \rightarrow0\)), the numerical scheme is said to be consistent with the partial differential equation.

Example: Godunov Scheme

To understand the calculation procedure of local truncation error, the Godunov scheme (in case of \(\lambda <0\)) is given by

\[H=\left(1+c\right)q^n_{n}-cq^n_{i+1}\]

Therefore, numerical operator and local truncation error can be expressed by

\[\begin{aligned} L_a\left(q_i^{n}\right)&=q_i^{n+1}-H\\ &=q_i^{n+1}-\left(1+c\right)q^n_{n}+cq^n_{i+1}=0 \end{aligned}\] \[\begin{aligned} \tau^n_i &= \frac{1}{\Delta t}L_a\left(q\left(x_i,t_n\right)\right) \\ &= \frac{1}{\Delta t}\left[q\left(i\Delta x, (n+1)\Delta t\right)-(1+c)q\left(i\Delta x, \Delta t\right) + cq\left((i+1)\Delta x, n\Delta t\right)\right] \end{aligned}\]

by Taylor expansions at the point \(\left(x_i, t_n\right)\),

\[\begin{aligned} q\left(i\Delta x, (n+1)\Delta t\right) &= q\left(i\Delta x, n\Delta t\right) + \Delta t\partial_t q\left(i\Delta x, n\Delta t\right) \\ &+ \frac{1}{2}\Delta t^2\partial_t^2q\left(i\Delta x, n\Delta t\right) + O(\Delta t^3), \end{aligned}\] \[\begin{aligned} q\left((i+1)\Delta x, \Delta t\right) &= q\left(i\Delta x, n\Delta t\right) + \Delta x\partial_x q\left(i\Delta x, n\Delta t\right) \\ &+ \frac{1}{2}\Delta x^2\partial_x^2q\left(i\Delta x, n\Delta t\right) + O(\Delta x^3), \end{aligned}\]

Finally, the local truncation error is given by

\[\begin{aligned} \tau^n_i&=\frac{1}{\Delta t}[q\left(i\Delta x, n\Delta t\right) + \Delta t\partial_t q\left(i\Delta x, n\Delta t\right) \\ &+ \frac{1}{2}\Delta t^2\partial_t^2q\left(i\Delta x, n\Delta t\right) + O(\Delta t^3) \\ &- (1+c)q\left(i\Delta x, n\Delta t\right) \\ & +c \{ q\left(i\Delta x, n\Delta t\right) + \Delta x\partial_x q\left(i\Delta x, n\Delta t\right) \\ & + \frac{1}{2}\Delta x^2\partial_x^2q\left(i\Delta x, n\Delta t\right) + O(\Delta x^3) \}] \\ \\ &= \partial_t q(i\Delta x, n\Delta t)+\lambda\partial_x q(i\Delta x, n\Delta t) \\ &+ \frac{1}{2}\Delta t \partial_t^2q(i\Delta x, n\Delta t)+\frac{1}{2} \lambda \Delta x \partial_x^2q(i\Delta x, n\Delta t) \\ &+ O(\Delta t^2) + O(\Delta x^2) \end{aligned}\]

Because the first term is exact differential equation, this term goes to zero. Then, the leading term of the local truncation error is given by

\[\tau^n_i=\frac{1}{2}\Delta t \partial_t^2q(i\Delta x, n\Delta t)+\frac{1}{2} \lambda \Delta x \partial_x^2q(i\Delta x, n\Delta t)=O(\Delta t) + O(\delta x)\]

Therefore, the Godunov scheme is the first-order accurate in space in time. Moreover, this equation shows that the local truncation error depend on the mesh parameter \(\left(\Delta t, \Delta x\right)\) and due to the derivatives \(\left(\partial^2\right)\), the local truncation error increase at the discontinuity region.

Using Cauchy-Kowalewskaya procedure, time derivatives can be converted to space derivatives as follows:

\[\tau^n_i = \frac{1}{2}\lambda\Delta x \left(1+c\right)\partial_x^2q(i\Delta x, n\Delta t)\]

where the term of \(\partial_x^2q\) means the diffusion term and \(\alpha = \frac{1}{2}\lambda\Delta x \left(1+c\right)\) is the diffusion coefficient, which shows that

  • when \(\Delta x\) increase, the diffusion coefficient raise, so the numerical diffusion become severe.
  • due to this numerical diffusion term, the first order method cause the high diffusive solution at the discontinuity region.
  • when \(c\) goes to \(-1\), the local truncation error decrease, thus the accuracy increase.

Modified Equation

In the local truncation error equation,

\[\begin{aligned} \tau^n_i &= \partial_t q(i\Delta x, n\Delta t)+\lambda\partial_x q(i\Delta x, n\Delta t) \\ &+ \frac{1}{2}\Delta t \partial_t^2q(i\Delta x, n\Delta t)+\frac{1}{2} \lambda \Delta x \partial_x^2q(i\Delta x, n\Delta t) \\ &+ O(\Delta t^2) + O(\Delta x^2) \end{aligned}\]

Originally, the order of accuracy is first order because the \(q(x,t)\) is the solution of the original PDE (here, the first line). However, if the \(q(x,t)\) is the solution of the following PDE, more accurate solution can be obtained because the order of accuracy become second order.

\[\partial_t q+\lambda\partial_x q + \frac{1}{2}\left[\Delta t \partial_t^2q+\lambda \Delta x \partial_x^2q\right] = 0\]

This equation called the modified equation. Here, if the Cauchy-Kovalevskaya procedure is applied to the modified equation to convert the second order time derivative to a space derivative, the final form of the modified equation is given by

\[\partial_t q+\lambda\partial_x q =\alpha\partial_x^2q\]

where \(\alpha=-\frac{1}{2}\lambda\Delta x \left(1+c\right)\) called the coefficient of numerical viscosity. This equation is a second-order PDE of the advection-diffusion type (parabolic), which shows that

  • due to the diffusion term (i.e., numerical viscous term), not physical viscosity, even the solutions of the inviscid equation behave as solutions of a viscous equation.
  • First-order scheme have modified equation that are of the advection-diffusion type with numerical viscosity (second-order)

  • Second-order scheme have modified equation of the dispersion type (third-order) as follows: \(\partial_t q+\lambda\partial_x q = \gamma \partial_x^3q\) where \(\gamma\) is coefficient of numerical dispersion. this dispersion error show up in the form of phase errors. For example, Lax-Wendroff method has this modified equation, \(\partial_t q+\lambda\partial_x q = \gamma \partial_x^3q;\quad\gamma=\frac{1}{6}\lambda\Delta x^2\left(c^2-1\right).\)

    Because \(\gamma\) depend on the power of grid size, when \(\Delta x\) decrease, \(\gamma\) rapidly decrease. This means that this scheme has a high convergence speed to exact solution.

It is noted that when dealing with physically viscosity situation, both numerical and physical viscosity coexist. Depending on the grid size, numerical viscosity can have higher values than physical viscosity, so it is essential to handle these two terms carefully.

Linear Stability Analysis

The von Neumann method for linear stability analysis of a numerical scheme introduces a trial function (i.e., test function) as follows:

\[q^n_i = A^n\exp\left(I\theta i\right)\]

where \(A\) represent an amplitude, \(A^n\) means \(A\) raised to the power \(n\) (not time notation), the complex unit is denoted by \(I\equiv\sqrt{-1}\), \(\theta=P\Delta x\) is a phase angle, and \(P\) is the wave number in the \(x\)-direction.

Representing a numerical solution with a trail function means that as the solution advances in time \((q^n, n\uparrow)\), the amplitude increases \((A^n, n\uparrow)\). Therefore, for stability, the stability condition has to be satisfied as follows:

\[\vert A \vert \le 1\]

Example: Godunov Scheme

In case of \(\lambda >0\),

\[q^{n+1}_i=q^n_i-c\left(q_i^n-q^n_{i-1}\right)\]

Using these notation,

\[\begin{aligned} q^{n+1}_i&=A^{n+1}\exp\left(I\theta i\right) \\ q^n_i &= A^{n}\exp\left(I\theta i\right)\\ q^n_{i-1} &= A^{n}\exp\left(I\theta (i-1)\right), \end{aligned}\] \[A^{n+1}\exp\left(I\theta i\right)=A^{n}\exp\left(I\theta i\right) - c\left[A^{n}\exp\left(I\theta i\right)-A^{n}\exp\left(I\theta (i-1)\right)\right]\]

After cancelling \(A^{n}\exp\left(I\theta i\right)\) and some mathematical manipulations, the amplitude can be obtained as follows:

\[\Vert A\Vert^2=(1-c)^2+2c(1-c)cos\theta+c^2\]

Therefore, by imposing \(\Vert A\Vert^2\le 1\), the stable condition for the Godunov scheme is given by \(c\le0\le1\).

Time Step

Now, the time step can be obtained from stability condition. From the condition \((c=\frac{\lambda \Delta t}{\Delta x}\le1)\), the time step is given by

\[\begin{aligned} \Delta t &\le \frac{\Delta x}{\lambda}, \\ \Delta t &= C_{cfl}\frac{\Delta x}{\lambda}, \\ \end{aligned}\] \[\left(0 < C_{cfl} \le C_{lim}\right).\]

where \(C_{lim}\) is the linear stability limit and \(C_{cfl}\) is the CFL number, which is typically defined as \(C_{cfl}=0.9C_{lim}\).

Shortcut

For the high-order schemes, applying the linear stability analysis mentioned previously become very complicated. To overcome it, the shortcut for stability is as follows:

Shortcut to Stability

If a numerical scheme can be expressed in three-point linear scheme as follows:

\[q^{n+1}_i = b_{-1}q^n_{i-1}+b_0q^n_i+b_iq^n_{i+1},\]

It is possible to derive a shortcut to stability, using the von Neumann method, giving these two conditions:

  • \(b_0=\left(b_{-01}+b_1\ge0\right)\),
  • \(b_0=\left(b_{-1}+b_1\right)+4b_{-1}b_1 \ge 0\).

Shortcut to Accuracy

In addition to shortcut for stability, there is shortcut for accuracy, which is more simpler than using LTE.

Linear scheme is \(p\)-th order accurate in space and time in the sense of LTE when satisfying

\[\sum^r_{k=-l}k^\eta b_k=\left(-c\right)^\eta,\quad\eta=0,1,...,p.\]

Various Forms

In addition to general form and conservative form, the numerical schemes can be written as follows:

Incremental Form

\[\begin{aligned} C_{i+1/2}&=C_{i+1/2}\left(q^n_{i-k},q^n_{i-k+1},...,q^n_{i},...,q^n_{i+k}\right) \\ D_{i+1/2}&=D_{i+1/2}\left(q^n_{i-k},q^n_{i-k+1},...,q^n_{i},...,q^n_{i+k}\right) \end{aligned}\]

Using these, a three-point scheme can be written as

\[q_i^{n+1}=q^n_i-C_{i-1/2}\Delta q_{i-1/2} +D_{i+1/2}\Delta q_{i+1/2}\]

where \(\Delta q_{i-1/2} = q^n_i-q^n_{i-1}\) and \(\Delta q_{i+1/2} = q^n_{i+1}-q^n_{i}\)

Viscous Form

\[d_{i+1/2}=d_{i+1/2}\left(q^n_{i-k},q^n_{i-k+1},...,q^n_{i},...,q^n_{i+k}\right)\]

Using this, a three-point scheme can be written as

\[\begin{aligned} q^{n+1}_i = q^n_i&-\frac{1}{2}\frac{\Delta t}{\Delta x}\left[f(q^n_{i+1})-f(q^n_{i-1})\right] \\ &+ \frac{1}{2}(d_{i+1/2}\Delta q_{i+1/2}-d_{i-1/2}\Delta q_{i-1/2}) \end{aligned}\]

where the function \(d_{i+1/2}\) is the coefficient of numerical viscosity, which is different from that of LTE (i.e., \(\alpha_{visc}\)). The second term in equation is called the central part and third term is called viscous part (diffusion part).

For simplicity, the three-point linear scheme can be represented by

\[q^{n+1}_i = b_{-1}q^n_{i-1}+b_0q^n_i+b_1q^n_{i+1}.\]

If this scheme is at least first-order, from the accuracy lemma, the coefficients are given by

\[\begin{aligned} \sum^r_{k=-l}k^0 b_k&=\left(-c\right)^0 \quad\rightarrow\quad b_{-1}+b_0+b_1=1,\\ \sum^r_{k=-l}k^1 b_k&=\left(-c\right)^1 \quad\rightarrow\quad b_{-1}-b_1=c \end{aligned}\]

There are two equations, but three parameters, giving a one-parameter family of solutions. From the first equation,

\[d=b_{-1}+b_1=1-b_0\] \[\begin{aligned} b_{-1}&=\frac{1}{2}(d+c)\\ b_0&=1-d\\ b_1&=\frac{1}{2}(d-c) \end{aligned}\]

Using these, the scheme can be expressed by

\[\begin{aligned} q^{n+1}_i=q^n_i&-\frac{1}{2}c(q^n_{i+1}-q^n_{i-1})\\ &+\frac{1}{2}d(q^n_{i+1}-2q^n_i+q^n_{i-1}) \end{aligned}\]

As mentioned previously, \(d\) is the coefficient of numerical viscosity, second term is called the central part, and third term is the diffusion part.

Using the coefficient of numerical viscosity, the stability condition become much simpler

\[c^2\le d\le 1,\]

and coefficient of numerical viscosity in the sense of truncation error analysis has following relation:

\[\alpha_{visc}=\frac{1}{2}\Delta x \lambda\left(\frac{d-c^2}{c}\right).\]

Particular values of \(d\) give particular schemes as follows:

\[d = \begin{cases} 1 & \rightarrow \text { Lax-Friedrichs } \\ \frac{1}{2}(1+c^2) & \rightarrow \text { FORCE } \\ \vert c\vert & \rightarrow \text { Godunov upwind } \\ c^2 & \rightarrow \text { Lax-Wendroff } \end{cases}\]

Proposition.

Godunov upwind scheme for the linear advection equations is the monotone scheme with the smallest truncation error.

Total Variation

Total Variation

The total variation of a mesh function \(q^n=\{q^n_i\}\), which is set of discretized values, is defined as

\[TV(q^n)=\sum^\infty_{-\infty}\vert q^n_{i+1}-q^n_{i}\vert\]

If the solution \(q\) is given by uniform value, total variation goes to zero, on the other hand, if there is oscillation in solution, total variation increase. Thus, total variation can measure the oscillatory character of the mesh function.

Total Variation Method

Numerical method is TVD (Total Variation Diminishing), if satisfying

\[TV(q^{n+1}) \le TV(q^n)\]

This means that total variation does not increase with time. For example, if the solution is flat at \(t=n\) but has an oscillatory solution at \(t=n+1\), the total variation will increase, therefore this scheme is not total variation method.

Example: Lax-Friedrichs Scheme

\[\begin{aligned} q^{n+1}_i &= \frac{1}{2}(1+c)q^n_{i-1}+\frac{1}{2}(1-c)q^n_{i+1} \\ q^{n+1}_{i+1} &= \frac{1}{2}(1+c)q^n_{i}+\frac{1}{2}(1-c)q^n_{i+2} \end{aligned}\] \[\begin{aligned} q^{n+1}_{i+1}-q^{n+1}_{i} &=\frac{1}{2}(1+c)(q^n_i-q^n_{i-1})+\frac{1}{2}(1-c)(q^n_{i+2}-q^n_{i+1}) \\ \vert q^{n+1}_{i+1}-q^{n+1}_{i}\vert&=\frac{1}{2}(1+c)\vert q^n_i-q^n_{i-1}\vert+\frac{1}{2}(1-c)\vert q^n_{i+2}-q^n_{i+1}\vert \end{aligned}\]

After summation over all integer \(i\),

\[\sum\vert q^{n+1}_{i+1}-q^{n+1}_{i}\vert = \sum\vert q^n_{i+1}-q^{n}_{i} \vert\] \[TV(q^{n+1}) \le TV(q^{n}).\]

Therefore, the Lax-Friedrichs scheme is TVD.

Harten’s Theorem

The set of monotone methods is a subset of the set of TVD methods, and for any scheme in incremental form, sufficient conditions for the scheme to be TVD are as follows:

\[C_{i+\frac{1}{2}}\ge0, \quad D_{i+\frac{1}{2}}\ge0, \quad 0\le C_{i+\frac{1}{2}}+D_{i+\frac{1}{2}}\le1\]

where coefficients are data dependent, i.e., \(C_{i+\frac{1}{2}}=C_{i+\frac{1}{2}}\left(q\right)\) and \(D_{i+\frac{1}{2}}=D_{i+\frac{1}{2}}\left(q\right)\), thus applies to non-linear schemes.

It is noted that monotone methods and TVD methods are only defined for scalar equations, not for systems.