Elsevier

Applied Mathematics and Computation

Volume 225, 1 December 2013, Pages 610-621
Applied Mathematics and Computation

The von Neumann analysis and modified equation approach for finite difference schemes

https://doi.org/10.1016/j.amc.2013.09.046Get rights and content

Abstract

The von Neumann (discrete Fourier) analysis and modified equation technique have been proven to be two effective tools in the design and analysis of finite difference schemes for linear and nonlinear problems. The former has merits of simplicity and intuition in practical applications, but only restricted to problems of linear equations with constant coefficients and periodic boundary conditions. The later PDE approach has more extensive potential to nonlinear problems and error analysis despite its kind of relative complexity. The dissipation and dispersion properties can be observed directly from the PDE point of view: Even-order terms supply dissipation and odd-order terms reflect dispersion. In this paper we will show rigorously their full equivalence via the construction of modified equation of two-level finite difference schemes around any wave number only in terms of the amplification factor used in the von Neumann analysis. Such a conclusion fills in the gap between these two approaches in literatures.

Introduction

The von Neumann (discrete Fourier) analysis and the modified equation technique have been proven to be two effective and fundamental tools in the design and analysis of finite difference schemes (even finite element or finite volume methods) for differential equations. The link between them is seemingly well-known and clear, but only in smooth regions of solutions or at low frequencies. The main purpose of this paper is to establish the full connection between them for all wave number modes, for which the modified equation is written in a very compact form just in terms of the amplification factor in the context of von Neumann analysis.

We consider time-dependent partial differential equations with constant coefficients,vt=Lxv,Lx=L(x)=s=1Kαsxs,where Lx represents linear differential operators with constant coefficients αs. A conventionally used two-level difference scheme takes the form,Bun+1=Aun,where the operator B is invertible; or explicitly (1.2) is expressed asp=-N1N2Bpuj+pn+1=p=-M1M2Apuj+pn.Here the discrete value ujn on the grid point (xj,tn) approximates the solution v(xj,tn) of (1.1), xj=jh, tn=nτ, jZ, nN,h represents the spatial mesh size and τ is the length of time step. This scheme has the stencil sizes of the implicit and explicit parts N=N1+N2+1 and M=M1+M2+1, respectively. The coefficients Ap and Bp are constant, satisfying a basic consistency condition,p=-N1N2Bp=p=-M1M2Ap.We are particularly concerned with explicit schemes, for which N1=N2=0 and so N=1. In fact, in our analysis, the implicit scheme (1.3) can be converted into an explicit one by a linear transformation.

To analyze and use the discrete algebraic equation (1.3) in the approximation of the solution of the continuous equation (1.1), four most basic and fundamental issues are the consistency, stability, convergence and wave propagation. The consistency describes whether (1.3) approximates (1.1). Once (1.3) is consistent with (1.1), the Lax-equivalence theorem states that the solution of (1.3) converges to the solution of (1.1) if and only if the scheme (1.3) is stable [10]. The wave propagation can be understood by using a so-called dispersion relation. All these issues can be addressed with the von Neumann analysis [2] and the modified equation technique [24]. The von Neumann analysis is relatively simple and intuitive, and the unique amplification factor implies the consistency, determines the stability of (1.3) and displays properties of wave propagations, despite the fact that it mostly applies to equations with constant coefficients and periodic boundary condition [16]. The modified equation can be interpreted as the actual partial differential equation that the scheme (1.3) solves [24]. This approach is originated in [4] and can be generalized for linear problems with variable coefficients or even nonlinear problems [8], the dissipation and dispersion effects of (1.3) can be understood from the partial differential equation point of view: Even order terms determines the dissipation and odd order terms represents the dispersion. This dissipative and dispersive information derived from the modified equation approach is very heuristic, unfortunately just valid for solutions in smooth regions or at low frequency modes [24]. Therefore the connection with the von Neumann analysis is only restricted there.

However, in many applications, for example, in the simulation of turbulent flows, aeroacoustic computations, or in the analysis of numerical boundary conditions, the understanding of numerical (e.g., oscillatory) phenomena and many others [3], [18], one is more concerned with solution behaviors for a large range of wave numbers, particularly, middle- or large-wave-number, which arises a question how to derive modified equations for any wave number Fourier mode, or establish the full connection between these two approaches. This kind of missions has started earlier though just sporadic. For oscillatory (large-wave-number) parts, a formal analysis is presented in [16] for some specific schemes (e.g., the box scheme) and in [12] for the stable three-point generalized Lax-Friedrichs schemes. The oscillatory behavior of solutions can be illustrated rather clearly via the associated modified equation.

As such, the modified equation approach is worth reinspecting. It was pointed out in [7]: “Because of the lack of any theoretical foundation, this use has been accompanied by constant difficulties and results derived from modified equations have sometimes been regarded with apprehension. As a result a situation has arisen where authors either disregard entirely the technique or have an unjustified faith in its scope.” Nevertheless, this approach is very heuristic in the design of numerical schemes and applicable to nonlinear problems [5], [17], [20], [6], [22], and therefore we hope to clarify the connection with the classical discrete Fourier analysis for the full range of wave number modes.

The full connection can be intuitively built as follows. Consider a Fourier mode ujn=λneijkh, i2=-1,k is real. The frequency ω and the amplification factor λ are linked with λ=eiωτ, and thereby the dispersion relation ω=ω(k;h,τ) for (1.3). Then the (discrete) frequency ω is a function of k, but also of h and τ. We will simply write ω as ω=ω(k), holding h,τ fixed later on. Then this Fourier mode can be expressed asujn=e-ωInτeik[jh-vp(k)nτ],where vp(k)=-ωR(k)/k is the phase speed, ωR(k) and ωI(k) are the real and imaginary parts of ω(k), respectively, i.e., ω(k)=ωR(k)+iωI(k). Hence this mode propagates with the phase speed vp(k) and is dissipated with the rate -ωI(k) at each step if ωI>0, which is just the von Neumann criterion for stability1|λ(k)|1.To understand how the overall shape of the wave amplitudes around a fixed wave number k0, let us think of k as the wave number around k0,k=k0+k̃ and thereby ω=ω0+ω̃, ω0=ω(k0),k̃ and ω̃ are small. Then with the smoothness assumption of ω̃ in terms of k̃, one obtainsω(k)=ω(k0)+ω̃(k0+k̃)=ω0+s=1ω̃(s)(k0)s!k̃s,ω̃(k0)=0.In this paper we always use the notation ω(s)(k) to denote the s-th derivative of ω in terms of k and similarly for others. Then this Fourier mode propagates, see [21], in the form,ujne-ωI(k0)nτ·eik0[jh-vp(k0)nτ]·eik̃[jh-G(k0)nτ]by ignoring higher order terms. Here G(k0)=-dω̃(k0)/dk is the group velocity at k0. Hence the principal part of this mode propagates with the phase speed vp(k0) and is damped with the rate -ωI(k0) at each time step. The dispersion of (1.3) and other dissipation effects are reflected through ω̃(k0+k̃).

We can depict the propagation of the overall shape of this wave amplitudes (wave packet) in terms of modified equations. Setujn=αjnũjn,αjn=ei(nω0τ+jk0h),ũjn=ei(nω̃τ+jk̃h),i2=-1.With the change of k̃,ũjn is the wave packet around (k0,ω0). By comparing this with (1.7), and substituting ik̃ by x and iω̃ by t, the modified equation for ũjn is expressed in a compact formtũ=1τs=1i-ss![lnλ](s)(k0)xsũ,[lnλ](s)(k0)=iω̃(s)(k0)τ.where λ(s)(k0)=dsλ(k0)dks. In terms of the group velocity G(k0), (1.10) is rewritten astũ+G(k0)xũ=-s=2i-(s-1)s!G(s-1)(k0)xsũ,which implies that the wave packet of ũ propagate with the group velocity G(k0), and the behavior is totally determined in terms of this unique quantity. This modified equation approach provides a PDE interpretation of the wave packet, through which we clarify its full connection with the von Neumann analysis: (i) Take the von Neumann analysis for (1.3) to obtain the amplification factor λ(k0) for a fixed wave number k0. Then we derive the modified equation (1.10) or (1.11) for the wave packet of ũ. (ii) On the other hand, as long as all coefficients of the modified equation at a fixed wave number k0 are known, all derivatives of the amplification factor in the von Neumann analysis are obtained and thus the amplification factor λ(k0) is known. Thus the von Neumann analysis and the modified equation technique are unified in a common framework and we can conveniently adopt either the von Neumann analysis or the modified equation to analyze the dissipation (stability) and dispersion (oscillations) of the scheme (1.3). We point out particularly that as (k0,w0)=(0,0), this modified equation (1.10) becomes that in [24] in a compact form,tu=1τs=1i-ss![lnλ](s)(0)xsu.This implies that the frequency defined by the dispersion relation for (1.3) be just the symbol of the modified equation (1.12).

We remark that although the above analysis is made only for (1.1) with constant coefficients, it works for the local stability or instability of the scheme (1.3) with variable coefficients or nonlinear problems. Indeed, we have developed a heuristic modified equation approach for nonlinear equations [13]. In a forthcoming contribution [14], we will investigate variable coefficient cases.

We organize this paper as follows. In Section 2 we will derive the modified equation rigorously. In Section 3 we present a throughout explanation about the full connection. In Section 4 we provide some simple discussions.

Section snippets

Modified equation for a wave packet around a fixed wave number

This section provides the derivation of a modified equation for a wave packet around a fixed wave number k0 and the corresponding frequency ω0=ω(k0). Hence we make an ansatz,ujn=αjnũjn,αjn=ei(nω0τ+jk0h),ω0=ω(k0).This formulation somewhat corresponds to the setting in the Fourier analysisujn=einωτeikjh,ũjn=einω̃τeik̃jh,ω=ω0+ω̃,k=k0+k̃,and extends the analysis in [12], [16] for highest frequency modes. As (k0,ω(k0))=(0,0), αjn1 and ujn=ũjn.

With the assumption that k̃1, it is shown in [21, p.

The Full connection between two approaches

In Section 2 we derived the modified equation (2.4) for the scheme (1.3). In this section we present a throughout explanation about its full connection with the von Neumann analysis in terms of consistency, stability and wave propagation, although it is quite clear from the previous discussions. We assume again that (1.1) takes an explicit formLxv=s=1Kαsxsv,where αs are constant, s=1,2,,K. The symbol of the operator Lx is a polynomial, denoted bySL(k)=L(ik).Then the elementary solution u(x,t)

Some extensions and discussion

The von Neumann analysis has the advantage of simplicity, intuition and easy implementation and it can supply at least local information about properties of dissipation and dispersion, even though there are the severe restriction over constant coefficients of PDE and periodic boundary conditions, and even counterexamples about the frozen coefficient method which leads to the instability of variable coefficients problems [15]. In contrast, the application of the modified equation approach may be

Acknowledgement

Jiequan Li is supported by NSFC (91130021, 11031001), BNU Innovation Funding (2012LZD08) and a special funding from Institute of Applied Physics and Computational Mathematics, Beijing. Both authors thank the reviewer for his helpful suggestions.

References (24)

  • D.F. Griffiths et al.

    On the scope of the method of modified equations

    SIAM J. Sci. Statist. Comput.

    (1986)
  • A. Harten et al.

    On finite-difference approximations and entropy conditions for shocks

    Commun. Pure Appl. Math.

    (1976)
  • Cited by (16)

    • On some time marching schemes for the stabilized finite element approximation of the mixed wave equation

      2015, Computer Methods in Applied Mechanics and Engineering
      Citation Excerpt :

      This sometimes cannot be achieved, thus one just aims to have low-dispersion low-dissipation schemes [1,2]. Dispersion and/or dissipation of numerical schemes can be evaluated using Fourier techniques (see [3–5]), energy methods (see [6]), or modified equation analysis (see [4,7]). Fourier analysis can be carried out for semi-discretizations or full discretizations.

    • Taylor series expansion using matrices: An implementation in MATLAB®

      2015, Computers and Fluids
      Citation Excerpt :

      Although the procedure for determining the modified equation is extremely laborious, the final results contains highly valued information regarding the nature of the truncation error for numerical schemes involving two or less time levels. Li and Yang [17] showed that the von Neumann analysis and the modified equation approach are fully equivalent, yet the latter is able to handle non-linear problems. The modified equation approach, originally developed by Warming and Hyett [18], starts by substituting the TSE of each term (to desired order) in the FDE, then a table is constructed to eliminate the higher time derivatives by means of a linear combination of derivatives of the same expanded FDE.

    View all citing articles on Scopus
    View full text