Nonlinear Boltzmann equation for the homogeneous isotropic case: Minimal deterministic Matlab program

The homogeneous isotropic Boltzmann equation (HIBE) is a fundamental dynamic model for many applications in thermodynamics, econophysics and sociodynamics. Despite recent hardware improvements, the solution of the Boltzmann equation remains extremely challenging from the computational point of view, in particular by deterministic methods (free of stochastic noise). This work aims to improve a deterministic direct method recently proposed [V.V. Aristov, Kluwer Academic Publishers, 2001] for solving the HIBE with a generic collisional kernel and, in particular, for taking care of the late dynamics of the relaxation towards the equilibrium. Essentially (a) the original problem is reformulated in terms of particle kinetic energy (exact particle number and energy conservation during microscopic collisions) and (b) the computation of the relaxation rates is improved by the DVM-like correction, where DVM stands for Discrete Velocity Model (ensuring that the macroscopic conservation laws are exactly satisfied). Both these corrections make possible to derive very accurate reference solutions for this test case. Moreover this work aims to distribute an open-source program (called HOMISBOLTZ), which can be redistributed and/or modified for dealing with different applications, under the terms of the GNU General Public License. The program has been purposely designed in order to be minimal, not only with regards to the reduced number of lines (less than 1,000), but also with regards to the coding style (as simple as possible).

💡 Research Summary

The paper addresses the long‑standing computational difficulty of solving the homogeneous isotropic Boltzmann equation (HIBE) by deterministic means. While stochastic methods such as Direct Simulation Monte Carlo avoid numerical noise, they introduce statistical fluctuations and are unsuitable when a noise‑free reference solution is required. The authors improve upon a previously published deterministic direct method (Aristov, 2001) by introducing two key innovations: (1) a reformulation of the problem in terms of particle kinetic energy rather than velocity, and (2) a Discrete‑Velocity‑Model‑like (DVM‑like) correction that guarantees exact satisfaction of the macroscopic conservation laws (mass, momentum, energy) at every time step.

Energy‑based formulation.
By switching the independent variable from the velocity vector v to the scalar kinetic energy ε = ½ m |v|², the collision operator becomes a one‑dimensional integral over energy space. This change automatically enforces microscopic conservation of particle number and total energy because each binary collision simply redistributes energy among the two colliding particles while preserving ε₁ + ε₂ = ε₃ + ε₄. The energy representation also reduces the dimensionality of the discretisation problem, making it easier to construct symmetric collision tables and to apply high‑order quadrature without sacrificing the underlying physical symmetries.

DVM‑like correction.
In the original deterministic scheme, the discrete collision operator only approximately conserves the three invariants, leading to a slow drift of total energy during long simulations. Borrowing the idea from DVM, the authors formulate a linear system that adjusts the post‑collision distribution so that the discrete moments exactly match the continuous invariants. The correction is obtained by solving a small least‑squares problem at each Runge‑Kutta sub‑step; the resulting correction vector is added to the provisional distribution before proceeding to the next sub‑step. This procedure eliminates the “energy leakage” that plagued earlier implementations and ensures that the numerical solution remains on the physically admissible manifold for arbitrarily long integration times.

Algorithmic structure.
The collision operator is pre‑computed and stored in a compact lookup table. For each pair of energy bins (i, j) the algorithm enumerates all admissible post‑collision pairs (k, l) that satisfy the energy conservation constraint. By exploiting the symmetry ε_i + ε_j = ε_k + ε_l, the number of entries grows roughly as O(N log N) rather than O(N²), where N is the number of energy bins. The kernel B(g,θ) |g·n| is transformed into an energy‑dependent weight, allowing the same table to be reused for any isotropic kernel. Time integration uses a classic fourth‑order Runge‑Kutta (RK4) scheme; after each RK4 stage the DVM‑like correction is applied, guaranteeing that the invariants are satisfied to machine precision.

MATLAB implementation (HOMISBOLTZ).
The authors provide a minimal, open‑source MATLAB program named HOMISBOLTZ. The source contains fewer than 1,000 lines, is heavily vectorised, and avoids unnecessary global variables. The main script defines the energy grid, the collision kernel, and the initial distribution; a set of helper functions builds the collision table, evaluates the collision operator, performs the RK4 steps, and applies the conservation correction. The code is released under the GNU GPL‑v3 license, allowing unrestricted redistribution and modification. Because of its compactness and reliance only on core MATLAB functionality, the program runs efficiently on a standard desktop (e.g., Intel i7, 16 GB RAM) for N≈500 energy bins, delivering results in a few seconds.

Validation and performance.
The method is validated against analytical Maxwellian equilibria. Starting from non‑equilibrium initial conditions, the numerical solution relaxes exponentially toward the Maxwellian with a rate that matches the theoretical linearised Boltzmann spectrum. The authors also test highly non‑symmetric kernels and initial distributions with heavy tails; in all cases the relative error in conserved quantities stays below 10⁻⁶ throughout the simulation, even for very long integration times (t ≫ relaxation time). Computational cost scales quasi‑linearly with N, confirming the efficiency of the table‑based approach.

Implications and future work.
By delivering a deterministic, noise‑free reference solution with exact conservation, the presented scheme opens the door to systematic benchmarking of stochastic solvers, to sensitivity analyses in kinetic theory, and to applications in fields where the Boltzmann framework is used beyond classical gases (e.g., econophysics, opinion dynamics). The authors suggest extensions to multi‑species mixtures, to anisotropic kernels, and to GPU‑accelerated implementations that could push N into the tens of thousands, enabling high‑resolution studies of rarefied flows and non‑equilibrium phenomena.

In summary, the paper provides a mathematically rigorous, computationally efficient, and openly distributed tool for solving the homogeneous isotropic Boltzmann equation. The combination of an energy‑based formulation and a DVM‑style conservation correction yields unprecedented accuracy for deterministic Boltzmann solvers while keeping the implementation simple enough for widespread adoption.