Giacomo Dimarco

G. Dimarco, L. Pareschi

We discuss Implicit-Explicit (IMEX) Runge Kutta methods which are particularly adapted to stiff kinetic equations of Boltzmann type. We consider both the case of easy invertible collision operators and the challenging case of Boltzmann collision operators. We give sufficient conditions in order that such methods are asymptotic preserving and asymptotically accurate. Their monotonicity properties are also studied. In the case of the Boltzmann operator, the methods are based on the introduction of a penalization technique for the collision integral. This reformulation of the collision operator permits to construct penalized IMEX schemes which work uniformly for a wide range of relaxation times avoiding the expensive implicit resolution of the collision operator. Finally we show some numerical results which confirm the theoretical analysis.

Giacomo Dimarco, Lorenzo Pareschi

We introduce a class of exponential Runge-Kutta integration methods for kinetic equations. The methods are based on a decomposition of the collision operator into an equilibrium and a non equilibrium part and are exact for relaxation operators of BGK type. For Boltzmann type kinetic equations they work uniformly for a wide range of relaxation times and avoid the solution of nonlinear systems of equations even in stiff regimes. We give sufficient conditions in order that such methods are unconditionally asymptotically stable and asymptotic preserving. Such stability properties are essential to guarantee the correct asymptotic behavior for small relaxation times. The methods also offer favorable properties such as nonnegativity of the solution and entropy inequality. For this reason, as we will show, the methods are suitable both for deterministic as well as probabilistic numerical techniques.

Giacomo Dimarco, Lorenzo Pareschi

We consider the development of high order asymptotic-preserving linear multistep methods for kinetic equations and related problems. The methods are first developed for BGK-like kinetic models and then extended to the case of the full Boltzmann equation. The behavior of the schemes in the Navier-Stokes regime is also studied and compatibility conditions derived. We show that, compared to IMEX Runge-Kutta methods, the IMEX multistep schemes have several advantages due to the absence of coupling conditions and to the greater computational efficiency. The latter is of paramount importance when dealing with the time discretization of multidimensional kinetic equations.

Giacomo Dimarco, Lorenzo Pareschi

In this note we discuss the construction of high order asymptotic preserving numerical schemes for the Boltzmann equation. The methods are based on the use of Implicit-Explicit (IMEX) Runge-Kutta methods combined with a penalization technique recently introduced in [F. Filbet, S. Jin: A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources,J. Comp. Phys. 229, (2010), pp. 7625-7648.].

Giacomo Dimarco, Russell Caflisch, Lorenzo Pareschi

We consider the development of Monte Carlo schemes for molecules with Coulomb interactions. We generalize the classic algorithms of Bird and Nanbu-Babovsky for rarefied gas dynamics to the Coulomb case thanks to the approximation introduced by Bobylev and Nanbu (Theory of collision algorithms for gases and plasmas based on the Boltzmann equation and the Landau-Fokker-Planck equation, Physical Review E, Vol. 61, 2000). Thus, instead of considering the original Boltzmann collision operator, the schemes are constructed through the use of an approximated Boltzmann operator. With the above choice larger time steps are possible in simulations; moreover the expensive acceptance-rejection procedure for collisions is avoided and every particle collides. Error analysis and comparisons with the original Bobylev-Nanbu (BN) scheme are performed. The numerical results show agreement with the theoretical convergence rate of the approximated Boltzmann operator and the better performance of Bird-type schemes with respect to the original scheme.

Giacomo Dimarco, Lorenzo Pareschi

The development of efficient numerical methods for kinetic equations with stochastic parameters is a challenge due to the high dimensionality of the problem. Recently we introduced a multiscale control variate strategy which is capable to accelerate considerably the slow convergence of standard Monte Carlo methods for uncertainty quantification. Here we generalize this class of methods to the case of multiple control variates. We show that the additional degrees of freedom can be used to improve further the variance reduction properties of multiscale control variate methods.

Giacomo Dimarco

In this work we propose a generalization of the Moment Guided Monte Carlo method developed in [11]. This approach permits to reduce the variance of the particle methods through a matching with a set of suitable macroscopic moment equations. In order to guarantee that the moment equations provide the correct solutions, they are coupled to the kinetic equation through a non equilibrium term. Here, at the contrary to the previous work in which we considered the simplified BGK operator, we deal with the full Boltzmann operator. Moreover, we introduce an hybrid setting which permits to entirely remove the resolution of the kinetic equation in the limit of infinite number of collisions and to consider only the solution of the compressible Euler equation. This modification additionally reduce the statistical error with respect to our previous work and permits to perform simulations of non equilibrium gases using only a few number of particles. We show at the end of the paper several numerical tests which prove the efficiency and the low level of numerical noise of the method.

Wei Chen, Giacomo Dimarco, Lorenzo Pareschi

Plasma kinetic equations exhibit pronounced sensitivity to microscopic perturbations in model parameters and data, making reliable and efficient uncertainty quantification (UQ) essential for predictive simulations. However, the cost of uncertainty sampling, the high-dimensional phase space, and multiscale stiffness pose severe challenges to both computational efficiency and error control in traditional numerical methods. These aspects are further emphasized in presence of collisions where the high-dimensional nonlocal collision integrations and conservation properties pose severe constraints. To overcome this, we present a variance-reduced Monte Carlo framework for UQ in the Vlasov--Poisson--Landau (VPL) system, in which neural network surrogates replace the multiple costly evaluations of the Landau collision term. The method couples a high-fidelity, asymptotic-preserving VPL solver with inexpensive, strongly correlated surrogates based on the Vlasov--Poisson--Fokker--Planck (VPFP) and Euler--Poisson (EP) equations. For the surrogate models, we introduce a generalization of the separable physics-informed neural network (SPINN), developing a class of tensor neural networks based on an anisotropic micro-macro decomposition, to reduce velocity-moment costs, model complexity, and the curse of dimensionality. To further increase correlation with VPL, we calibrate the VPFP model and design an asymptotic-preserving SPINN whose small- and large-Knudsen limits recover the EP and VP systems, respectively. Numerical experiments show substantial variance reduction over standard Monte Carlo, accurate statistics with far fewer high-fidelity samples, and lower wall-clock time, while maintaining robustness to stochastic dimension.

Giacomo Dimarco, Cory Hauck, Raphaël Loubère

In this paper we genealize the fast semi-Lagrangian scheme developed in [J. Comput. Phys., Vol. 255, 2013, pp 680-698] to the case of high order reconstructions of the distribution function. The original first order accurate semi-Lagrangian scheme is supplemented with polynomial reconstructions of the distribution function and of the collisional operator leading to an effective high order accurate numerical scheme for all regimes, from extremely rarefied gas to highly collisional siuation. The main idea relies on updating at each time step the extreme points of the distribution function for each velocity of the lattice instead of updating the solution in the cell centers, these extremes points being located at different positions for any fixed velocity of the lattice. The result is a class of scheme which permits to preserve the structure of the solution over very long times compared to existing schemes from the literature. We propose a proof of concept of this new approach along with numerical tests and comparisons with classical numerical methods.

Giacomo Dimarco, Axel Klar, Theresa Köfler et al.

This work introduces a novel high-order numerical framework for solving kinetic equations, designed to remain uniformly valid across all regimes of the mean free path, spanning from the rarefied kinetic scale to the incompressible hydrodynamic limit. The method is built upon a micro-macro decomposition, which reformulates the underlying kinetic equation into a coupled system consisting of a macroscopic part, representing the fluid-dynamic evolution, and a microscopic part, describing the non-equilibrium deviations. The proposed framework ensures high-order temporal accuracy through the use of Implicit-Explicit Runge-Kutta methods, which provide stability and efficiency in stiff regimes, while spatial resolution is enhanced by combining finite-difference WENO reconstructions with high-order central difference approximations. A key feature of the proposed methodology is its Asymptotic-Preserving (AP) property. We demonstrate that, in the appropriate asymptotic limit as the mean free path tends to zero, the scheme consistently reduces to a high-order finite-difference formulation of the incompressible Navier-Stokes equations. To support the theoretical findings, a set of numerical experiments are performed on one- and two-dimensional benchmark problems, which confirm the accuracy, stability, and versatility of the method across different flow regimes.

Wei Chen, Giacomo Dimarco, Lorenzo Pareschi

The high dimensionality of kinetic equations with stochastic parameters poses major computational challenges for uncertainty quantification (UQ). Traditional Monte Carlo (MC) sampling methods, while widely used, suffer from slow convergence and high variance, which become increasingly severe as the dimensionality of the parameter space grows. To accelerate MC sampling, we adopt a multiscale control variates strategy that leverages low-fidelity solutions from simplified kinetic models to reduce variance. To further improve sampling efficiency and preserve the underlying physics, we introduce surrogate models based on structure and asymptotic preserving neural networks (SAPNNs). These deep neural networks are specifically designed to satisfy key physical properties, including positivity, conservation laws, entropy dissipation, and asymptotic limits. By training the SAPNNs on low-fidelity models and enriching them with selected high-fidelity samples from the full Boltzmann equation, our method achieves significant variance reduction while maintaining physical consistency and asymptotic accuracy. The proposed methodology enables efficient large-scale prediction in kinetic UQ and is validated across both homogeneous and nonhomogeneous multiscale regimes. Numerical results demonstrate improved accuracy and computational efficiency compared to standard MC techniques.

Giacomo Dimarco, Federica Ferrarese, Lorenzo Pareschi

In this work, we propose a data augmentation strategy aimed at improving the training phase of neural networks and, consequently, the accuracy of their predictions. Our approach relies on generating synthetic data through a suitable compartmental model combined with the incorporation of uncertainty. The available data are then used to calibrate the model, which is further integrated with deep learning techniques to produce additional synthetic data for training. The results show that neural networks trained on these augmented datasets exhibit significantly improved predictive performance. We focus in particular on two different neural network architectures: Physics-Informed Neural Networks (PINNs) and Nonlinear Autoregressive (NAR) models. The NAR approach proves especially effective for short-term forecasting, providing accurate quantitative estimates by directly learning the dynamics from data and avoiding the additional computational cost of embedding physical constraints into the training. In contrast, PINNs yield less accurate quantitative predictions but capture the qualitative long-term behavior of the system, making them more suitable for exploring broader dynamical trends. Numerical simulations of the second phase of the COVID-19 pandemic in the Lombardy region (Italy) validate the effectiveness of the proposed approach.

G. Dimarco, L. Pareschi, G. Samaey

We develop a new Monte Carlo method that solves hyperbolic transport equations with stiff terms, characterized by a (small) scaling parameter. In particular, we focus on systems which lead to a reduced problem of parabolic type in the limit when the scaling parameter tends to zero. Classical Monte Carlo methods suffer of severe time step limitations in these situations, due to the fact that the characteristic speeds go to infinity in the diffusion limit. This makes the problem a real challenge, since the scaling parameter may differ by several orders of magnitude in the domain. To circumvent these time step limitations, we construct a new, asymptotic-preserving Monte Carlo method that is stable independently of the scaling parameter and degenerates to a standard probabilistic approach for solving the limiting equation in the diffusion limit. The method uses an implicit time discretization to formulate a modified equation in which the characteristic speeds do not grow indefinitely when the scaling factor tends to zero. The resulting modified equation can readily be discretized by a Monte Carlo scheme, in which the particles combine a finite propagation speed with a time-step dependent diffusion term. We show the performance of the method by comparing it with standard (deterministic) approaches in the literature.

Giacomo Dimarco, Lorenzo Pareschi, Mattia Zanella

Kinetic equations play a major rule in modeling large systems of interacting particles. Recently the legacy of classical kinetic theory found novel applications in socio-economic and life sciences, where processes characterized by large groups of agents exhibit spontaneous emergence of social structures. Well-known examples are the formation of clusters in opinion dynamics, the appearance of inequalities in wealth distributions, flocking and milling behaviors in swarming models, synchronization phenomena in biological systems and lane formation in pedestrian traffic. The construction of kinetic models describing the above processes, however, has to face the difficulty of the lack of fundamental principles since physical forces are replaced by empirical social forces. These empirical forces are typically constructed with the aim to reproduce qualitatively the observed system behaviors, like the emergence of social structures, and are at best known in terms of statistical information of the modeling parameters. For this reason the presence of random inputs characterizing the parameters uncertainty should be considered as an essential feature in the modeling process. In this survey we introduce several examples of such kinetic models, that are mathematically described by nonlinear Vlasov and Fokker--Planck equations, and present different numerical approaches for uncertainty quantification which preserve the main features of the kinetic solution.

Giacomo Dimarco, Raphaël Loubère, Jacek Narski et al.

In this paper we deal with the extension of the Fast Kinetic Scheme (FKS) [J. Comput. Phys., Vol. 255, 2013, pp 680-698] originally constructed for solving the BGK equation, to the more challenging case of the Boltzmann equation. The scheme combines a robust and fast method for treating the transport part based on an innovative Lagrangian technique supplemented with fast spectral schemes to treat the collisional operator by means of an operator splitting approach. This approach along with several implementation features related to the parallelization of the algorithm permits to construct an efficient simulation tool which is numerically tested against exact and reference solutions on classical problems arising in rarefied gas dynamic. We present results up to the $3$D$\times 3$D case for unsteady flows for the Variable Hard Sphere model which may serve as benchmark for future comparisons between different numerical methods for solving the multidimensional Boltzmann equation. For this reason, we also provide for each problem studied details on the computational cost and memory consumption as well as comparisons with the BGK model or the limit model of compressible Euler equations.