Numerical analysis of a Lax-Wendroff type lattice Boltzmann scheme for the advection equation

Thomas Bellotti; Filipa Caetano; François Dubois; Benjamin Graille

Source-linked AI summary

Numerical analysis of a Lax-Wendroff type lattice Boltzmann scheme for the advection equation

Thomas Bellotti, Filipa Caetano, François Dubois, Benjamin Graille

arXiv:2610.32810v1Oct 23, 2026math.NAcs.NA

TL;DR

Stability and convergence are central questions for lattice Boltzmann schemes for evolutionary equations. This paper analyzes a three-velocity scheme for linear advection, deriving explicit stability conditions and proving convergence under specified smooth-data and initialization conditions. The analysis shows that the relaxation parameters govern stability without changing the scheme’s accuracy order.

Problem
Stability is an essential question for numerical schemes for evolutionary equations, and the paper also seeks to establish convergence of the numerical solution to the target equation.
Source: §1 · p. 2
Method
The paper combines von Neumann analysis with eigenvalue-structure analysis to characterize stability and uses finite-difference results to study convergence.
Source: §8 · pp. 36–37
Results
The relaxation parameters do not affect the scheme’s order of accuracy but determine its stability properties, and convergence is established for smooth initial data with non-conserved moments initialized at equilibrium.
Source: §8 · pp. 36–37
Takeaways & Limitations
The scheme recovers Lax-Wendroff at unit relaxation parameters, connecting the lattice Boltzmann formulation to a classical finite-difference scheme.
Source: §8 · p. 36
The convergence analysis considers relaxation parameters in (0, 2)^2 and requires the corresponding finite-difference scheme to be L2-stable.
Source: §6 · p. 32

Abstract

from arXiv · show

A linear monodimensional lattice Boltzmann scheme with three velocities tackling the advection equation is investigated. We prove that it is second-order accurate and provide necessary and sufficient stability conditions through an explicit formula involving the two relaxation parameters. Relying on existing results for finite difference schemes, we prove convergence for smooth initial data with initialization of the non-conserved moments at equilibrium. Several numerical simulations illustrate these theoretical findings.

1. Introduction

The paper studies a three-velocity lattice Boltzmann scheme for linear advection, with second-order accuracy and a Lax-Wendroff special case. It gives explicit stability conditions and examines convergence under stated initialization and smooth-data conditions.

1. Introduction

The scheme achieves second-order accuracy and coincides with the Lax-Wendroff method when both relaxation parameters equal one.The scheme retains its centered character for general parameter choices.Source: §1 · pp. 1–2
The two relaxation parameters determine stability while leaving accuracy uniform in their values.The two linear equilibria are fixed to obtain second-order accuracy, leaving the relaxation parameters free.Source: §1 · pp. 1–2
Explicit necessary and sufficient stability conditions are derived in terms of the Courant number and both relaxation parameters.The analysis supplements von Neumann analysis with eigenvalue-multiplicity considerations.Source: §1 · p. 2
Convergence of the conserved moment is established for smooth initial data when non-conserved moments are initialized at equilibrium.The proof reuses known finite-difference results.Source: §1 · p. 2
Numerical experiments illustrate wave propagation, convergence for smooth and non-smooth data, and empirical second-order convergence for an inviscid Burgers adaptation.The Burgers experiment concerns an adaptation of the scheme.Source: §1 · p. 3

2. Description of the scheme

The scheme uses three velocities on a one-dimensional grid to approximate a linear scalar advection equation. Its moments and relaxation parameters define the lattice Boltzmann update, which is also represented in Fourier space by an amplification matrix.

2. Description of the scheme

The target is the linear scalar advection equation with constant velocity, discretized on uniform space and time grids.The grid ratio is defined as λ := ∆x/∆t and assumed fixed.Source: §2 · p. 3
The three discrete velocities are 0, λ, and −λ, with particle distributions f0, f+, and f− and moments u, v, and w.The moment definitions include u = f0 + f+ + f− and v = λf+ − λf−.Source: §2 · p. 3
The scheme introduces two relaxation parameters for the non-conserved moments, while the equilibrium parameters are selected later for consistency.The centered relaxation parameters are tv = sv − 1 and tw = sw − 1.Source: §2 · p. 3
A spatial Fourier transform expresses the scheme through an amplification matrix G, which can be rewritten using centered relaxation parameters.The matrix is independent of the second-order moment choice under the stated conditions on that moment and unchanged particle equilibria.Source: §2 · p. 4

3. Consistency

The consistency analysis connects the lattice Boltzmann dynamics to an equivalent finite-difference scheme and derives modified equations. Choosing the equilibrium parameters appropriately gives second-order accuracy, subject to relaxation-parameter conditions.

3. Consistency

An exact correspondence with an equivalent multistep finite-difference scheme provides access to classical finite-difference analysis.The correspondence supports systematic derivation of associated macroscopic equations.Source: §3 · p. 4
Taking λc = V gives first-order accuracy, and additionally choosing d = c² yields second-order accuracy regardless of sv.The modified equation is used to analyze consistency.Source: §3 · p. 5
The product of centered relaxation parameters is called the magic parameter; taking it equal to 1/12 is associated with third-order accuracy.The paper notes that this accuracy result is known in the lattice Boltzmann community.Source: §3 · p. 5
The consistency analysis assumes nonzero sv and sw; when sv = 0 the scheme is inconsistent, while sw = 0 generally leaves only first-order accuracy.For sw = 0, second-order accuracy is recovered when sv = 2.Source: §3 · p. 5
The consistency analysis motivates linking the initial moments to the Cauchy problem’s initial datum so the conserved moment can approximate the PDE solution.The paper describes this as the intended connection for subsequent analysis.Source: §3 · p. 5

4. Stability

The section develops an ℓ2-stability analysis for the three-velocity scheme and derives a necessary-and-sufficient condition in its relaxation parameters. It uses polynomial tests to handle stability and eigenvalue multiplicities, then states the resulting parameter conditions and special cases.

4.1. Reminder on complex polynomials.

The stability analysis uses Schur, von Neumann, and simple von Neumann polynomial tests, whose recursive criteria reduce polynomial degree.The section defines Schur polynomials by roots strictly inside the unit circle and von Neumann polynomials by roots on or inside it.Source: §4.1 · pp. 5–6

4.2. Definition of the stability.

Stability requires the minimal polynomial of the amplification matrix to be simple von Neumann; that condition is necessary, and the characteristic polynomial supplies a sufficient test.When the characteristic polynomial is von Neumann but not simple, the case requires further investigation.Source: §4.2 · pp. 7–8
The characteristic polynomial also recasts the conserved variable’s dynamics as a multistep finite-difference scheme, but the two schemes’ stability is not generally equivalent.The paper notes a counterexample in which the lattice Boltzmann scheme is stable while the resulting finite-difference scheme is not.Source: §4.2 · p. 8

4.3. Main result: the stability condition.

For stable schemes other than the fully degenerate relaxation pair (tv, tw) = (−1, −1), the centered parameters lie in [−1, 1] and |c| ≤ 1.The fully degenerate pair is stable for any c.Source: §4.3 · p. 9
The D1Q3 scheme is stable if and only if its relaxation parameters satisfy the stated conditions, which depend on the velocity c.For 0 < |c| < 1, stability requires 0 ≤ sv < 2 and 0 ≤ sw ≤ s⋆w(sv); the theorem gives separate conditions for c = 0, |c| = 1, and |c| > 1.Source: §4.3 · p. 10
For 0 < |c| < 1, the stable region in (sv, sw) does not depend on c.Figure 1 marks stable points in green and unstable points in purple.Source: §4.3 · p. 10, Figure 1 · p. 11
With a single relaxation parameter, stability requires s ≳ 1.1010 when 0 < |c| < 1, and 0 ≤ s < 2 when |c| is 0 or 1.The threshold for 0 < |c| < 1 is defined through a fixed-point condition.Source: §4.3 · p. 10
For the scheme introduced in with sw = 2 − sv, stability requires 0 ≤ sv < 2 when 0 < |c| < 1, and 0 < sv < 2 when c ∈ {0, 1}.Source: §4.3 · p. 11

5. Proof of the stability condition

This section proves the stability criterion by analyzing interior relaxation parameters and classifying stability on parameter edges and vertices. It also identifies an admissible range for the second centered relaxation parameter and considers stability of specific parameter pairs across advection speeds.

5.1. Decomposition of the studied domain.

Symmetries in c and ξ preserve the characteristic and minimal polynomials, reducing the analysis to −1 ≤ tv, tw ≤ 1, 0 ≤ c ≤ 1, and 0 ≤ ξ ≤ π.The proof exploits invariance under changing the signs of c or ξ.Source: §5.1 · p. 12

5.2. Technical lemmas.

The proof characterizes the interior stability region through the sign of the quadratic Υ(t, s), whose roots define its stability boundary.For −1 ≤ t < 3−2√2, Υt has a single root t⋆w(t); for 3−2√2 < t ≤ 1, Υt(s) ≥ 0 exactly when s ≤ t⋆w(t) or s ≥ t⋄w(t).Source: §5.2 · pp. 16, 19
For |tv| < 1 and |tw| < 1, the scheme is stable at c = 0 or 1; for 0 < c < 1, it is stable exactly when −1 < tw ≤ t⋆w(tv).The threshold function t⋆w is decreasing and bijective from [−1, 1] onto [−1, 1].Source: §5.2 · pp. 18, 20
For (tv, tw) ∈ Itvtw, 1(ξ,c,tv,tw) ≥ 0 is equivalent to −1 < tw ≤ t*_w(tv) at c=0 and for 0<c<1.Source: §5.2 · pp. 21–22

5.4. Stability on Etvtw

On edges, the scheme is stable for all c when tw=-1 or tv=-1; when tw=1 or tv=1, it is stable only at c=0 or c=1.Source: §5.2 · pp. 22–24

5.5. Stability in Vtvtw

At the vertices, the scheme is stable for all c at (tv, tw)=(-1, -1), only at c=0 at (1, -1), for 0<c<1 at (-1, 1), and unstable for all c at (1, 1).Source: §5.2 · pp. 25–26

5.2. Technical lemmas.

With (tv, tw) = (−1, 1), the scheme is stable for c ∈ (0, 1) and unstable at c = 0 and c = 1.For 0 < c < 1, the amplification matrix is diagonalizable and its iteration matrix is uniformly power bounded.Source: §5.2 · pp. 27–29
With (tv, tw) = (1, 1), the scheme is unstable for every c ∈ [0, 1].The proof establishes instability separately at c = 0 and c = 1, and for interior c by finding a frequency with β < 0.Source: §5.2 · pp. 29–31
For the technical parameter α, the proof bounds its values by a triangle inside the radius-3 circle, with |α| = 3 only at c = 0 and ξ = π.This geometric bound follows from the stated triangle constraint on α and identifies the sole point where it reaches the circle.Source: §5.2 · pp. 31–32

6. Convergence and discrete waves

The scheme converges for smooth data under stated stability conditions and equilibrium initialization of non-conserved moments; its three modes have distinct propagation and damping behavior.

Convergence and discrete waves

For |c| ≤1, 0 < sv < 2, and 0 < sw ≤ s⋆w(sv), the conserved moment converges for smooth initial data with non-conserved moments initialized at equilibrium.The convergence theorem assumes u◦∈Hσ(R) with σ ≥3.Source: §6 · pp. 32–33
The scheme is second-order accurate, while full second-order convergence for non-smooth data requires sufficient initial smoothness; the expected L2 order is O(∆x^min(2,2/3σ)).The reported order follows from the Sobolev or Besov regularity stated in the passage.Source: §6 · pp. 32–33
The three-mode representation expresses the numerical solution as a superposition whose initialization-dependent coefficients shape propagation and damping, especially at low frequencies.The physical eigenvalue zφ(ξ) satisfies zφ(ξ)=1−icξ+o(ξ) near ξ=0.Source: §6 · pp. 32–33
The first parasitic mode is geometrically damped for 0 < sv < 2 and propagates at velocity −c; the second is damped for 0 < sw < 2 and does not propagate.When sw = 2−sv, the second mode can carry only constants.Source: §6 · pp. 33–34

7. Numerical experiments

The experiments illustrate physical and parasitic wave behavior, test convergence across data regularities, and examine a nonlinear Burgers adaptation. They report empirical second-order convergence, including for the Burgers scheme before blowup.

7.1. Discrete waves

Wave simulations show the physical wave dispersing, a slowly damped parasitic mode traveling at −c, and a stationary spurious mode visible in an unstable configuration.Relaxation parameters near one damp the parasitic modes quickly, making the scheme behave essentially like Lax-Wendroff.Source: §7.1 · p. 34, Figure 14 · p. 34

7.2. Convergence and order

For the tested initial data, empirical convergence agrees with the predicted regularity trend and reaches full second-order accuracy starting at σ ≥3.The comparison uses simulations with 2000 and 4000 grid points.Source: §7.2 · p. 35, Figure 15 · p. 35

7.3. Empirical convergence for the Burgers’ equation

Nonlinear equilibria give a formally second-order scheme for inviscid Burgers, but the experiments do not establish stability.The scheme can be unstable even when its linearization is stable.Source: §7.3 · p. 36
Empirical second-order convergence holds for all tested relaxation parameters on inviscid Burgers before blowup; sv=sw=11/10 gives smaller asymptotic errors than other choices, including Lax-Wendroff.The simulations use periodic boundaries and end at time 1/2, slightly before blowup.Source: Figure 16 · p. 35, §7.3 · p. 36

8. Conclusion

The paper characterizes stability for a three-velocity lattice Boltzmann scheme and establishes convergence under smooth-data and initialization conditions. Numerical experiments support the analysis, while extensions to broader models and settings remain future directions.

Conclusion

The stability conditions characterize when the amplification matrix is uniformly power-bounded across Fourier frequencies, using the Courant number and two relaxation parameters.The analysis combines von Neumann analysis, eigenvalue multiplicities, and the Cohn-Schur algorithm.Source: §8 · pp. 36–37
Convergence of the conserved moment is established for smooth initial data with equilibrium initialization of non-conserved moments, and experiments support the theoretical findings.The experiments also illustrate wave propagation and second-order convergence, with additional tests on non-smooth data and an inviscid Burgers adaptation.Source: §8 · p. 37
Future work is proposed for schemes with more velocities, higher spatial dimensions, and more complex target equations.These extensions aim to broaden the framework for stability analysis and transport applications.Source: §8 · p. 37

Appendix A. Expression of the polynomials

The appendix expresses the amplification matrix characteristic polynomial and related polynomials used in the stability analysis.

Appendix A. Expression of the polynomials

The characteristic polynomial φ3 of the amplification matrix G is written as a cubic in X, with coefficients depending on ξ, c, tv, and tw.The passage gives the polynomial as a sum of coefficient terms through X^3.Source: Appendix A · p. 37
Applying the Cohn-Schur algorithm to φ3 yields reduced polynomials, including quadratic forms and related derivative expressions.The appendix provides expressions for φ◦2, φ•2, and related terms.Source: Appendix A · pp. 37–39
The appendix lists first-derivative polynomial expressions and coefficient terms, including Re(α3,2).The displayed coefficient expressions depend on ξ, c, tv, and tw.Source: Appendix A · pp. 37–39