1 Introduction
Lattice Boltzmann schemes are a class of computational methods used to simulate systems of conservation laws under the form of Partial Differential Equations (PDEs). Their basic way of working is the following: instead of taking
PDEs and directly discretize them, a lattice Boltzmann scheme enlarges the size of the problem from to and treats it in a kineticlike fashion. This means that the new variables undergo, at each time step, a local collision phase where different particle distribution functions interact, followed by a latticeconstrained stream phase where no interaction is possible. The advantage of such idiosyncratic approach compared to more traditional numerical methods (e.g. Finite Difference, Finite Volume, Finite Elements, etc.) is that the local nature of the collision phase allows for massive parallelization of the method and the latticeconstrained stream can be computationally implemented as a pointer shift. Although this way of proceeding is highly beneficial from a computational perspective, it yields a deficient structure to construct a clear and rigorous notion of consistency with respect to the target equations, as well as a rigorous theory of stability. Indeed, only formal procedures, either based on the ChapmanEnskog expansion [8] or on the equivalent equations by Dubois [16, 18]are currently available to study the consistency of lattice Boltzmann schemes. As far as stability is concerned, most of the studies rely on the linear stability analysis of the eigenvalues of the system, see
[4, 37].In order to bridge the gap between the lattice Boltzmann methods and the traditional approaches known to numerical analysts, the aim of the present contribution is to show that any lattice Boltzmann scheme can be rewritten as a corresponding multistep Finite Difference scheme on the conserved variables, regardless of the linearity of the equilibria. This is made possible by developing an appropriate formalism based on commutative algebra and therefore yields a proper notion of consistency with respect to the target equations, which is that of Finite Difference schemes (see any standard textbook such as [38]). Furthermore, we confirm that the customary von Neumann analysis used for lattice Boltzmann schemes is equivalent to performing the same analysis on the corresponding Finite Difference scheme and is consequently particularly relevant. The price to pay for passing from an explicit scheme with variables and utilizing information only at the previous timestep to a method with variables is to increase the number of previous timesteps the new solution depends on, yielding a multistep Finite Difference scheme.
In the past, few authors have noticed that for some particular lattice Boltzmann schemes, one has a corresponding (sometimes called “equivalent”) Finite Difference formulation on the conserved variables. Despite this, no general theory has been formulated. For instance: Suga [39] derives by direct computations a threestages Finite Difference scheme from a unidimensional threevelocities scheme,^{1}^{1}1It is customary to call a scheme in a dimensional space using discrete velocities. limiting the computations to a linear framework with one relaxation parameter (SRT). Dellacherie [13] derives a twostages Finite Difference scheme for the lattice Boltzmann scheme. Again, this is limited to one spatial dimension and to a linear framework. A higher level of generality has been reached by the works of Ginzburg and collaborators, see [22] for a recap. They succeeded, using a link formalism, in writing a class of Lattice Boltzmann schemes as Finite Difference schemes [12]
. With their highly constrained link structure to be enforced, the resulting Finite Difference scheme with three stages is valid regardless of the spatial dimension and the choice of discrete velocities. The limitations are that the choice of moments is heavily constrained and only the case of one conserved moment is handled. Moreover, the evolution equation of the moving particles can depend on the distribution of the still particles only
via the conserved moment the equilibria depend upon and the schemes must be tworelaxation time (TRT) models with “magic parameter” equal to onefourth for any link. The difficulty in establishing a general result comes from the coupling between spatial operators and time shifts. We must mention that during the drafting of the present contribution, an interesting work by Fuc̆ik and Straka [21] has been published covering the very same subject and essentially coming to the same conclusion as our paper. Their focus is different than ours since they adopt a purely algorithmic approach rather than a precise algebraic characterization of lattice Boltzmann schemes. We actually provide more insight into the bound on the number of time steps of the corresponding Finite Difference scheme and our formalism, based on polynomials, aims at providing a direct link with the classical tools for the stability analysis and allows to establish a link with the Taylor expansions from [18], as introduced in [3]. In [21], the authors rely on a decomposition of the scheme using an hollow matrix^{2}^{2}2Matrix with zero entries on the diagonal. yielding an equivalent form of the scheme with the diagonal nonequilibrium part, after a finite number of steps of their algorithm. However, to the best of our understanding, the origin of such algorithm is not fully clear. In their work, the spatial shifts of data introduced by the stream phase are taken into account using a rather cumbersome system of indices, whereas we rely on an straightforward algebraic characterization of the stream phase.Our paper is structured as follows: in Section 2
, we introduce – in guise of friendly introduction – the link of our problem with Ordinary Differential Equations (ODEs). The right formalism to make lattice Boltzmann schemes looking very close to a system of ODEs is provided in
Section 3 and allows to prove the main results of the work showcased in Section 4. We devote Section 5 to discuss examples, possible simplifications of the problem and particular cases deserving particular attention. In Section 6, we prove the equivalence of the von Neumann analysis for lattice Boltzmann and Finite Difference schemes. In Section 7, we show how the wellknown tools for Finite Difference schemes can be used to prove convergence theorems for lattice Boltzmann schemes. We corroborate our claim via numerical simulations. We eventually conclude in Section 8.2 The example of Ordinary Differential Equations
Since our way of reducing any lattice Boltzmann scheme to a multistep Finite Difference scheme has been originally inspired by an analogy with systems of ODEs, let us introduce this way of reasoning with the following example. Consider the system of ODEs of size with matrix given by
(1) 
Transforming a single equation of higher order into a system of first order equations like Equation 1 by considering the companion matrix is a current practice, which unsurprisingly makes the problem more handy from the computational standpoint. Though being the analogous of what we aim at doing of lattice Boltzmann schemes, the other way around, passing from a system of first order to a single equation of higher order, seems to be seldom considered. We proceed like in [10]. By iterating, we have that for .^{3}^{3}3We shall consistently use the notation for and . Let be real coefficients, then write . Taking as the coefficients of the characteristic polynomial^{4}^{4}4In the whole work, the indeterminate of any polynomial shall be denoted by . of , by virtue of the CayleyHamilton theorem, we deduce the corresponding equation on the first variable (playing the role of the conserved moment), given by
(2) 
This provides a systematic way of performing the transformation without having to rely on hand computations and substitutions. To give an example, consider
Hence, the corresponding ODE on the first variable is given by .
3 Algebraic form of lattice Boltzmann schemes
Now that the reader is familiar – through a simple example – with the main idea and the final aim of the present contribution, we introduce the general framework of lattice Boltzmann schemes and the right formalism to treat them almost as systems of ODEs.
3.1 Spatial and temporal discretization
We set the problem in any spatial dimension considering the whole space , because we are not interested in studying boundary conditions. The space is discretized by a dimensional lattice of constant step in all direction. The time is uniformly discretized with step . The discrete instants of time shall be indexed by the integer indices so that the corresponding time is . We finally introduce the socalled lattice velocity defined by . Observe that the developing theory is totally discrete and thus fully independent from the scaling between and .
3.2 Discrete velocities and shift operators
The first choice to be made when devising a lattice Boltzmann scheme concerns the discrete velocities with , which are multiples of the lattice velocity, namely for any with . Therefore, particles are stuck to move – at each time step – on the lattice . We denote the distribution density of the particles moving with velocity by for every . The shift operators associated with the discrete velocities are an important element of the following analysis.
Definition 1 (Shift operator).
Let , then the associated shift operator on the lattice , denoted , is defined in the following way. Take be any function defined on the lattice,^{5}^{5}5The function could take values in any ring, see [30]. then the action of is
We also introduce .
The shift yields information sought in the upwind direction with respect to the considered velocity. Let us introduce the natural binary operation between shifts.
Definition 2 (Product).
Let the “product” be the binary operation defined as , for any .
Henceforth, the product is understood whenever no ambiguity is possible. This operation provides an algebraic structure to the shifts, directly inherited from that of .
Proposition 1.
forms an Abelian group.
Moreover, there is only “one movement” for each Cartesian direction which “generates” the shifts. More precisely
(3)  
where is the customary notation for the generating set of a group. We can add one more binary operation, which is noninternal to . This yields the cornerstone of this work, namely the set of Finite Difference operators, finite combinations of weighted shifts operators via a sum. It is defined as follows, see Chapter 3 of [30].
Definition 3 (Finite Difference operators).
The set of Finite Difference operators on the lattice is defined as
the group ring (or group algebra) of over . The sum the product^{6}^{6}6Which interestingly corresponds to the discrete convolution product. of two elements are defined by
Furthermore, the product of with elements of is given by
With the two binary operations, behaves closely to , or as stated by the following result, see [30].
Proposition 2 (Ring of Finite Difference operators).
is a commutative ring.^{7}^{7}7It also an (Hopf) algebra over and can also be viewed as a free module where the scalars belong to and the basis are the elements of the group .
Observe that is not a field: not every element of has multiplicative inverse, take for example the centered approximation of the derivative along : and see for instance the concept of indefinite sum in the calculus of Finite Differences [33, 32]. The elements having inverse are called “units” and divide all the other elements. It can be easily seen that the units are the product of a nonzero real number and a shift in . Indeed for any and . The inverse of a unit shall also be denoted by a bar.
Remark 1.
One can see as the ring of Laurent polynomials of variables over the field , where the indeterminates are , and . For example, for , the identification holds. This automatically implies that is more than a commutative ring, namely a unique factorization domain.
Remark 2.
The reals can be identified with the subring .
3.3 Lattice Boltzmann algorithm: collide and stream
Any lattice Boltzmann scheme consists in an algorithm made up of two phases: a local collision phase performed on each site of the lattice and a stream phase, where particles are exchanged between different sites of the lattice. Let us introduce each of them.
3.3.1 Collision phase
We adopt the point of view of the multiplerelaxationtimes (MRT) schemes, where it is customary to consider the collision written as a diagonal relaxation in the moments basis, see [11]. For this reason, we introduce a change of basis called moment matrix . The entries of can depend on and/or on but cannot be a function of the space and time variables. Gathering the distributions into , the moments are recovered by . We also introduce

the matrix is the relaxation matrix which is a singular with , where is the number of conserved moments:
where the first entries are zero^{8}^{8}8This is not always the case in literature but shall be used consistently in this paper. We put them at the beginning for the sake of presentation. and correspond to the conserved moments, the following are such that for , see [16].

We employ the notation for , where are possibly nonlinear functions of the conserved moments. Since these equilibria are then multiplied by , the first components do not need to be defined.
The collision phase reads, denoting by any postcollision state
(4) 
In the collision phase Equation 4, the entries of can depend on or , but not on space and time. The equilibria are allowed to follow the same dependencies plus those on space and time and can also depend on some “external variable” like in the case of vectorial schemes [23].
3.3.2 Stream phase
The stream phase is diagonal in the space of the distributions. It can be written as
(5) 
where for the first time, the matrices have entries in a commutative ring, see [19] and [6], instead than in the field . The set of square matrices of size with entries belonging to forms a ring under the usual operations between matrices. Even if is commutative from Proposition 2, is not commutative for , as for real matrices and matrices of firstorder differential operators [18].
3.3.3 Monolithic scheme
The stream phase Equation 5 can be rewritten in a nondiagonal form in the space of moments as done by [18, 20] by introducing the matrix and merged with the collision phase Equation 4 to obtain the scheme
(6) 
where and . In the sequel, we shall not indicate the spatial variable for the sake of readability.
We observe that the operators are the eigenvalues of the matrix . However, they are not the eigenvalues of the matrix . Indeed, it is general false that the eigenvalues of belong to the space . It is interesting to interpret the lattice Boltzmann scheme under the form Equation 6 as discretetime linear control system with matrices on a commutative ring [6]. The moments are the state of the system evolving via the matrix , whereas the equilibria are the control via being a feedback observing only a part of the state, namely the conserved moments.
We introduce our example of choice, which shall be used through the whole paper.
Example 1 ( scheme with one conserved moment).
Consider the scheme with one conserved moment [15] by taking , and . We have , and with and
taking and where has been introduced in Equation 3. It can be used to simulate the nonlinear conservation law under the acoustic scaling . The matrices and are
4 Main result of the paper
With a new way of writing any lattice Boltzmann scheme using Definition 3 and thanks to Proposition 2, which provides the ideal setting to generalize the CayleyHamilton theorem, we can proceed like in Section 2 to prove the main result of the paper: any lattice Boltzmann can be viewed as a multistep Finite Difference scheme on the conserved variables.
4.1 Characteristic polynomial and CayleyHamilton theorem
Polynomials with coefficients in and matrices with entries in play a central role in what we are going to develop.
Definition 4 (Characteristic polynomial).
Let be a commutative ring and for some . The characteristic polynomial of , denoted , is given by , where is the determinant and is the identity matrix.
The naive computation of the characteristic polynomial using its definition via the determinant could be computationally expensive, especially when dealing with symbolic computations like in our case. For this reason, we employed the FaddeevLeverrier algorithm [25] which is of polynomial complexity, generally lower than that of the pivot method.
The process is detailed in Algorithm 1 and only uses matrixmatrix multiplications and the computation of the trace, denoted by .
Example 2.
Coming back to Example 1, it is easy to show either by manual computations or by using Algorithm 1 that with
We see that if either or are equal to one, this shall be discussed in Section 5.2. On the other hand if we have .
A central result used in this work is the CayleyHamilton theorem for matrices over a commutative ring, see [6] for the proof, generalizing the same result holding for matrices on a field utilized in Section 2.
Theorem 3 (CayleyHamilton).
Let be a commutative ring and for some . Then is a monic polynomial in the ring in the indeterminate , under the form with . Then^{9}^{9}9Sometimes, we shall indulge to the notation . .
This result states that any square matrix with entries in a commutative ring verifies its characteristic equation.
4.2 Corresponding Finite Difference schemes
The previous Theorem 3 is the key for proving the following results, whose backbone is essentially the same than in Section 2.
4.2.1 One conserved moment
We first analyze the case of one conserved moment, namely , to keep the presentation as simple as possible. We shall eventually deal with once the principles are established.
Proposition 4 (Corresponding Finite Difference scheme for ).
Let , then the lattice Boltzmann scheme Equation 6 corresponds to a multistep explicit Finite Difference scheme on the conserved moment under the form
where are the coefficients of , the characteristic polynomial of .
This result means that the conserved moment satisfies an explicit multistep Finite Difference scheme with at most steps, thus involving discrete time instants, see Figure 1. The maximal size of spatial influence at each past time step can be deduced by looking at Algorithm 1, derived from the Newton’s identities.
It is interesting to observe that also the nonconserved moments satisfy a Finite Difference numerical scheme, see the following proof. However, these schemes would depend on the conserved moment via the equilibria and are therefore not independent from the rest of the system.
Proof.
Let . Then for any , applying Equation 6 recursively we have
We perform a temporal shift in order to fix the first term on the right hand side regardless of the value of . Introduce , therefore
This holds true, in particular, for any . We can then consider the coefficients of the characteristic polynomial of and write
Applying the CayleyHamilton Theorem 3 by virtue of Proposition 2, we know that . Using the monicity of the characteristic polynomial and coming back by setting gives
The last sum can start from . Performing a change of indices in the last double sum yields the result.
(7) 
∎
Example 3.
We come back to Example 1. Using Proposition 4, we have the corresponding Finite Difference scheme given by
(8) 
One can easily check its consistency – under the acoustic scaling – with the target conservation law.
Remark 3.
One could think of allowing and/or to depend on the space and time variables. This would imply to consider weights made up of functions instead of the real numbers in Definition 3. However, would no longer be commutative, because the multiplication by a function does not commute with shifts (not shiftinvariant according to [35]). For example, take and a function , then
for every and for any function . The righthand sides are not equal in general, except if is constant.
4.2.2 Several conserved moments and vectorial schemes
Consider now to deal with multiple conservation laws, namely . We select a conserved moment and we consider the other conserved moments as “slave” variables as the equilibria have been until so far, for , because they imply variables that we eventually want to keep. In particular, we utilize different polynomials for different conserved moments to obtain the Finite Difference schemes. To formalize this concept, for any square matrix , consider for any , corresponding to the matrix where only the rows and columns of indices are conserved and the remaining ones are set to zero. We can also consider the matrix obtained by keeping only the rows and the columns indexed in . A useful corollary of Theorem 3 and of the Laplace formula for the determinant is the following.
Corollary 5.
Let and , then one has that . Moreover, the polynomial annihilates .
This means that the characteristic polynomial of is directly linked to that of the smaller matrix , which is thus faster to compute, and that the latter is an annihilator for the first matrix.
For any conserved moment indexed by we introduce the matrix and . Notice that we have the decomposition . Indeed, we “save” the conserved moments other than the by placing them into , which shall not participate in the computation of the characteristic polynomial. With this notations, we have generated a family of problems from Equation 6 under the form
(9) 
It is useful to stress that the term in Equation 9 does not involve any conserved moment other than the . Conversely, does not involve any function except the conserved moments other than the . Then, the corresponding Finite Difference schemes come under the form stated by the following Proposition.
Proposition 6 (Corresponding Finite Difference scheme for ).
Let , then the lattice Boltzmann scheme Equation 6 rewritten as Equation 9 corresponds to the multistep explicit Finite Difference schemes on the conserved moments under the form
for any where are the coefficients of the characteristic polynomial of .
This Proposition states that for each conserved moment, the corresponding Finite Difference scheme has at most steps, thus involves discrete times. This result encompasses and generalizes Proposition 4. The proof is the same than that of Proposition 4 by taking advantage of Corollary 5. We show in another contribution [3] that the result of Proposition 6 is the right one to bridge between the consistency analysis of Finite Difference schemes and the Taylor expansions on the lattice Boltzmann schemes for proposed by [18].
Example 4 ( for two conservation laws).
Consider the scheme [2] with , and , and
(10) 
thus having . This scheme can be used to simulate the system of conservation laws and under the acoustic scaling . Using Proposition 6 we have
One could remark that the linear part is different from one scheme to the other, since we have used different polynomials for each conserved moment.
4.3 Initialization schemes
In the corresponding Finite Difference schemes in Proposition 4 and Proposition 6, the only remaining freedom is to devise the initialization schemes for the multistep schemes at regime, analogously to Equation 2 in Section 2. This is the counterpart of the freedom of choice on the initial data for the original lattice Boltzmann scheme, which are not necessarily taken at equilibrium, see [27]. By applying Equation 6 to the initial data as many times as needed, one progressively obtains the initialization schemes, as function of the initial datum. It is worthwhile observing that the choice of initial datum does not play any role in the previous procedure and does not influence the stability analysis of Section 6. It only comes into play during the consistency analysis of the numerical method, which is not investigated in this paper, in particular, as far as time boundary layers are concerned, see [40, 34].
5 Examples, simplifications and particular cases
Now that the main results of the paper, namely Proposition 4 and Proposition 6, have been stated and proved, we can analyze and comment some particular cases which deserve a closer look. More examples are available in the Appendices.
Example 5 (ODEs).
To illustrate some basic peculiarities that easily transpose to lattice Boltzmann schemes, we introduce the following matrices extending the discussion of Section 2.
For , we have , corresponding to . However, contrarily to in Section 2, the characteristic polynomial does not correspond to the minimal polynomial . Thus in this case, we could use the latter to obtain Equation 2 having . This phenomenon is studied in Section 5.1. It indicates that we can achieve a more compact corresponding ODE by using the annihilating polynomial of smallest degree on every variable. This does not change the core of the strategy.
For , we obtain , corresponding to . However, by inspecting , one notices that the first two equations do not depend on the last variable . For this reason, we could have considered the matrix obtained from by removing the last row and column. In this case , corresponding to the equation . This kind of situation for lattice Boltzmann schemes is investigated in Section 5.2. It is interesting to observe once more that divides . This shows that an initial inspection of the matrix can yield a reduction of the size of the problem that can be achieved by a simple trimming operation, which eliminates some variable from the problem but treats the remaining ones as usual.
Finally, consider . In this case the characteristic polynomial and the minimal polynomial coincide corresponding to the equation . However, if we take the polynomial such that divides
Comments
There are no comments yet.