A stability property for a mono-dimensional three velocities scheme with relative velocity

by   François Dubois, et al.
Paris-Sud University

In this contribution, we study a stability notion for a fundamental linear one-dimensional lattice Boltzmann scheme, this notion being related to the maximum principle. We seek to characterize the parameters of the scheme that guarantee the preservation of the non-negativity of the particle distribution functions. In the context of the relative velocity schemes, we derive necessary and sufficient conditions for the non-negativity preserving property. These conditions are then expressed in a simple way when the relative velocity is reduced to zero. For the general case, we propose some simple necessary conditions on the relaxation parameters and we put in evidence numerically the non-negativity preserving regions. Numerical experiments show finally that no oscillations occur for the propagation of a non-smooth profile if the non-negativity preserving property is satisfied.



There are no comments yet.


page 1

page 2

page 3

page 4


Quasi-linear analysis of dispersion relation preservation for nonlinear schemes

In numerical simulations of complex flows with discontinuities, it is ne...

Stabilized integrating factor Runge-Kutta method and unconditional preservation of maximum bound principle

Maximum bound principle (MBP) is an important property for a large class...

Energy Stability of Explicit Runge-Kutta Methods for Non-autonomous or Nonlinear Problems

Many important initial value problems have the property that energy is n...

Maximum bound principle preserving integrating factor Runge-Kutta methods for semilinear parabolic equations

A large class of semilinear parabolic equations satisfy the maximum boun...

Energy plus maximum bound preserving Runge-Kutta methods for the Allen-Cahn equation

It is difficult to design high order numerical schemes which could prese...

Type-Preserving Matrices and Security of Block Ciphers

We provide a new property, called Non-Type-Preserving, for a mixing laye...

An upwind DG scheme preserving the maximum principle for the convective Cahn-Hilliard model

The design of numerical approximations of the Cahn-Hilliard model preser...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

Studying stability of lattice Boltzmann schemes is a non-trivial problem. Classically for this purpose, the scheme is linearized around a constant state and a Fourier analysis is performed. We refer to the work of Lallemand and Luo [9] for the D2Q9 scheme applied to hydrodynamics. Note also the work of Ginzburg et al. [8] extending the Fourier analysis to a wide variety of different two and three dimensional lattice Boltzmann schemes. Instabilities and their interpretattion in terms of bulk viscosity has been proposed by Dellar [2]. But no mathematical analysis has been performed.

A new way of improving stability is proposed by Geier [7], who proposed a new generalized lattice Boltzmann scheme with the approach of relative velocities and utilized it for hydrodynamics applications [4]. The tentative of analysis of this method for a two-dimensional scalar linear equation has been also proposed [5].

Even if it is a difficult task, it is well known that Fourier analysis is not the best method for analyzing nonlinear hyperbolic equations. Total variation diminishing schemes, developed for suppressing oscillations in higher order CFD algorithms, provide an alternative nonlinear stability analysis tool for analysing the schemes for nonlinear wave propagation. The convergence of such schemes is well established [10]. The underlying stability notion concerns the maximum principle. This notion can be extended to nonlinear cases and a first attempt has been proposed in [1] for lattice Boltzmann schemes for the D1Q2 scheme for scalar nonlinear hyperbolic equations.

In this contribution, we propose to investigate the stability in the maximum sense, of a linear mono-dimensional lattice Boltzmann scheme with three velocities. More precisely, we look to a positivity constraint for a particle distribution function in the context of a relative velocities.

In Section 2, we describe the scheme and the underlying advection model. More precisely, the local relaxation step is written as a linear operator on the particle distribution functions. If all the coefficients of the underlying matrix are nonnegative, the non-negativity of the distribution is maintained during this step. Because the transport step is just a change of locus, the non-negativity is maintained for the whole time step of the scheme. The question is then to find appropriate conditions to handle this property. In Section 3, a necessary and sufficient condition is derived on the parameters in order to ensure that the scheme has the stability property. In Section 4, we completely describe the classical case where the relative velocity is reduced to zero. In Section 5, the general case is presented. With an analytical study for necessary conditions and a numerical one for a complete description of the stability zones. In Section 6, numerical experiments show the correlation of the positivity constraint for a particle distribution and the presence of oscillations for discontinuous profiles.

2. Description of the framework

2.1. Description of the scheme

In this contribution, we investigate a mono-dimensional three velocities linear lattice Boltzmann scheme with relative velocity [4]. Denoting the spatial step, the time step, and the lattive velocity, this scheme can be described in a generalized D. d’Humière’s framework [3]:

  1. [(1)]

  2. the 3 velocities , , and ;

  3. the 3 associated distributions , , and ;

  4. the 3 moments

    , , and given by

    where is a given scalar representing the relative velocity;

  5. the equilibrium value of the 3 moments

    where and are given scalars (without loose of generality, we assume that );

  6. the 2 relaxation parameters and such that the relaxation phase reads

In this formalism, the moments are defined as polynomial functions of the discrete velocities and the discrete distribution functions. Indeed, introducting , , and , the three moments read

The equilibrium values are chosen such that the equilibrium distributions do not depend on the relative velocity . Indeed, we have:

Note that this scheme can be used (see e.g. [6]) to simulate a scalar transport equation with constant velocity given by

Denoting , , the discret spatial points and , , the discret time points, one time step of the scheme reads

2.2. Matrix notation for the relaxation step

In this contribution, we are concerned with the non-negativity of the particle distribution functions. This property can be viewed as refering to a weak maximum principle for linear schemes. Indeed, it is always possible, by adding a constant, to assume that all the particle distribution functions are initialy nonnegative. If the scheme ensures that this property of non-negativity remains as time marches, each particle distribution function is then bounded, as their total sum is conserved. Moreover, as the transport step consists simply in exchanging the position of the particle distribution functions, we focus on the relaxation step.

We use a matrix notation for the relaxation step as it can be read as a multiplication by a matrix. As this step is local in space, we omit the dependancy on time and space. We define the vector of the distribution functions

One relaxation step then reads

where the matrix is defined by


The coefficients of the matrix are obtained by the relations , and those of the matrix by the change of basis formula: is the coefficient of the -element in the definition of according to

The matrix is then the change of basis that transforms the vector into the vector :

The matrix can then be viewed as the change of basis matrix from the classical moments without relative velocity toward the moments with relative velocity.

2.3. Remark on the choice of the moments

Note that the last moment that is chosen in this contribution is not the energy but a moment that is orthogonal to the two first ones, and . In this section, we show that all the results of the contribution would be identical by choosing the last moment as the energy: the relaxation matrix would still be the same.

We consider two schemes with two different choices of polynomials: the moments of the first scheme are defined by while the moments of the second scheme by . The first moment is the same in both schemes to be able to simulate the same transport equation. We then have . We define the change of basis matrix associated to the tranformation into :

The first line of is then .

We assume that the equilibrium values of the distribution functions are the same, that is , and that the relaxation parameters are the same, that is . Then, we have for all iff


First, we immediately obtain the following relations by identifying the coefficients of and of :

We deduce that

Moreover, as the second and the third columns of are zero (the equilibrium values depend only on the first moment ), and as the first line of is then , we have

We have


As the matrices , , and are invertible, the condition is equivalent to . Denoting

a straightforward calculation yields

Then the property for all values of and is equivalent to , that ends the proof. ∎

2.4. Positivity of the matrix

The velocity being fixed, we propose to give a full description of the sets

Indeed, the non-negativity of the matrix imposes that all the distributions , , remain non-negative if they are so at the initial time. These sets are first described by a set of nine inequalities that can be joined into just one. Numerical illustrations are then given to visualize it in the characteristic cases including SRT, MRT, and relative velocity scheme.

3. Positivity of the iterative matrix

The nine inequalities obtained from the matrix can be combined neatly into one formula.

The inequalities are


We prove now that the previous nine inequalities can be written in a much more lucid way.

We introduce the reduced parameters and according to


Then the nine previous inequalities displayed in (1) are equivalent to


Consider first the two inequalities associated with and :

They can be synthetized in the following form

and we can write this relation as


Write now the inequalities (1) associated with and :

In other words, . Then


We now focus on the inequalities (1) associated with and :

We have . In consequence,


Considering the inequalities (1) with and , we have

and . In consequence,


The last inequality can be written as and this inequality is equivalent to . In other terms,


The inequalities (4), (5), and (6) establish a triple majoration of whereas the inequalities (7) and (8) show a double minoration of the same quantity. The proof is completed. ∎

4. The particular case

In this section, we suppose that the relative velocity is reduced to zero. Then the necessary and sufficient conditions (3) for the stability can be written as


To fix the ideas, we suppose that the advection velocity is non-negative:


The case follows directly. When , the reduced stability conditions


are equivalent to the following conditions for the relaxation parameters


joined with a natural Courant type condition for explicit schemes on the advection velocity


Of course, the conditions (9) have still to be imposed for the equilibrium parameter when the pair is given. In particular,



Figure 1. Necessary and sufficient stability regions described by the inequalities (11) for a null relative velocity . Illustration proposed for .

We first observe that . Then . Secondly, we have and because both and are positive, we have . We have also and (13) is established.

Moreover, and implies . Due to the positivity of the parameter , we deduce from (9) a new set of inequalities: . Then . Moreover, because .

Conversely, if the relations (12) and (13) are satisfied, we have , and . Then . Moreover, , because and . Thus .

Finally the inequalities (11) are established and the proposition is proven. ∎

5. The general case

We analyse the general case of non-zero in this section, with suitable illustrations of the stability region for various ranges of the parameters.

5.1. Necessary conditions for stability

In this subsection, we prove the following proposition.

We suppose that the necessary and sufficient stability conditions (3) are satisfied:

with the notations (2):

We suppose also that the advection velicity is positive. Then if the scheme is stable, the point satisfies the following inequalities


With these necessary stability conditions, the parameter has been eliminated. We have also, necessarily




We start from the inequalities (3). Then we have and . With the same type of inequalities, and .

If we remark that , we deduce that and the inequality (18) is proven.

We have the triangular inequality . Then from the genaral stability conditions (3), we deduce and . In a similar way, , from (3) and finally . We put together the two inequalities and we have .

In consequence, we have and the third point is proven. Since we made the choice of then . The first inequality of the two first points are established. From we have and the relation (17) is true. Moreover, and the last inequality of (16) is true.

Consider now the inequalities . We deduce and . Due to the positvity of the absolute values, we have also and . A part of the fifth inequality of (16) is proven.

Finally, due to the triangular inequality, . We add this inequality with the two following ones: and . Then and . The fifth inequality of (16) is completely established. Because , we have also and the first inequality of (16) is also established. ∎

We illustrate the zones of necessary stability in the Figure 2 for five particular velocities: .

Figure 2. Necessary stability regions given by inequalities (16) for (top-left), (top-right), (middle-left), (middle-right), (bottom-left), (bottom-right).

5.2. Numerical study of necessary and sufficient conditions for stability

We now illustrate the necessary and sufficient stability regions for , , , , , and for various ranges of in the Figure 3.

Figure 3. Nunerical study of necessary and sufficient stability regions

For each value of the velocity , the necessary and sufficient stability regions given by the relations (3) are displayed for several values of the relative velocity: . In dotted line, the necessary stability region of the Proposition 5.1 is added for remind. The particular value is enhanced by filling the region in gray.

Some analysis can be drawn from these figures:

  • the stability region changes with the relative velocity;

  • the maximal value of the first relaxation parameter is obtained (not only) for ;

  • the stability region is not clearly more interessant (larger or including greater value of ) for ;

  • the point is always in the stability region.

To conclude this section, this notion of stability allows a large set of values for the relaxation parameters. If the scheme is used to simulate the hyperbolic advection equation without diffusive term, we can try to minimize the numerical diffusion while maintaining this stability property. This task is complicated and out of the scope of the paper as the second-order term reads as a non-linear formula that links all the parameters , , and .

6. Numerical illustrations

In this section, we illustrate the stability property with numerical simulations involving the D1Q3 model with and without relative velocity used to simulate the linear advection equation. The stability is demonstrated for the parameters chosen according to the analysis presented in the previous sections. Oscillations are seen whenever the parameters go beyond the stability limits presented, as highlighted in the results.

The parameters chosen for the simulations are the following:

left (stable) 0.25 0.0 1.6 1.3 0.3076923076923076
0.25 0.25 1.6 1.3 -0.17548076923076938
right (unstable) 0.25 0.0 1.9 1.4 0.14285714285714302
0.25 0.25 1.9 1.4 -0.10491071428571441

For the left figures, the parameters are chosen in order to satisfy the stability property. We observe numerically also a maximum principle, even if it is not formally the stability notion that we investigate. For the right figures, the parameters are chosen in order not to break the stability property.

Figure 4. Smooth profile with a continuous derivative. The parameters for the D1Q3 are tuned in order to have (left) or not (right) the non-negativity property
Figure 5. Continuous profile. The parameters for the D1Q3 are tuned in order to have (left) or not (right) the non-negativity property
Figure 6. Discontinuous profile. The parameters for the D1Q3 are tuned in order to have (left) or not (right) the non-negativity property

If the profile is smooth, we have observed that no numerical oscillations occur even if the non-negativity property of the matrix is not satisfied. If the profile is just continuous, small negative values of the macroscopic quantity are observed when our non-negativity property is not satisfied. Last but not least, classical oscillations are visible for discontinuous profiles if our non-negativity property is not satisfied. These oscillations are eliminated when the non-negativity property of the matrix is realized.

7. Conclusion

In this contribution, we have investigated a stability property for a classical mono-dimensional linear three velocities lattice Boltzmann scheme with relative velocity. This property ensures that the particle distribution functions remain nonnegative if they are at initial time. We then give a necessary and sufficient condition to describe the stability region. The case without relative velocity is completely describe and simplest necessary conditions are given for the general case. We finally propose some numerical simulations that illustrate the stability property: even if the stability notion that we investigate is not exactly a cosntraint of convexity, a numerical maximum principle is observed if the parameters are in the stability region whereas numerical oscillations appear (in particular for non smooth profils) if the parameters are not there.

Moreover, relative velocities modify the stability array in a nontrivial manner. For instance, intuition might have suggested that the stability region for the relative velocity equal to the advection velocity contains all the others but it is realy not the case. For a given advection velocity, relative velocities cannot be used to increase the value of the first relaxation parameter, the one involved in the numerical diffusion.

The non-negativity of the relaxation matrix could be extended to nonliear schemes. The theoretical study will be much more technical and has not been performed. Nevertheless, numerical experiments for the Burgers equation show that the behavior of the D1Q3 scheme, and in particular the possibility or not of oscillations, is analogous to the linear case.


  • [1] F. Caetano, F. Dubois, and B. Graille. A result of convergence for a mono-dimensional two-velocities lattice Boltzmann scheme. arXiv:1905.12393 [math], May 2019. arXiv: 1905.12393.
  • [2] P. J. Dellar. Nonhydrodynamic modes and a priori construction of shallow water lattice Boltzmann equations. Physical Review E, 65(3):036309, February 2002.
  • [3] D. d’Humières. Generalized Lattice-Boltzmann equations. In Rarefied Gas Dynamics: Theory and simulation, volume 159, pages 450–458. AIAA Progress in astronomics and aeronautics, 1992.
  • [4] F. Dubois, T. Février, and B. Graille. Lattice Boltzmann schemes with relative velocities. Communications in Computational Physics, 17(4):1088–1112, 2015.
  • [5] F. Dubois, T. Février, and B. Graille. On the stability of a relative velocity lattice Boltzmann scheme for compressible Navier–Stokes equations. Comptes Rendus Mécanique, 343(10-11):599–610, 2015.
  • [6] F. Dubois and P. Lallemand. Towards higher order lattice Boltzmann schemes. Journal of Statistical Mechanics: Theory and Experiment, 2009:P06006, 2009.
  • [7] M. Geier, A. Greiner, and J. C. Korvink. Cascaded digital lattice Boltzmann automata for high reynolds number flow. Physical Review E, 73:066705, 2006.
  • [8] I. Ginzburg.

    Truncation Errors, Exact And Heuristic Stability Analysis Of Two-Relaxation-Times Lattice Boltzmann Schemes For Anisotropic Advection-Diffusion Equation.

    Communications in Computational Physics, 11(5):1439–1502, May 2012.
  • [9] P. Lallemand and L.-S. Luo. Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability. Physical Review E, 61:6546–6562, 2000.
  • [10] S. Osher. Convergence of generalized muscl schemes. SIAM Journal on Numerical Analysis, 22(5):947–961, 1985.