In this talk, we begin with the nonlinear Schrödinger equation (NLSE) with different low regularity potentials and nonlinearities arising from modeling and simulation for quantum physics and chemistry, nonlinear/quantum optics, and quantum information and computation, etc. Optimal error bounds for time-splitting methods and exponential wave integrators are established for the NSLE under the proper regularity assumption on its solution determined by the low regularity potential and nonlinearity. Then we propose a novel symmetric and explicit Gautschi-type exponential wave integrator (sEWI) for the NLSE with low regularity potential and nonlinearity, and establish its optimal error estimates under various regularity assumptions on potential and nonlinearity. Extensions to the NLSE with singular potentials and nonlinearities are presented. Finally, extensions to other dispersive PDEs with low regularity potential and nonlinearity are discussed. 

This talk is based on joint works with Remi Carles, Yue Feng, Bo Lin, Ying Ma, Chunmei Su, Qinglin Tang and Chushan Wang.

Gradient Random Walk methods are particle-based techniques originally developed for the simulation of reaction–diffusion systems, in which particles evolve to represent the spatial derivatives of the solution. Inspired by vortex methods for the Navier–Stokes equations, these approaches attracted considerable interest in the 1990s because of several appealing features: they operate without a computational grid, adapt naturally by concentrating particles in regions of steep gradients, and can significantly reduce statistical variance compared with classical Random Walk schemes.

In this talk, we revisit this framework and show how its underlying ideas can be extended beyond reaction–diffusion systems to hyperbolic conservation laws. We first formulate a classical Monte Carlo approach for relaxation approximations of hyperbolic problems. We then introduce a new particle dynamics based on the evolution of spatial derivatives, leading to a gradient-based representation of the solution. When coupled with an asymptotic-preserving splitting strategy, this construction yields a new class of Gradient-Based Monte Carlo methods specifically designed for hyperbolic equations. We will present the main analytical and computational ideas behind the method and illustrate how the gradient formulation enhances adaptivity, reduces statistical fluctuations, and improves the representation of sharp solution features.

In the talk I will review stochastic particle approximations of Boltzmann-like equations and present a novel quantitative error analysis. The convergence result is based on a coupling tailored for the collisional dynamics and on Wasserstein concentration inequalities. I will also discuss extensions to particle-based methods in optimization and sampling. Joint work with Lorenzo Pareschi.

Several systems of evolutionary partial differential equations may contain stiff terms, which require an implicit treatment. Typical examples are hyperbolic systems with stiff hyperbolic or parabolic relaxation characterized by a relaxation parameter ε. In the hyperbolic-to-hyperbolic relaxation (HSHR) a natural treatment consists in adopting implicit-explicit (IMEX) schemes, in which the relaxation is treated by an implicit scheme, while the hyperbolic part is treated explicitly [1]. 

In the hyperbolic-to-parabolic relaxation (HSPR) standard IMEX methods relax to an explicit scheme for the parabolic limit, thus suffering from parabolic CFL restriction. In [2, 3] this drawback has been overcame by a penalization method, so that the system appears as the limit relaxed system plus a small perturbation. Furthermore, in [4] a unified IMEX approach has been introduced for systems which may admit both limits. This latter approach generalizes the two ones: HSHR and HSPR. All these approaches are capable to capture the correct asymptotic limit of the system when ε → 0, i.e., the scheme is asymptotic preserving (AP) independently of the scaling used. However, the AP property guarantees only the consistency of the scheme in the stiff limit ε → 0, but it does not not imply in general that the scheme preserves the order of accuracy in time in the limit and the order of accuracy

may drop to low orders. In the literature of hyperbolic system with stiff relaxation this order reduction phenomenon is extensively studied and cured, see for example [1]. A scheme that preserves the order of accuracy in time in the limit is said asymptotically accurate (AA). In this talk we show that under several assumptions on the IMEX scheme, all the numerical approaches used to solve hyperbolic systems with stiff relaxation are both AP and AA, i.e., they maintain the correct order of accuracy of the original IMEX scheme in the limit of the relaxation parameter ε.

References

[1] S. Boscarino, L. Pareschi, G. Russo, Implicit-Explicit methods for evolutionary partial differential equations, 

https://doi.org/10.1137/1.9781611978209 SIAM books, Book Series Name: Mathematical Modeling and Computation, 2024.

[2] S. Boscarino, L. Pareschi, G. Russo, Implicit-Explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit, SIAM J. SCI. COMPUT., Vol. 35, No. 1, pp. A22-A51.

[3] S. Boscarino, G. Russo, Flux-Explicit IMEX Runge-Kutta schemes for hyperbolic to parabolic relaxation problems, SIAM J. NUMER. ANAL., Vol. 51, No. 1, pp. 163-190.

[4] S. Boscarino, L. Pareschi, G. Russo, A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation SIAM J. NUMER. ANAL., Vol. 55, No. 4, pp. 2085-2109.

The Stein Variational Gradient Descent method is a variational inference method in statistics that has recently received a lot of attention. The method provides a deterministic approximation of the target distribution, by introducing a nonlocal interaction with a kernel. Despite the significant interest, the exponential rate of convergence for the continuous method has remained an open problem, due to the difficulty of establishing the related so-called Stein-log-Sobolev inequality. Here, we prove that the inequality is satisfied for each space dimension and every kernel whose Fourier transform has a quadratic decay at infinity and is locally bounded away from zero and infinity. Moreover, we construct weak solutions to the related PDE satisfying exponential rate of decay towards the equilibrium. The main novelty in our approach is to interpret the Stein-Fisher information, also called the squared Stein discrepancy, as a duality pairing between H−1(Rn) and H−1(Rn), which allows us to employ the Fourier transform. We also provide several examples of kernels for which the Stein-log-Sobolev inequality fails, partially showing the necessity of our assumptions.

Systems of collectively moving particles exhibit intringuing emergent phenomena. For instance, in the Vicsek model (in which particles moving with unit speed try to align with each other up to some noise), large-density bands of highly-aligned particles may form and travel through a low density background of disordered particles. The mechanisms underpinning this phenomenon are still poorly understood. In this talk, we show that a macroscopic model rigorously derived from the Vicsek model, the self-organized hydrodynamic model (SOH), exhibits travelling band formation. We rigorously show that travelling-bands are among the possible travelling-wave solutions of this model in one-dimension, which opens the way to a new understanding of travelling band formation in the Vicsek model.

Numerous physical, chemical, biological, and social dynamic processes are characterised by convergence to long-time equilibria. These are often described as PDEs of kinetic type, whereby `position' and `velocity' are independent variables. These may also arise when modelling multi-agent interacting processes of particles, individuals, etc. In many important cases the diffusion/dissipation required to arrive to such equilibria is explicitly present in some of the spatial directions only, that is there exist evolution PDEs with degenerate diffusion yet converging to equilibrium states as time goes to infinity. This, somewhat counter-intuitive at first, state of affairs suggests that decay to equilibrium is due to finer hidden structure, which allows for the transport terms to also `propagate dissipation’  into the spatial directions where no dissipation appears explicitly in the PDE model. Villani coined the term ‘hypocoercivity' to describe this phenomenon in his celebrated 2009 AMS Memoir, in which he presented abstract sufficient conditions for this property for general classes of PDEs. The question of preserving hypocoercivity structures is highly relevant also in the context of numerical methods for kinetic PDEs, since the trend to equilibrium is often non-monotone and, typical numerical methods are not designed specifically for `long-time' simulations. Taking as a representative example the inhomogeneous (degenerate) Fokker-Planck equation, I will present a framework on how to design numerical methods that provably preserve the same hypocoercivity structures as the equation itself. Of key importance in the construction is the treatment of lack of sufficiently high regularity of standard Galerkin/finite element spaces and how this is circumvented by careful definition of consistent stabilisation mechanisms, thereby, arriving naturally to the proof of hypocoercivity of the numerical solution in appropriate exponentially weighted function spaces. Time-permitting, some comments on the computational complexity of the new class of methods for full kinetic models will be given. I plan to conclude with a series of numerical experiments showcasing the practicality of the proposed numerical methodology.

Fractional differential equations arise in many applications to model systems with memory effects. Their nonlocal nature, however, prevents the direct use of standard adaptive ODE and DAE solvers. In this talk, we present a general framework that transforms a broad class of fractional differential equations into enlarged systems of ordinary differential equations by approximating the fractional kernel as a sum of exponentials.

The resulting systems are typically very large and stiff. We show how their special arrow-type Jacobian structure can be exploited to solve the associated linear systems with computational complexity that grows only linearly with the number of auxiliary variables. This makes it possible to employ existing efficient adaptive stiff integrators, such as Radau5 and DASSL, without modifying their numerical core.

We discuss both existing and newly developed sum-of-exponentials approximations of the fractional kernel and investigate parameter choices that lead to efficient, practical implementations. Numerical experiments demonstrate the accuracy and efficiency of the proposed approach on representative test problems. Finally, we illustrate how widely available stiff ODE solvers can be turned into effective tools for the numerical solution of fractional differential equations through this reformulation.

Many numerical methods relying on particle approximation most prominently optimization, Monte-Carlo, high-dimensional ODE integration methods, or filtering methods. In this talk we investigate their kinetic and mean field limits. The later have been used successfully to investigate convergence and stability properties. In this talk we aim to cover Metropolis Monte Carlo methods, Ensemble Kalman Filter methods, and consensus-based optimization methods. Numerical examples illustrate the theoretical findings. 

In this talk, I will review the development of fast Fourier spectral methods for a broad class of nonlinear Boltzmann equations, including models with general collision kernels, multispecies interactions, and quantum collisions, highlighting the pioneering contributions of Professor Lorenzo Pareschi that have shaped this area of research. Finally, I will present recent advances in extending these methods to kinetic models arising in fusion reaction systems. 

In this talk, we discuss a numerical framework for solving diffusively scaled kinetic equations from the rarefied kinetic scale to the incompressible hydrodynamic limit. The method is built upon a micro-macro decomposition reformulating the kinetic equation into a coupled system consisting of a macroscopic and a microscopic part, describing the non-equilibrium deviations. The proposed framework ensures high-order temporal accuracy through the use of GSA IMEX Runge–Kutta schemes, 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. In the asymptotic limit, as the mean free path tends to zero, the scheme consistently reduces to a high-order finite-difference formulation of the incompressible thermal Navier–Stokes equations. 

The theoretical findings are supported by a variety of numerical examples.

Joint work with G. Dimarco, T. Köfler, L. Pareschi, S. Tiwari.

PDE-constrained optimization has become a powerful framework for addressing inverse and control problems in systems governed by partial differential equations (PDEs). Kinetic theory encompasses a broad class of equations that describe the non-equilibrium dynamics of interacting particles, and applying PDE-constrained optimization to these systems reveals a variety of unique and intriguing behaviors.
In this talk, I will discuss two representative cases. In the first, we explore the use of PDE-constrained optimization for stabilizing plasma instabilities—a challenge central to achieving controlled fusion. Dispersion-relation-based linear analysis and landscape analysis are employed to identify optimal stabilization strategies. In the second case, we highlight an essential yet often overlooked adjustment required for gradient computation when PDEs are solved using particle methods.

The study of dynamical systems defined on networks (graphs) has a long history of results and applications. Early work linked graph spectra and Laplacian operators to diffusion andstability analyses; the Kuramoto model and subsequent reductions (e.g., Ott–Antonsen) catalyzed rigorous study of synchronization in coupled oscillators. The emergence of complex-network models (e.g. Watts–Strogatz, Barabási–Albert models) and epidemic-network studies highlighted how topology shapes collective dynamics and thresholds. Methods for stability and control, notably the master stability function, clarified when and how networks synchronize or desynchronize. More recently, graph limit theory (graphons) provided a principled passage from discrete networks to nonlocal continuum equations, enabling mean-field descriptions and rigorous convergence results for large dense graphs. Complementary frameworks (graphexes, sparse limits) extend these ideas to sparse regimes. These advances support contemporary work connecting continuum limits to graph neural networks, reservoir computing, neuroscience, ecology, and epidemic modeling. We will review some of these applications and highlight several open challenges, including inferring graphons from data, extending continuum techniques to sparse and heterogeneous graphs, and translating theoretical limits into practical tools for the control and learning of large-scale networks.

This work was developed in collaboration with G. Aletti (UniMI), G. Patanè (IMATI.CNR), A. Marta (UniMI), and M. Calabrò (PoliMI).

We are interested in the development of hydrodynamic equations from kinetic models, including the evolution of entropy. It turns out that it is possible to derive from kinetic equations the structure of symmetric hyperbolic thermodynamically compatible (SHTC) models proposed by Godunov-Romenski and used by Dumbser and coworkers. This derivation connects the entropy dissipation provided by the SHTC model with non equilibrium effects in kinetic theory. 

Another development in the same line concerns the derivation of Baer-Nunziato like models for two phase flows from kinetic equations.

Joint work with: Thomas Rey and Tommaso Tenna

We present a Multilevel Kinetic-Diffusion Monte Carlo (ML-KDMC) method for the simulation of kinetic Boltzmann transport equations with a BGK collision operator, with particular motivation from neutral-particle transport in the plasma edge of nuclear fusion reactors. The underlying Kinetic-Diffusion Monte Carlo method remains accurate in both low- and high-collisional regimes while avoiding the prohibitive cost growth typically associated with the diffusive limit. We show that embedding KDMC in a Multilevel Monte Carlo framework, based on a hierarchy of increasingly large time steps, further reduces the computational cost. A central challenge in fusion applications is the accurate estimation of moments of the particle distribution through event scoring along particle trajectories. Since velocity information is lost during the diffusive step of KDMC, this traditionally requires dedicated estimators. We demonstrate that such estimators can be avoided by reinterpreting KDMC as a distribution-composition hybrid fluid-kinetic method.

This work is joint with Emil Loevbak, Bert Mortier, Pieterjan Robbe, and Zhirui Tang.

The BGK equation was formulated to be a numerically tractable model of the Boltzmann equation. More than 70 years after its publication, and in spite of the enormous increase in computational power since then, it is still being used and attracts quite some attention from the mathematical community. It is phenomenological in the sense that it preserves many of the properties of the Boltzmann equation, but there is no direct connection with the underlaying particle systems. The BGK has been obtained as the limit of random particle systems, but here I will show how the BGK equation arises from a deterministic, reversible particle system, a modification of the hard sphere dynamics. This leads to a hard sphere version of the BGK equation, but the particle system can be modified so as to yield the usual BGK equation in the limit. In the talk I will not attempt to present any proof, but the results are illustrated with numerical simulations.