ABSTRACT. This talk focuses on the construction and analysis of numerical schemes for computing rough solutions of nonlinear dispersive equations, such as the Korteweg–De Vries (KdV) equation and the nonlinear Schrödinger (NLS) equation, with regularity below H¹. To establish stability estimates under low-regularity conditions, we introduce continuous formulations of the numerical schemes, wherein the discrete numerical solution is recast as the solution of a perturbed continuous equation. This approach reduces the discrete stability analysis to analyzing the stability of the continuous equation with respect to perturbations. Consequently, stability can be established by utilizing continuous-level Bourgain or Strichartz estimates, thereby circumventing the restrictive CFL conditions typically required when studying the stability of numerical schemes at the discrete level. Furthermore, a high- and low-frequency splitting technique may be further employed to improve the convergence rates of the numerical solutions. For the KdV equation, we prove that the proposed method converges in L² with order β−ϵ under the regularity condition u∈C([0,T];Hᵝ) for β∈(0,1]. Similarly, for the NLS equation, the method converges in L² with order 1.5β−ϵ under the regularity condition u∈C([0,T];Hᵝ) for β∈(0,1/2].

ABSTRACT. Regularisation of turbulent flows is essential for computationally efficient forecasting. In this talk, we highlight a geometric mechanism for regularising two-dimensional turbulence that preserves key conserved quantities. We begin with a high-level introduction to how energy is distributed across scales in turbulent flows, using energy spectra to quantify turbulence intensity and to illustrate energy cascades that govern flow evolution. We then discuss the conserved quantities unique to two-dimensional turbulence and how they give rise to the inverse energy cascade. This sets the stage for geometric regularisation, which allows for control over the energy distribution to enable reduced-complexity computations while retaining fundamental invariants and the inverse energy cascade.

ABSTRACT. I will give an overview of the matrix hydrodynamics approach to fluids and plasmas. Discovered by V. Zeitlin in the 90s, it is the only known method to obtain discrete models for 2D fluids fully capturing the non-canonical Hamiltonian structure of fluid equations. Further, I will discuss structure preserving time integration methods for such models, in particular, Lie-Poisson integrators for Zeitlin's matrix equations. MHD-Zeitlin equations arising in plasma physics are an example of a non-isospectral matrix flow. We will see the connection of these equations to isospectral flows and why the theory developed for isospectral flows is insufficient to develop Lie-Poisson integrators for MHD-Zeitlin equations.

ABSTRACT. We discuss an evolving low-rank approximation of the vorticity solution of the Euler equations describing the flow of a two-dimensional incompressible ideal fluid on the sphere. Such problem can be approximated by the so-called Zeitlin model, an isospectral Lie--Poisson flow on the Lie algebra of traceless skew-Hermitian matrices. We propose an approximation of Zeitlin's model based on a time-dependent low-rank factorization of the vorticity matrix and evolve a basis of eigenvectors according to the Euler equations. In particular, we discuss the properties of the approximate flow and we show that the error in the solution, in the approximation of the Hamiltonian and of the Casimir functions only depends on the approximation of the vorticity matrix at the initial time. The computational complexity of solving the approximate model is shown to scale quadratically with the order of the vorticity matrix and linearly if a further approximation of the stream function is introduced.

ABSTRACT. Two-dimensional (2-D) incompressible, inviscid fluids produce fascinating patterns of swirling motion. How and why the patterns emerge are long-standing questions, first addressed in the 19th century by Helmholtz, Kirchhoff, and Kelvin. Countless researchers have since contributed to innovative techniques and results. Yet, the overarching problem of swirling 2-D motion and its long-time behavior remains largely open. Here we shed light on this problem via a link to isospectral matrix flows. The link is established through V. Zeitlin's beautiful model for the numerical discretization of Euler's equations in 2-D. When considered on the sphere, Zeitlin's model offers deep connections between 2-D hydrodynamics and unitary representations of the rotation group. Consequently, it provides a dictionary that maps hydrodynamical concepts to matrix Lie theory, which in turn gives connections to matrix factorizations, random matrices, and integrability theory, for example. Results about finite-dimensional matrices can then be transferred to infinite-dimensional fluids via quantization theory, which is here used as an analysis tool (albeit traditionally describing the limit between quantum and classical physics). We demonstrate how the dictionary is constructed and how it unveils techniques for 2-D hydrodynamics. We also give accompanying convergence results for Zeitlin's model on the sphere.

ABSTRACT. Sampling from quasi-stationary distributions (QSDs) is a key component of accelerated MD algorithms, which simulate long trajectories from metastable systems. Their efficiency, and in particular the parallel efficiency of the parallel replica algorithm, is often hindered by the QSD sampling efficiency. A crucial mechanism is the estimation of decorrelation times, measuring convergence to the QSD.

Traditional methods use either theoretically-inspired heuristics, which are only valid in certain physical regimes, or MCMC-inspired diagnostics applied to Fleming-Viot processes, which tend to display poor sensitivity, and be as a result overly conservative. Both these limitations imply that efficient deployments of previous parallel replica methods were restricted to highly metastable systems, and in particular were not applicable to moderately metastable systems common in biology.

The goal of this work is to explore an alternative to these approaches using parametric diagnostics defined by recurrent neural networks (RNN). We fit the RNN diagnostic on a synthetic ensemble of Fleming-Viot processes, driven by dynamics representative of effective molecular collective variable dynamics.

In this talk, we explain the theoretical basis for this method, discuss its implementation, and demonstrate it on increasingly complex examples, comparing its performance to that of the traditional MCMC diagnostics.

ABSTRACT. The narrow escape problem is a prototypical example to study entropic metastability, inspired by concrete applications in biology and chemistry. The problem is the determination of the exit time and position of a Brownian particle trapped in a domain with narrow holes. Unlike the usual double well potential, this form of metastability is not of energetic form as there is no potential involved, hence no Eyring-Kramers formula. Our goal is to study this problem in any dimension for arbitrary shapes of hole, and thanks to the quasi-stationary distribution approach to metastability, to rigorously derive the counterpart Eyring-Kramers type formulas in this context. Our numerical predictions are illustrated by numerical simulations, both Monte Carlo like and deterministic ones based on solving appropriate partial differential equations.

ABSTRACT. We study the mathematical validity of the recently (re-)proposed Norton formulation of non--equilibrium molecular dynamics, (Evans et al 1983, Blassel & Stoltz 2024). In the standard formulation of non--equilibrium molecular dynamics (NEMD), one perturbs a reference system by an external forcing of fixed intensity and then measures the resulting response of an observable. In the Norton formulation, the reference system is perturbed by an external forcing with varying intensity so that the value of a given observable is fixed. Numerical simulations suggest that the variance of the resulting estimator has better scaling with dimension compared to the standard NEMD estimator. Numerical evidence and a non--rigorous proof in (Evans 1993) suggest that the two formulations are equivalent, however rigorous proof of this equivalence remains open. We use a propagation of chaos to rigorously prove the equivalence of the Norton and standard formulation of NEMD for a toy model of mean--field interacting particle systems.

Joint work with Gabriel Stoltz

ABSTRACT. This talk addresses the numerical approximation of quantum statistical observables using mean-field molecular dynamics and machine learning-based surrogate models. We first discuss the approximation of canonical quantum correlation observables by classical molecular dynamics. In particular, we consider a mean-field molecular dynamics formulation in which the nuclei Hamiltonian is obtained by taking the partial trace over the electron degrees of freedom. Under this formulation, the resulting classical dynamics approximates the quantum correlation observables with error of order $\mathcal{O}(M^{-1}+\epsilon^2 t)$, where $M$ denotes the nucleus-electron mass ratio, $t$ is the correlation time, and $\epsilon^2$ is a parameter measuring the variance of the mean-field approximation.

We then present a machine-learning approach for approximating the mean-field potential using random Fourier feature neural networks. For a network with $K$ nodes trained on $J$ data samples, we show that the resulting molecular dynamics approximation of correlation observables has expected error of order $\mathcal{O}((K^{-1}+J^{-1/2})^{1/2}). Numerical experiments support the theoretical analysis and demonstrate the accuracy of the proposed framework.

ABSTRACT. We report recent results on the numerical analysis of one-dimensional linear dispersive differential equation with concentrated potential, featuring a competition between weak dispersion and localization induced by the concentrated potential. After obtaining precise regularity estimates of the exact solution in terms of the small parameters of the problem, we apply a natural first-order exponential integrator to discretize the equation, and establish optimal error bounds in terms of these small parameters. The analysis combines iterated Duhamel's expansions and cancellations in oscillatory phases that cannot be obtained directly from regularity estimates of the exact solution. We also show that other classical numerical schemes, such as Lie or centered splitting schemes and low regularity integrators, fail to display optimal rates of convergence. Numerical results confirm the theoretical error estimates. This is joint work with Chushan Wang.

ABSTRACT. We propose a fast algorithm for solving the fractional Fokker–Planck equation in arbitrary spatial dimensions. The method reformulates the nonlocal problem into a sequence of tractable subproblems, enabling efficient computation without relying on high-dimensional grid-based discretizations. In contrast to conventional approaches, the proposed algorithm avoids the curse of dimensionality and significantly reduces the cost associated with the fractional operator. The overall computational complexity scales linearly with the spatial dimension. Numerical experiments demonstrate the efficiency and scalability of the method in high-dimensional settings.

ABSTRACT. In this talk, we begin with the nonlinear Schroedinger 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.

ABSTRACT. Some of the most computationally demanding problems in fluid mechanics can be formulated as large-scale and high-dimensional matrix and tensor differential equations (MDEs and TDEs). These arise in applications ranging from turbulent combustion and probability density function transport to uncertainty quantification, sensitivity analysis, and kinetic descriptions of complex flows. Even in three-dimensional turbulent reacting flows, the multiscale structure and extreme dimensionality of the governing equations render direct numerical simulation prohibitively expensive. The underlying challenge is structural: the number of degrees of freedom grows combinatorially with dimension, resolution, or parametric complexity. Low-rank matrix approximations and their tensor-network extensions provide a principled framework for addressing this challenge by exploiting intrinsic low-dimensional structure and multidimensional correlations in evolving flow fields. In this presentation, we focus on computational strategies for solving nonlinear MDEs and TDEs directly in adaptive low-rank form. We introduce a novel formulation based on adaptive cross interpolation and develop CUR-type decompositions grounded in the Discrete Empirical Interpolation Method (DEIM) for matrix, tensor-train, and Tucker representations. These DEIM-CUR algorithms enable on-the-fly reduced-order modeling with near-optimal accuracy for broad classes of nonlinear dynamical systems. The effectiveness of these approaches is demonstrated in applications involving turbulent combustion, transport-dominated flows, uncertainty quantification, sensitivity analysis, and high-dimensional partial differential equations arising in fluid mechanics.

ABSTRACT. The linear Boltzmann equation provides a physically accurate model for charged particle transport. This is essential for applications such as radiation therapy, where dose deposition must be carefully controlled. However, solving it numerically is challenging especially for protons, due to multiscale effects, a high-dimensional phase space and strongly forward-peaked scattering. The resulting need for finely resolved phase space discretizations or specialized methodological modifications often leads to prohibitive memory requirements and computational costs.

The dynamical low-rank approximation has been shown to drastically reduce dimensionality in a variety of kinetic problems, including electron and proton transport. While orders of magnitude in cost reduction of runtime and allocated memory were achieved compared to full-rank deterministic computations, the spatial and angular resolutions required to compete with gold standard Monte Carlo reference codes exceed what is customary in application and small negative regions remain. To tackle these issues, we explore a combination of dynamical low-rank approximation with locally refined spatial discretizations. For this, we compare naive rectilinear refinement with the use of octrees and discuss refinement strategies according to the expected dose distribution as well as tissue density. We study the effect on the required rank and number of spatial cells as well as restrictions imposed on the energy discretization. First results indicate that the possibility of strong refinement in regions with steep gradients, such as near the Bragg peak of a proton beam, is very effective in eliminating unphysical behavior and allows for a significantly smaller number of spatial cells overall.

ABSTRACT. The adjoint state method is an essential tool for PDE-constrained optimization, but for high-dimensional PDEs it often suffers from the curse of dimensionality, where memory and computational cost grow exponentially with the number of dimensions. This is even a more severe problem if the computation of the gradient requires storing the full high-dimensional forward trajectory to compute the adjoint.

In this talk, we present a novel numerical scheme that overcomes this bottleneck. By employing a time-reversible dynamical low-rank integrator, we compute the forward problem and then solve both the forward and adjoint equations in tandem by evolving them backward in time. The low-rank approximation reduces the dimensionality of the forward and the adjoint solution, while the time-reversibility avoids storing the full forward trajectory. Our algorithm only requires the storage of non-linear quantities, which are typically lower-dimensional (such as the self-consistent electric field for kinetic equations in plasma physics).

We propose a strategy on how to ensure approximate time-reversibility in the presence of rank-deficiency or chaotic dynamics. Finally, we demonstrate the efficiency of our numerical method with examples from plasma physics, such as the stabilization of beam-heated plasmas or plasma beam focusing.

ABSTRACT. Dynamical low-rank approximation (DLRA) has emerged as a highly efficient framework for solving high-dimensional time-dependent partial differential equations, such as the Schrödinger equation. This talk focuses on Tucker tensors as a data-sparse tensor decomposition that enables tractable computations in high-dimensional settings. While recent developments have produced the class of Basis Update and Galerkin (BUG) integrators, these methods are generally limited by a first-order error bound in time.

We present a novel time integration scheme that achieves a proven second-order error bound. An additional key feature of the approach is its inherent parallelism: the new scheme enables the fully parallel update of all nodes within the tensor decomposition, significantly reducing computational complexity. Generalizations to tree tensor networks are discussed and illustrated by numerical experiments for many-body quantum spin systems.

ABSTRACT. I first present the Singular Euler—Maclaurin expansion, an extension of the 300-year-old classical Euler-Maclaurin summation formula to long-range interactions on high-dimensional lattices with applications in spin systems [1,2]. This method allows for the exact representation of a discrete lattice in terms of its continuous analog, with corrections given in terms of a generalization of the Riemann zeta function, the so-called Epstein zeta function. With the properties and efficient computation of this function analyzed in [4] and a high-performance implementation available in our library EpsteinLib [5], a new toolset is provided for studying arbitrary long-range interacting lattices. I briefly discuss recent extensions to systems with boundaries [6] as well as to micromagnetics systems [7]. Building on this framework, I subsequently study 2D unconventional BCS superconductors with long-range interactions, finding a rich phase diagram with topologically non-trivial and mixed-parity phases as well as stabilization of Higgs modes in the non-equilibrium dynamics [3].

In the second part of the talk, I present how generalized zeta functions, built from the Epstein zeta function, can allow for the precise evaluation of n-body interaction energies in chemistry given by (n-1) d-dimensional lattice sums, reducing the scaling from exponential to linear in the number of interaction partners n [8]. For cuboidal lattices with a two-body Lennard–Jones potential coupled to a three-body Axilrod–Teller–Muto potential, we demonstrate that increasing the three-body coupling can drive a structural transition from fcc to bcc [8,9].

In long-range-interacting quantum lattices, current graph decomposition methods, such as pCUT, transform the problem of an exponentially growing Hilbert space dimension into the computation of high-dimensional lattice sums associated with graphs, which are usually computed with low-precision Monte Carlo methods. In the final part of the talk, I present ongoing work on computing the arising graph zeta functions efficiently and precisely, laying the foundation for exploring quantum systems in regimes inaccessible to other methods.

References: [1] Singular Euler-Maclaurin expansion on multidimensional lattices, Andreas A. Buchheit and Torsten Keßler, Nonlinearity 35 3706 (2022) [2] On the Efficient Computation of Large Scale Singular Sums with Applications to Long-Range Forces in Crystal Lattices, Andreas A. Buchheit and Torsten Keßler, J. Sci. Comput. 90, 53 (2022) [3] Exact Continuum Representation of Long-range Interacting Systems and Emerging Exotic Phases in Unconventional Superconductors, Andreas A. Buchheit, Torsten Keßler, Peter K. Schuhmacher, and Benedikt Fauseweh, Phys. Rev. Research 5, 043065 (2023) [4] Computation and properties of the Epstein zeta function with high-performance implementation in EpsteinLib, Andreas A. Buchheit, Jonathan Busse, and Ruben Gutendorf, arXiv preprint 2412.16317 (2025) [5] Github Repository: github.com/epsteinlib/epsteinlib, pip install epsteinlib [6] On the computation of lattice sums without translational invariance, Andreas A. Buchheit, Torsten Keßler, and Kirill Serkh, Math. Comp. 94 (2025), 2533-2574 [7] Zeta expansion for long-range interactions under periodic boundary conditions with applications to micromagnetics, J. Comp. Phys. 559, 114885 (2026) [8] Exact lattice summations for Lennard-Jones potentials coupled to a three-body Axilrod–Teller–Muto term applied to cuboidal phase transitions, Andres Robles-Navarro, Andreas A. Buchheit, et. al., J. Chem. Phys. 163, 094104 (2025) [9] Epstein zeta method for many-body lattice sums, Andreas A. Buchheit, Jonathan K. Busse, arXiv:2504.11989 (2025), accepted in IMA Journal of Numerical Analysis ————

ABSTRACT. Computing excited states of many-body quantum systems remains a central challenge in quantum chemistry and condensed matter physics, with existing approaches broadly divided between variational and linear response methods. In this talk, I present a unified geometric framework based on Kähler manifolds that naturally encompasses both strategies for a wide class of electronic structure models, including Hartree–Fock and low-rank tensor methods in Tucker (CASSCF) or tensor-train (DMRG) formats [1,2]. I will also revisit the foundations of TDDFT through constrained Schrödinger dynamics, highlighting alternative geometric formulations and their implications for nonadiabatic modeling [3,4].

[1] L. Grazioli, Y. Hu, E. Cancès, Critical point search and linear response theory for computing electronic excitation energies of molecular systems. Part I: General framework, application to Hartree-Fock and DFT, J. Chem. Phys. [2] L. Grazioli, Y. Hu, T. Nottoli, F. Lipparini, Eric Cancès, Critical point search and linear response theory for computing electronic excitation energies of molecular systems. Part II. CASSCF, arXiv:2604.13753 [3] E. Cancès, T. Duez, J. van Gog, A.B. Lauritsen, M. Lewin, J. Toulouse, Geometric Time-Dependent Density Functional Theory, arXiv:2601.07724 [4] E. Cancès, T. Duez, J. van Gog, A.B. Lauritsen, M. Lewin, J. Toulouse, Geometric theory of constrained Schrödinger dynamics with application to time-dependent density-functional theory on a finite lattice, arXiv:2601.07719

ABSTRACT. The Ginzburg–Landau model describes complex superconducting phenomena, including the formation of vortex lattices under applied magnetic fields. In recent years, the numerical approximation of its energy minimizers has regained significant attention, yielding new insights into their structure and behavior. In this talk, we review the current state of the art in the approximation of Ginzburg–Landau minimizers using finite element methods. We also discuss the computation of minimizers via adaptive gradient descent schemes, which crucially depend on the choice of initial states. To address this challenge, we explore the use of neural networks to generate initial guesses, thereby enhancing the robustness of the optimization process.

ABSTRACT. We aim to accurately approximate ground states of the stationary Gross–Pitaevskii eigenvalue problem with rotation. These states are characterized as minimizers of the associated energy functional and may exhibit vortices at a priori unknown locations, thereby calling for adaptive resolution.

We employ a novel variational adaptivity framework that combines a discrete energy-decreasing scheme (based on an appropriate gradient flow for the GPE) with adaptive finite element discretizations. More specifically, we design an effective hp-adaptive refinement strategy driven solely by energy reduction, rather than by classical a posteriori error estimators. Numerical results indicate exponential convergence for the ground-state energy and wave function.

ABSTRACT. We consider interacting particle dynamics with Vicsek type interactions, and their macroscopic PDE limit, in the non-mean-field regime; that is, we consider the case in which each particle/agent in the system interacts only with a prescribed subset of the particles in the system (for example, those within a certain distance). In this non-mean-field regime the influence between agents (i.e. the interaction term) can be normalised either by the total number of agents in the system (global scaling) or by the number of agents with which the particle is effectively interacting (local scaling). We compare the behaviour of the globally scaled and the locally scaled systems in many respects, considering for each scaling both the PDE and the corresponding particle model. In particular we observe that both the locally and globally scaled particle system exhibit pattern formation (i.e. formation of travelling-wave-like solutions) within certain parameter regimes, and generally display similar dynamics. The same is not true of the corresponding PDE models. Indeed, while both PDE models have multiple stationary states, for the globally scaled PDE such (space-homogeneous) equilibria are unstable for certain parameter regimes, with the instability leading to travelling wave solutions, while they are always stable for the locally scaled one, which never produces travelling waves. This observation is based on a careful numerical study of the model, supported by further analysis.

Based on joint work with P. Butta, T. Hodgson, B. Goddard and K. Painter

ABSTRACT. Interacting particle systems with pairwise interactions provide a flexible framework for modeling clustering dynamics, but their direct simulation becomes computationally prohibitive at large scales. In this talk, I present two complementary coarse-graining strategies for the efficient description of clustering dynamics in such systems. First, I demonstrate how stochastic partial differential equations (SPDEs) arise as effective mean-field descriptions that capture stochastic effects beyond deterministic mean-field limits and accurately reproduce cluster formation and merging dynamics. Second, I introduce a data-driven reduction of transfer operators that yields low-dimensional Markov models describing metastable configurations and transition pathways of the particle system. Together, these approaches provide computationally efficient reduced models for interacting particle systems, enabling the analysis of long-time dynamics and large-scale behavior.

ABSTRACT. The study of opinion dynamics is essential for understanding how collective consensus and societal polarization emerge from individual interactions. In this talk, we introduce a stochastic agent-based model (ABM) that captures the co-evolution of opinion formation and social dynamics. By taking the limit of large agent numbers, we derive reduced PDE and SPDE models that characterize macroscopic phenomena, such as clustering into echo chambers, while significantly reducing computational complexity. We demonstrate our model's performance using real-world survey data on political identity and public opinion regarding governmental issues.

ABSTRACT. Mean-field stochastic differential equations, also known as McKean-Vlasov equations, arise as the limiting equations of interacting particle systems with symmetric interaction potentials in the infinite-particle limit. Such systems play an important role in a variety of fields ranging from biology and physics to sociology and economics. Global insights into the behavior of complex dynamical systems can be obtained by analyzing associated transfer operators such as the Perron-Frobenius, Koopman, and forward-backward operators and their infinitesimal generators. We extend transfer operator theory to the McKean-Vlasov setting and show how variants of the extended dynamic mode decomposition can be used to compute finite-dimensional approximations of transfer operators and generators. This framework enables us to identify metastable sets, spatiotemporal patterns, and time-evolving clusters, and also to learn interaction kernels and global potentials of the interacting particle systems.

ABSTRACT. Interacting particle systems provide a powerful modeling framework for collective dynamics in nature and engineering. While prior methods have primarily addressed symmetric interactions using various learning techniques, many real-world systems exhibit asymmetric interactions, which demand more general and flexible modeling tools. In this talk, I will present a new Sparse Bayesian Learning (SBL) framework for identifying asymmetric interaction kernels in the Motsch–Tadmor model. By reformulating the nonlinear inverse problem as a subspace identification task, we establish identifiability guarantees and enable robust kernel recovery. Incorporating informative priors, the proposed SBL algorithm offers principled model selection and uncertainty quantification, achieving reliable inference from noisy trajectory data.

ABSTRACT. Depth is widely viewed as a central contributor to the success of deep neural networks, while standard approximation theory typically provides guarantees only for the final output and leaves the role of intermediate layers unclear. This talk presents recent theoretical results that address this gap through a layer-wise approximation framework. We construct a fixed-width mixed-activation architecture in which each layer or grade produces an intermediate readout that is itself a certified approximant to the target function, with approximation errors controlled at progressively finer geometric scales. Thus, depth admits a precise mathematical interpretation as coarse-to-fine approximation refinement. Motivated by the multigrade deep learning viewpoint, these results provide a mathematical foundation for viewing deep networks as nested approximation schemes, with potential relevance to scientific machine learning for differential equations.

ABSTRACT. We propose a novel training method, Sigma-Active Neuron Least Squares (ANLS), for multivariate two-layer neural networks of smooth activation functions. Sigma-ANLS is obtained from a general iterative training framework whose core mechanism is single-neuron optimization (SNO). In particular, -ANLS is designed to exploit both locality and non-locality for SNO. The locality features are derived from the gradient-flow analysis of training. For the nonlocality, we propose stochastic heuristic strategies. Numerical examples are provided that demonstrate the effectiveness of ANLS compared with existing first- and second-order optimization algorithms on various learning tasks ranging from function approximation and dynamical systems to solving PDEs.

ABSTRACT. We present a Cartesian grid method for homogeneous and inhomogeous scattering problems on complex domains. The method is a generalization of the traditional boundary integral method. It solves the scattering problem in the framework of boundary integral method but avoids direct evaluation of boundary and volume integrals. The evaluation is done by indirectly solving equivalent interface probelms on Cartesian grids with fast solvers. For (exterior) problems on unbounded domains, we introduce an artificial circle or sphere to accelerate the solution of interface problems, while preserving accuracy. In the talk, we shall also present numerical examples to demonstrate the method.

ABSTRACT. We consider the harmonic map heat flow problem for a radially symmetric case. For discretization of this problem we apply a H1-conforming finite element method in space combined with a semi-implicit Euler time stepping. The semi-implicit Euler method results in a linear problem in each time step. We restrict to the regime of smooth solutions of the continuous problem and present an error analysis of this discretization method. This results in optimal order discretization error bounds. Key ingredients of the analysis are a discrete energy estimate, that mimics the energy dissipation of the continuous solution, and a convexity property that is essential for discrete stability and for control of the linearization error. We also present numerical results that validate the theoretical ones.

ABSTRACT. The Stefan problem is the prototypical model of phase transitions, dating back to work in the 1800s. In this talk we discuss an evolving surface finite element method (ESFEM) for solving the two-phase Stefan problem when posed on a surface evolving with prescribed evolution. In particular, we discuss the error analysis for an ESFEM discretisation of the enthalpy formulation of the two-phase Stefan problem. We will also discuss some qualitative properties of solutions to the enthalpy formulation on an evolving surface.

This talk is based on joint work with Thomas Sales and Chandrasekhar Venkataraman (University of Sussex)

ABSTRACT. The discrete de Rham (DDR) complex is an arbitrary-order discretisation of the continuous de Rham complex that can be defined on general polytopal meshes. After rewriting Einstein's equations using exterior calculus, giving them a form reminiscent of Maxwell's equations, one can apply the DDR method to obtain a numerical scheme that preserves certain auxiliary constraints. The main momentum constraint's, however, remain elusive, and we check numerically their growth in some simulations.

ABSTRACT. Dispersive nonlinear partial differential equations can be used to describe a range of physical systems, from water waves to spin states in ferromagnetism. The numerical approximation of solutions with limited differentiability (low regularity) is crucial for simulating fascinating phenomena arising in these systems, including emerging structures in random wave fields and dynamics of domain wall states, but it poses a significant challenge to classical algorithms. Recent years have seen the development of tailored low-regularity integrators to address this challenge.

Inherited from their description of physical systems, many such dispersive nonlinear equations possess a rich geometric structure, such as a Hamiltonian formulation and conservation laws. To ensure that numerical schemes lead to meaningful results, it is vital to preserve this structure in numerical approximations. This, however, results in an interesting dichotomy: the rich theory of existent structure-preserving algorithms is typically limited to classical integrators that cannot reliably treat low-regularity phenomena, while most prior designs of low-regularity integrators break geometric structure in the equation. In this talk, we will outline recent advances incorporating structure-preserving properties into low-regularity integrators. Starting from simple discussions on the non-linear Schrödinger equation, we will showcase a unified framework that allows for the very first systematic construction of symmetric resonance-based integrators to approximate a wide class of nonlinear dispersive equations at low regularity. This unified framework builds on so-called decorated tree series and a novel formalism based on forest formulae. This is joint work with Yvonne Alama Bronsard, Yvain Bruned, and Katharina Schratz.

ABSTRACT. The use of machine learning for modeling, understanding, and controlling large-scale physics systems is quickly gaining in popularity, with examples ranging from electromagnetism over nuclear fusion reactors and magneto-hydrodynamics to fluid mechanics and climate modeling. These systems – governed by partial differential equations – present unique challenges regarding the large number of degrees of freedom and the complex dynamics over many scales both in space and time, and additional measures to improve accuracy and sample efficiency are highly desirable. We present an end-to-end equivariant surrogate model consisting of an equivariant convolutional autoencoder and an equivariant convolutional LSTM using G-steerable kernels. As a case study, we consider the three-dimensional Rayleigh-Bénard convection, which describes the buoyancy-driven fluid flow between a heated bottom and a cooled top plate. While the system is E(2)-equivariant in the horizontal plane, the boundary conditions break the translational equivariance in the vertical direction. Our architecture leverages vertically stacked layers of D₄-steerable kernels, with additional partial kernel sharing in the vertical direction for further efficiency improvement. We demonstrate significant gains in sample and parameter efficiency, as well as better scaling to more complex dynamics.

ABSTRACT. Many dynamical systems in physics and other fields possess some form of geometric structure, such as Lagrangian or Hamiltonian structure, symmetries, and conservation laws. Preserving these structures in the course of discretisation typically leads to numerical algorithms with improved stability properties and reduced errors for long-time and strongly nonlinear simulations.

Most geometric integrators known and used today are based on linear approximation spaces such as splines or finite elements. When approximating differential equations with a strongly nonlinear solution manifold, these methods tend to require a large number of degrees of freedom. The reason is simple: Even if the intrinsic dimension of the nonlinear solution space is very low, an accurate approximation by a linear space requires a high dimensional approximation space. It thus seems desirable to use nonlinear approximation spaces instead, potentially with drastically reduced dimension.

In this talk, we will explore the construction of geometric integrators based on nonlinear representations of the solution such as neural networks or symbolic expressions. We will address some of the challenges such as initialisation and solution of the resulting discrete Euler-Lagrange equations, and show first, promising results for ODEs.

ABSTRACT. In this talk, we address the inverse problem of simultaneously recovering free energy defined on the density manifold from discretized observations of the density flow, which is generated by a Wasserstein gradient or a Hamiltonian flow. We formulate the problem as an optimization task that minimizes a loss function specifically designed to enforce the underlying variational structure of Wasserstein flows, ensuring consistency with the geometric properties of the density manifold. Our framework employs a kernel-based operator approach within the associated Reproducing Kernel Hilbert Space (RKHS), yielding a closed-form representation of the unknown components. Furthermore, a comprehensive error analysis is conducted, providing convergence rates under adaptive regularization parameters as the temporal and spatial discretization mesh sizes tend to zero. Finally, a stability analysis is presented to bridge the gap between discrete trajectory data and continuous-time flow dynamics for the Wasserstein Hamiltonian flow.

ABSTRACT. Partial Differential Equations (PDEs) are fundamental tools for modeling phenomena across physics, biology, and engineering, yet their numerical approximation becomes particularly challenging in the presence of singularities, discontinuities, or complex geometries. While classical numerical methods provide strong theoretical guarantees, they typically rely on high regularity assumptions and encounter significant limitations in non-smooth or high-dimensional settings. Recent advances in scientific machine learning introduce new approximation paradigms that may offer promising ways to overcome these challenges, primarily by leveraging variational formulations to construct optimal nonlinear surrogate models.

In this talk, we present a novel framework based on variational formulations and their associated gradient flows to study the computability and complexity of continuous PDE solutions from an optimization perspective. Our approach approximates PDE solution operators via discrete gradient flows and links structural properties of PDEs—such as coercivity, ellipticity, and convexity—to the intrinsic complexity of their solutions. Furthermore, we outline possible extensions to discontinuous functions using Hybrid Surrogate Models and discuss the potential of analog hardware as a platform for developing new computability frameworks for non-smooth problems.

ABSTRACT. Sampling invariant distributions from an It\^o diffusion process presents a significant challenge in stochastic simulation. Traditional numerical solvers for stochastic differential equations require both a fine step size and a lengthy simulation period, resulting in biased and correlated samples. The current deep learning-based method solves the stationary Fokker--Planck equation to determine the invariant probability density function in the form of deep neural networks, but they generally do not directly address the problem of sampling from the computed density function. In this work, we introduce a framework that employs a weak generative sampler (WGS) to directly generate independent and identically distributed (iid) samples induced by a transformation map derived from the stationary Fokker--Planck equation. Our proposed loss function is based on the weak form of the Fokker--Planck equation, integrating normalizing flows to characterize the invariant distribution and facilitate sample generation from a base distribution. Our randomized test function circumvents the need for min-max optimization in the traditional weak formulation. Our method necessitates neither the computationally intensive calculation of the Jacobian determinant nor the invertibility of the transformation map. A crucial component of our framework is the adaptively chosen family of test functions in the form of Gaussian kernel functions with centers related to the generated data samples. Experimental results on several benchmark examples demonstrate the effectiveness and scalability of our method, which offers both low computational costs and excellent capability in exploring multiple metastable states. This is based on arXiv 2509.12841 (to appear in SISC 2026, joint work with Z. Cai, Y. Huang and Y. Cao)

ABSTRACT. Learning dynamical systems that respect physical symmetries and constraints remains a fundamental challenge in data-driven modeling. Integrating physical laws with graph neural networks facilitates principled modeling of complex N-body dynamics and yields accurate and permutation-invariant models. However, training graph neural networks with iterative, gradient-descent-based optimization algorithms (e.g., Adam, RMSProp, LBFGS) often leads to slow training, especially for large, complex systems. In comparison to 15 different optimizers, we demonstrate that Hamiltonian Graph Networks (HGN) can be trained 150-600x faster--but with comparable accuracy--by replacing iterative optimization with random feature-based parameter construction. We show robust performance in diverse simulations, including N-body mass-spring and molecular systems in up to 3 dimensions and 10,000 particles with different geometries, while retaining essential physical invariances with respect to permutation, rotation, and translation. We reveal that even when trained on minimal 8-node systems, the model can generalize in a zero-shot manner to systems as large as 4096 nodes without retraining. Our work challenges the dominance of iterative gradient-descent-based optimization algorithms for training neural network models for physical systems.

ABSTRACT. We review the Parker Problem, i.e., whether a braided magnetic field can relax under ideal plasma dynamics to a smooth force-free equilibrium. Several variations of this problem exist [1], but in this talk, we consider only the case of line-tied and periodic boundary conditions for magnetic braids. We revisit older results obtained with a Lagrangian code using mimetic operators [2,3] and explain the advantages and weaknesses of this approach. We contrast these results with more recent analytical results by Enciso and Peralta-Salas [4] for periodic boundary conditions.

[1] Pontin and Hornig, The Parker problem: existence of smooth force-free fields and coronal heating, Living Reviews in Solar Physics, Volume 17(5), 2020

[2] Candelaresi,Pontin, and Hornig, Mimetic Methods for Lagrangian Relaxation of Magnetic Fields, SIAM J. on Sci. Comp., 36(6) 2014.

[3] Candelaresi, Pontin and Hornig, Magnetic field relaxation and current sheets in an ideal plasma, Astrophysical Journal 808:134, 2015.

[4] Enciso and Peralta-Salas, Obstructions to Topological Relaxation for Generic Magnetic Fields, Arch. Rational Mech. Anal. 249:6 (2025)

ABSTRACT. The magneto-frictional equations are used in solar physics to compute both static and quasi-static models of the Sun’s coronal magnetic field. Here, we examine how accurately magneto-friction (without fluid pressure) is able to predict the relaxed state in a one dimensional test case containing two magnetic null points. Firstly, we show that relaxation under the full ideal magnetohydrodynamic equations in the presence of nulls leads necessarily to a non-force-free state, which could not be reached exactly by magneto-friction. Secondly, the magneto-frictional solutions are shown to lead to breakdown of magnetic flux conservation, whether or not the friction coefficient is scaled with magnetic field strength. When this coefficient is constant, flux is initially conserved, but only until discontinuous current sheets form at the null points. In the ensuing weak solution, we show that magnetic flux is dissipated at these current sheets. The breakdown of flux conservation does not occur for an alternative viscous relaxation scheme.

ABSTRACT. In this talk, I will present an enriched Galerkin (EG) method for the two-dimensional incompressible Navier–Stokes equations that preserves both kinetic energy and enstrophy in the inviscid limit. The method is based on an H(gradrot)-conforming EG velocity space enriched by the lowest-order Raviart–Thomas space, together with a piecewise-constant pressure space. Using the rotational form of the nonlinear term, we propose both a nonlinear scheme and a linearized variant. Both schemes exactly conserve the discrete kinetic energy and enstrophy, and the Picard iteration for the nonlinear scheme preserves these invariants at each iteration. I will also discuss the stability analysis and a priori error estimates for the nonlinear scheme. Several numerical experiments will be presented to demonstrate the accuracy and conservation properties of the proposed method.

ABSTRACT. Fusion energy has great potential as a long-term solution to society’s growing energy needs. However, realising fusion requires confining highly ionised gases (i.e. plasmas) at extreme temperatures and pressures. Through the combined interaction of long-range electromagnetic forces and short-range collisional interactions, the dynamics of fusion plasmas are extremely multi-scale and highly nonlinear, which makes modelling these systems a problem at the leading edge of high-performance scientific computing. In the continuum limit, fusion plasmas are described by variations of the magnetohydrodynamic (MHD) equations, ranging from ideal to visco-resistive and so-called extended-MHD. In this presentation, we provide an overview of recent progress in extended-MHD modelling for so-called tokamaks and stellarators, which are fusion energy systems based on toroidal confinement of high-temperature plasmas. Generically, these systems are extremely ill-conditioned and stiff and have long posed challenges for high-fidelity simulations of dynamical fusion phenomena including instabilities, relaxation and loss-of-control events. We will discuss recent advancements, the state-of-the-art and current challenges and bottlenecks experienced within the fusion community, with the view of finding opportunities for collaboration.

ABSTRACT. In this talk, we present a two-phase model for incompressible viscoelastic flow. The model is governed by the incompressible Navier-Stokes equations with an additional elastic stress tensor in the two fluid phases. We use the Oldroyd-B constitutive law to describe the viscoelastic behaviour in both phases. We first motivate the study by outlining recent advancements in the numerical analysis of viscoelastic fluids, specifically addressing critical challenges such as energy decay, positivity preservation, and convergence. Subsequently, we explain how to couple these equations with two distinct interface descriptions: the Cahn-Hilliard phase-field model and a sharp-interface approach utilizing linear parametric finite elements. We discuss recent progress for both frameworks, emphasizing techniques to ensure energy stability and volume conservation. Finally, we demonstrate the robustness of the proposed methods through numerical experiments, including convergence rate studies and typical benchmark problems.

ABSTRACT. Many important phenomena in material science, cell biology, and geometry can be formulated as curvature driven interface problems, including phenomena such as reshaping of cell membranes and melting and solidification processes. In this talk, we study cut finite element method (CutFEM)–based discretizations as a flexible computational strategy for evolving interface problems without requiring remeshing.

We review the key methodological ingredients of the approach, including level set evolution, reinitialization procedure, and the overall computational workflow for simulations on moving domains. The framework is designed to accommodate a spectrum of interface coupled multiphysics problems. Examples include two phase flow, incompressible flow with embedded rigid body motion, flow coupled to complex interface dynamics, and situations involving topological changes. The approach is implemented on a structured octree based background mesh within the Gridap ecosystem \and supports distributed memory parallel execution.

ABSTRACT. One of the main causes of degradation of solid insulators is electrical treeing: a typically thin, elongated, and highly branched gas-filled fracture formed by progressive erosion of the polymeric surface. This defect is caused by the prolonged action of intense electric fields that trigger Partial Discharges. The underlying physical problem can be modeled by a system of Partial Differential Equations [3], taking into account the movement of charges in the fracture and the evolution of the electric field and potential in both materials.

The geometrical complexity of the gas domain introduces a prohibitive computational cost for 3D simulations of this phenomenon, which we overcome with a dimensional reduction of the treeing geometry, approximated as a 1D graph embedded in the 3D dielectric domain, and the derivation of a mixed-dimensional model [1,2]. An accurate solution of the electrostatic equation and a proper representation of all the components of the electric field within the defect are fundamental to properly capture the discharge dynamics.

We propose a numerical solution of the electrostatic equation based on mixed Finite Elements (FEM) in 3D and FEM on the 1D graph, with a coarse 3D mesh, independent of the 1D discretization, resulting in a reduced number of degrees of freedom. The transverse components of the electric field in the 1D domain are then computed a posteriori as the superimposition of different contributions.

The reduced model is tested on simplified geometries to assess the performance and accuracy of the method, and then on a realistic tree geometry, achieving a significant reduction in the computational cost. Finally, considering independent meshes on the two domains allows simulations in contexts where the generation of a fully conforming 3D mesh would be unfeasible.

References: [1] B. Crippa, A. Scotti, and A. Villa, “A mixed-dimensional model for the electrostatic problem on coupled domains,” Journal of Computational Physics, p. 114 015, 2025. [2] B. Crippa, A. Scotti, and A. Villa, “A one-dimensional reduced plasma model for the electrical treeing,” arXiv preprint arXiv:2512.07900, 2025. [3] A. Villa, L. Barbieri, M. Gondola, A. R. Leon-Garzon, and R. Malgesini, “A PDE-based partial discharge simulator,” Journal of Computational Physics, vol. 345, pp. 687–705, 2017.

ABSTRACT. We present a fast and accurate method for evaluating layer and volume potentials for the modified Helmholtz equation in two dimensions. The method is based on a decomposition of the Green’s function into a short-range local part and a smooth long-range part. The long-range contribution is evaluated efficiently using the nonuniform fast Fourier transform (NUFFT), while the local contribution is further decomposed into an analytically treated singular component and a remainder represented by a telescoping sum over dyadic refinement levels. This decomposition enables high-order accurate evaluation of the local contribution while keeping the computational cost of the NUFFT low.

The volumetric discretization is constructed by embedding the computational domain in a uniform triangular background grid. Elements intersected by the boundary are cut and retriangulated to obtain a boundary-fitted cut–cell mesh with curved boundary elements represented through isoparametric mappings. Unlike finite element and finite volume methods, the present approach does not discretize differential operators on the mesh. Instead, the solution is represented entirely through integral operators, and the mesh is used solely for numerical quadrature. Consequently, the accuracy depends primarily on the quadrature rule rather than on element shape, making the method largely insensitive to degenerate, anisotropic, or highly skewed triangles. No stabilization or cell-merging procedures are therefore required.

We apply the method to interior Dirichlet and Neumann problems on complex geometries, including translated and rotated domains that produce cut cells with extreme aspect ratios. Numerical experiments demonstrate third-order convergence with respect to mesh size and show that the accuracy remains stable even in the presence of severely distorted elements. The results indicate that the proposed framework provides an efficient, geometrically flexible, and robust approach for high-order potential evaluation on complex domains.

ABSTRACT. Recently, continuous optimization methods based on the RSAV (relaxed scalar auxiliary variable) approach have been proposed, which ensure modified dissipation of gradient flows. Existing studies, however, rely solely on gradient information, and often exhibit degraded convergence when applied to ill-conditioned problems. In this work, we incorporate Hessian information into the linear operator arising in the RSAV scheme to accelerate convergence. Directly using the Hessian, however, presents two key challenges: (i) increased computational cost due to both the construction of the Hessian and the solution of the associated linear systems at each iteration; and (ii) potential loss of the modified dissipation property of the RSAV scheme when the Hessian is indefinite. To overcome these challenges, we propose using an approximate Hessian as the linear operator. To mitigate the first issue, we apply the Nyström method, a randomized low-rank approximation technique, to approximate the Hessian efficiently. To address the second, we enforce positive semidefiniteness of the approximate Hessian via eigenvalue decomposition. Numerical experiments on the training of PINNs (physics-informed neural networks), a class of typically ill-conditioned optimization problems, demonstrate that the proposed method yields accelerated convergence.

ABSTRACT. This paper deals with Physics-Informed Gaussian Processes (GPs), a probabilistic method to solve linear parametric Partial Differential Equations (PDEs). In fact, we use GP priors to infer statistical estimators to the PDE solution given possibly noisy observations. More specifically, this paper is devoted to bridge a gap in the literature and provide convergence rates of such estimators to the PDE solution in terms of Sobolev norms. Unlike previous works that consider using the PDE's Green function as covariance kernel for the GP, or that assume the stability of the differential operator and the solution map of the PDE between corresponding Sobolev spaces, we highlight sufficient assumptions on the target function, differential operator and covariance kernels to obtain such stability and convergence rates. In particular, we show that the Matérn kernels fall into these conditions. We illustrate these points with different numerical experiments.

ABSTRACT. Physics-Informed Machine Learning (PhiML) provides a structured framework for combining data-driven models with classical representations of physical systems, particularly in the study of dynamical systems and partial differential equations. In this session, we examine a set of established PhiML methodologies through illustrative examples, with an emphasis on modeling assumptions, mathematical structure, and practical trade-offs.

We begin with surrogate modeling via neural operators, focusing on Fourier Neural Operators (FNOs) and related architectures such as Transolver. These approaches are presented in the broader context of operator learning, where mappings between infinite-dimensional function spaces are approximated using finite parameterizations. We discuss the effectiveness of these surrogate modeling techniques for parametric PDEs and offer insights into when and how they are useful.

We then revisit Physics-Informed Neural Networks (PINNs) as a constrained learning framework in which differential equations are incorporated as soft penalties in the loss function. Rather than emphasizing novelty, we focus on the variational interpretation of PINNs, their expressivity, and well-known optimization and consistency challenges.

Throughout, MATLAB-based implementations are used to ground the discussion while emphasizing unifying theoretical ideas. The techniques presented — including FNOs, Transolver, PINNs — are now established tools within the PhiML landscape, and the goal of this talk is to clarify how they can be systematically understood and applied.

ABSTRACT. We propose a hybrid method for the steady \mbox{Darcy--Forchheimer} equation in which the coarse nonlinear solve of the two-level mixed finite element scheme is replaced by a physics-informed neural network (PINN) surrogate. The network is trained to minimize the residual of the governing equations, and its output is projected onto the coarse finite element space and used as the linearization point for the fine-grid correction. We introduce a computable dual residual norm~$\eta_{H}^{NN}$ that bounds the PINN coarse error and serves as a practical quality indicator before the fine solve is run. Under the mesh relation $h=O(H^{2})$, the hybrid scheme achieves first-order convergence in $h$, matching the rate of the classical two-level method. Numerical experiments confirm the theoretical rates at $\beta=1$ and reveal that pressure accuracy degrades at higher Forchheimer numbers, consistent with the growth of~$\eta_{H}^{NN}$.

ABSTRACT. Modern deep neural networks achieve strong performance across domains but often rely on highly overparameterised dense layers. This motivates the study of structured, low-rank operators that retain expressive power while controlling computational complexity. Tensor network formats, such as the tensor train (TT) decomposition, provide compressed representations of high-order linear maps, but their use in deep learning has largely been limited to linear layers or post-compression of pre-trained layers. We introduce the double haircomb (DH) operator, a nonlinear tensor network layer that can be composed into multi-layer architectures and trained directly in compressed form. The DH consists of two TTs coupled through a single connecting bond, forming a tree topology that can represent both multilinear and nonlinear maps between high-order tensor spaces through the incorporation of activation functions within the tensor network structure. We study both the computational and representational properties of this architecture. Empirically, small DH networks achieve competitive performance on hyperspectral image classification while using an order of magnitude fewer parameters than comparable dense models. We further characterise the architectural properties underlying this performance through theoretical analysis and synthetic validation experiments.

ABSTRACT. Despite decades of advances in computational science, predictive models of complex systems—from turbulence and climate to battery ageing—remain fundamentally incomplete. Physics‑based equations typically omit unresolved processes and external perturbations, while purely data‑driven models often lack interpretability and fail to generalize. Here we present a physics‑guided stochastic modelling framework that bridges these paradigms. Using Alternating Neural Integrators (ANI), we decompose the system dynamics into a deterministic component derived from domain knowledge and a stochastic residual inferred from observational data. The residual is represented by a conditional normalizing flow enhanced with sinh–arcsinh transformations, which captures non‑Gaussian features including skewness and heavy tails. We validate the framework across multiple domains, demonstrating that it recovers missing stochastic forcing and reproduces empirical trajectory statistics even when the deterministic model is deliberately simplified. As a practical application, we deploy the approach for lithium‑ion battery health prognostics using NASA’s ageing dataset. Conditioned on an equivalent‑circuit model, the learned stochastic compensator captures cycle‑to‑cycle variability and model mismatch during charge–discharge operation, yielding probabilistic forecasts that align with observed voltage evolution and capacity fade. Our results establish that imperfect physics models, when systematically augmented with data‑driven stochastic structure, yield predictions that are not only more accurate but also retain physical interpretability—offering a principled route toward reliable forecasting in complex systems where first‑principle knowledge is partial but indispensable.

ABSTRACT. Tensor decompositions are a fundamental tool in scientific computing and data analysis. In many applications — such as simulation data on irregular grids, surrogate modeling for parameterized PDEs, or spectroscopic measurements — the data has both discrete and continuous structure, and may only be observed at scattered sample points. The CP-HIFI (hybrid infinite-finite) decomposition generalizes the Canonical Polyadic (CP) tensor decomposition to settings where some factors are finite-dimensional vectors and others are functions drawn from infinite-dimensional spaces — a natural framework when the underlying data has continuous structure. The decomposition can be applied to a fully observed tensor (aligned) or, when only scattered observations are available, to a sparsely sampled tensor (unaligned). Current methods compute CP-HIFI factors by solving a sequence of dense linear systems arising from regularized least-squares problems to fit reproducing Kernel Hilbert space (RKHS) representations to the data, but these direct solves become computationally prohibitive as problem size grows. We propose new algorithms that achieve the same accuracy while being orders of magnitude faster. For aligned tensors, we exploit the Kronecker structure of the system to efficiently compute its eigendecomposition without ever forming the full system, reducing the solve to independent scalar equations. For unaligned tensors, we introduce a preconditioned conjugate gradient method applied, exploting the problem's structure for fast matrix-vector products and efficient preconditioning. We analyze the computational complexity and memory requirements of the new methods and demonstrate their effectiveness on problems with smooth functional modes. I will also discuss the “First Proof” project, which aims to understand the capabilities of AI systems on problems that come up in math research, and the role that results from that experiment played in this project.

ABSTRACT. Langevin-type and interacting-particle dynamics arise naturally in statistical mechanics and molecular dynamics via an SDE description. In the mean-field limit, the associated particle density evolves according to a Fokker-Planck or McKean-Vlasov equation. For the examples we are interested in, these equations can be viewed as Wasserstein gradient flows. These dynamics are also closely connected to sampling problems, since their equilibria encode the probability distributions of interest. When the associated free energy is non-convex, or when the linearised dynamics has a small spectral gap, convergence to equilibrium may be slow, and metastable behaviour can occur.

In this talk, we discuss a feedback-control framework obtained by linearising the dynamics around a target equilibrium and solving an algebraic Riccati equation to construct controlled potential perturbations acting on the slow eigenmodes of the linearised operator. The resulting closed-loop dynamics is exponentially stable at any prescribed rate. In addition, the feedback can be interpreted geometrically as a local convexification of the free energy via a finite-rank quadratic correction targeting slow or unstable modes. For Gibbs targets of the form exp(-V/sigma) that can be sampled from Langevin dynamics, the transformed linearised operator is self-adjoint in a suitable space, and both the operator and the Riccati feedback can be assembled using only the unnormalised weight exp(-V/sigma), bypassing the unknown partition function. This yields an explicit feedback law in the eigenbasis of the linearised operator. In this way, the original sampling problem is recast as a stabilisation problem, which in the Gibbs setting, reduces in practice to the numerical computation of some eigenvalues and eigenfunctions of the linearised operator.

At the particle level, the mean-field feedback induces a controlled interacting particle system implementable through empirical coordinates over the slow modes. We prove well-posedness of the closed-loop McKean-Vlasov SDE and establish finite-time propagation-of-chaos estimates with error bounds in Wasserstein distance. Therefore, the feedback control implemented at the PDE level can promote accelerated convergence to the target distribution.

ABSTRACT. Langevin dynamics is a diffusion process extensively used, in particular in molecular dynamics simulations, to sample Gibbs measures. Some alternatives based on piecewise deterministic kinetic velocity jump processes have gained interest over the last decade. One interest of the latter is the possibility to split forces at the continuous-time level, reducing the numerical cost for sampling the trajectory. Motivated by this, we introduce a numerical scheme based on hybrid dynamics combining velocity jumps and Langevin diffusion, numerically more efficient than their classical Langevin counterparts. Beyond computational speedup and thermodynamical sampling, dynamical properties such as diffusion coefficients are also well preserved. Furthermore, by leveraging the random nature of jumps, our framework can be integrated to classical multi-time-step (MTS) integrators, which helps mitigating the resonance issues of MTS. Finally, we will discuss the extension of this framework to "distilled" multi-time-step methods, recently developed for molecular dynamics using neural network potentials.

ABSTRACT. Computing the probability of rare events is a challenging task. Indeed, in that case, performing direct molecular dynamics simulations is often computationally prohibitive. Hence, we use the Adaptive Multilevel Splitting (AMS) algorithm. However, AMS struggles to converge when the initial velocity distribution corresponds to low values of the committor (i.e. the probability to observe the rare event): this results in a high variance. Therefore, we introduce an importance sampling method on the initial velocities, which allows to reduce the variance of the probability estimator. The method is first tested on a 1D double well potential before being applied to a 2D potential and a real life example: the alanine dipeptide.

Joint work with Pierre Marmey, Laura J. S. Lopes, Tony Lelièvre, Gabriel Stoltz, Hadrien Vroylandt

ABSTRACT. In this talk, we introduce a novel kinetic Langevin sampler based on a specific splitting scheme using the exact harmonic Langevin integrator. For strongly log-concave target measures the sampler makes use of a splitting in a quadratic and a convex $L_G$-$\nabla$Lipschitz perturbation of the strongly convex potential. For this splitting we consider both exact gradients and stochastic gradient approximations and provide convergence rates in $L^2$ Wasserstein distance and a detailed complexity analysis. In particular we obtain the same order for the contraction rate as for the continuous dynamics. To achieve $\varepsilon$-accuracy the required order of the steps size is comparable to that of the standard splitting schemes as OBABO or UBU.

ABSTRACT. We consider the nonlinear Schrödinger equation with competing cubic-quintic power law nonlinearities in two space dimensions. This equation is globally well-posed and admits a rich variety of standing, solitary wave solutions. In particular, when considered on wave-guide domains, there exist line solitary waves whose transverse (in-)stability is numerically shown to depend on the existence of a critical length in the transverse direction.

ABSTRACT. The numerical approximation of low-regularity solutions to the nonlinear Schrodinger equation (NLSE) is notoriously difficult and even more so if structure-preserving schemes are sought. Recent works have been successful in establishing symmetric low-regularity integrators for NLSE. However, so far, all prior symmetric low-regularity algorithms are fully implicit, and therefore require the solution of a nonlinear equation at each time step, leading to significant numerical cost in the iteration. In this work, we introduce the first fully explicit (multi-step) symmetric low-regularity integrators for NLSE. We demonstrate the construction of an entire class of such schemes which notably can be used to symmetrise (in explicit form) a large amount of existing low-regularity integrators. We provide rigorous convergence analysis of our schemes and numerical examples demonstrating both the favourable structure preservation properties obtained with our novel schemes, and the significant reduction in computational cost over implicit methods.

ABSTRACT. In this talk, we present several numerical methods for solving the nonlinear Schrödinger equation (NLSE) within a low-regularity framework, including logarithmic and semi-smooth nonlinearities. Compared to the cubic NLSE, low-regularity models exhibit distinct dynamical behaviors and may involve singular or semi-smooth nonlinear terms, posing significant challenges for both analysis and computation. In particular, the limited regularity of the nonlinearities may lead to order reduction in classical error estimates.

To address these issues, we develop some techniques and analyze a class of numerical schemes, including time-splitting Fourier spectral methods, and exponential wave integrator (EWI) Fourier spectral methods, which are well suited to low-regularity settings. We then establish rigorous error bounds that avoid order degradation and ensure the accuracy of the proposed methods. This is joint work with Prof Li-Lian Wang(NTU,Singapore).

ABSTRACT. We introduce a novel approach to the numerical approximation of the nonlinear Schrodinger equation with white noise dispersion in the regime of low-regularity solutions. Approximating such solutions in the stochastic setting is particularly challenging due to randomized frequency interactions and presents a compelling challenge for the construction of tailored schemes. In particular, we design the first resonance-based schemes for this equation, which achieve provable convergence for solutions of much lower regularity than previously required. A crucial ingredient in this construction is the Wong–Zakai approximation of a stochastic dispersive system, which introduces piecewise-linear phases that capture nonlinear frequency interactions and can subsequently be used to construct resonance-based schemes. We prove the well-posedness of the Wong–Zakai approximated equation and establish its proximity to the original full stochastic dispersive system. Based on this approximation, we demonstrate an improved strong convergence rate for our new scheme, which exploits the stochastic nature of the dispersive terms. Finally, we present numerical experiments that underscore the favorable performance of our novel method in practice. This is a joint work with Dr. Georg Maierhofer.

ABSTRACT. We present efficient iterative methods for computing low-rank solutions to nonlinear matrix and tensor equations arising from discretizations of partial differential equations. The central challenge is handling nonlinearity while maintaining sublinear computational cost through compressed low-rank representations.

For matrix equations G(X) = X, we introduce low-rank Anderson acceleration (lrAA), which adapts the classical Anderson acceleration framework to the SVD format. A key algorithmic component is Cross-DEIM, an adaptive cross approximation that uses DEIM-based index selection with warm-starting to efficiently evaluate nonlinear functions in low-rank form. An adaptive scheduling strategy controls the truncation tolerance to keep intermediate ranks low throughout the iteration.

We then extend these ideas to Tucker tensor format for solving nonlinear tensor equations in moderate dimensions, introducing Cross²-DEIM, a fiber-sampling Tucker cross approximation guided by DEIM, and Tucker-AA, an Anderson acceleration method operating directly in Tucker format. Numerical experiments on benchmark PDE problems in 2D and 3D demonstrate the efficiency and robustness of both methods.

ABSTRACT. Neural networks have delivered impressive performance across many tasks, but their memory footprint and computational demands can limit deployment on resource-constrained devices. Existing compression strategies often rely on sparsification, identifying subnetworks that nearly preserve baseline accuracy, while recent fine-tuning methods leverage low-rank updates for efficient adaptation of pretrained models. In this talk, we introduce an in-training low-rank compression strategy for neural networks with tensor-structured weights. Building on dynamical low-rank approximation techniques, our method adaptively updates the rank in each tensor mode during training, yielding an automatic and layer-dependent compression. We demonstrate the effectiveness and robustness of the approach through numerical experiments on a wide range of benchmarks, highlighting the favorable accuracy–compression tradeoff.

ABSTRACT. Federated Learning (FL) is a distributed Machine Learning (ML) paradigm that enables multiple local devices, that is, clients, and a central server to collaboratively train an ML model using data stored locally on the clients without transferring the data. Common challenges in FL include data privacy preservation, communication efficiency, and convergence guarantees. This presentation introduces FeDLRT: Federated Dynamical Low-Rank Training, which adopts a dynamical low-rank training scheme that constrains model training to low-rank factors of the model parameters, resulting in reductions in both communication and memory costs during the training process. By incorporating a variance-reduction scheme, convergence results are established for FeDLRT, even when the data distribution is heterogeneous, that is, non–independent and identically distributed, among the local clients. The advantages of FeDLRT are demonstrated through numerical experiments on training neural networks with various architectures in an FL setting.

ABSTRACT. The aim of this talk is to present novel global space-time methods for the approximation of the time-dependent Schrödinger equation on low-complexity manifolds. The backbone of the approach is the use of a least-square formulation of the time-dependent Schrödinger equation with many-body Coulomb interaction potentials, which can be obtained using Kato theory. The latter can be used in conjunction with low-rank tensor formats (such as Tensor Trains for instance) gaussian wavepackets to derive new variational principles to compute dynamical low-complexity approximations of the solution. These new approximations are different from the ones obtained through the classical Dirac-Frenkel principle. One significant advantage of this new variational formulation is that the existence of a dynamical low-rank approximation for any finite-time horizon can be proved with low-rank tensor formats, whereas dynamical low-rank approximations constructed with the Dirac-Frenkel principle can usually be proved to exist only locally in time. Illustrative numerical results will be presented to highlight the differences between the dynamical low-complexity approximations obtained with these different approaches.

ABSTRACT. In solid state physics, the electronic properties of an L-periodic crystal are often described in terms of an L-translation invariant energy functional defined on a set of L-periodic, locally trace class operators. As a consequence of Bloch's theorem, electronic quantities of interest can be expressed as integrals over the Brillouin zone of so-called energy bands, which are piecewise smooth, Lipschitz continuous periodic functions obtained by solving a parametrised family of periodic, non-linear elliptic eigenvalue problems. The practical computation of these quantities of interest therefore typically involves the use of uniform trapezoidal quadrature in the Brillouin zone together with spectral plane-wave discretisations of the parameterised family of non-linear eigenvalue problems.

In this talk, I will discuss the numerical analysis of such discretisations in the case of the finite-temperature Kohn-Sham LDA pseudo-potential energy functional. I will also briefly discuss the case of the Hartree-Fock energy functional for insulating crystals.

This presentation is based on joint-work with Eric Cancès (on Kohn-Sham DFT) and Alfred Kirsch (on Hartree-Fock).

ABSTRACT. In this talk, I discuss the convergence in Zhidkov spaces of the first-order Lie-Trotter and the second-order Strang splitting schemes for the time integration of the Gross-Pitaesvkii equation with a time-dependent potential and non-zero boundary conditions at infinity. In particular, after adapting the Lie derivatives framework for the numerical analysis of splitting methods to this setting, we can show the conservation of the generalized mass and the near-preservation of the Ginzburg-Landau energy balance law. Finally, numerical accuracy tests performed on a one-dimensional dark soliton corroborate our theoretical findings and we investigate the nucleation of quantum vortices in two experimentally relevant settings. This is joint work with Quentin Chauleur (Inria and Université de Lille).

ABSTRACT. We develop structure-preserving time integration schemes for Gaussian wave packet dynamics associated with the magnetic Schrödinger equation. The variational Dirac--Frenkel formulation yields a finite-dimensional Hamiltonian system for the wave packet parameters, where the presence of a magnetic vector potential leads to a non-separable structure and a modified symplectic geometry. By introducing kinetic momenta through a minimal substitution, we reformulate the averaged dynamics as a Poisson system that closely parallels the classical equations of charged particle motion. This representation enables the construction of Boris-type integrators adapted to the variational setting. In addition, we propose explicit high-order symplectic schemes based on splitting methods and partitioned Runge--Kutta integrators. The proposed methods conserve the quadratic invariants characterizing the Hagedorn parametrization, preserve linear and angular momentum under symmetry assumptions, and exhibit near-conservation of the averaged Hamiltonian over long time intervals. Rigorous error estimates are derived for both the wave packet parameters and observable quantities, with bounds uniform in the semiclassical parameter. Numerical experiments demonstrate the favorable long-time behavior and structure preservation of the integrators.

ABSTRACT. Generative models, including variational autoencoders, normalizing flows, generative adversarial networks, and diffusion models, have dramatically advanced the realism and quality of generated images, text, and audio. Beyond these tasks, generative models hold great promise as powerful tools for probability density estimation and high-dimensional sampling, which are central to uncertainty quantification (UQ) tasks such as amortized Bayesian inference and data assimilation. However, while research on image synthesis emphasizes producing high-quality individual samples, UQ applications require accurate approximation of statistical quantities of interest rather than visually realistic samples. As a result, direct application of existing generative models to UQ problems can lead to biased approximations or unstable training. In this talk, we will introduce several new generative approaches tailored to UQ. These include training-free diffusion models for density estimation, a score-based nonlinear filter for data assimilation, and training-free conditional diffusion models for amortized Bayesian inference. We will demonstrate their effectiveness across a range of tasks, including density estimation for unimodal and multimodal distributions, learning stochastic dynamical systems, parameter estimation via amortized inference, and scalable data assimilation for atmospheric models.

ABSTRACT. Neural operators are neural network-based surrogate models for approximating solution operators of parametric partial differential equations, enabling efficient many-query computations in science and engineering. Many applications, including engineering design, involve variable and often nonparametric geometries, for which generalization to unseen geometries remains a central practical challenge. In this work, we adopt a kernel integral perspective motivated by classical boundary integral formulations and recast operator learning on variable geometries as the approximation of geometry-dependent kernel operators, potentially with singularities. This perspective clarifies a mechanism for geometric generalization and reveals a direct connection between operator learning and fast kernel summation methods. Leveraging this connection, we propose a multiscale neural operator inspired by Ewald summation for learning and efficiently evaluating unknown kernel integrals, and we provide theoretical accuracy guarantees for the resulting approximation. Numerical experiments demonstrate robust generalization across diverse geometries for several commonly used kernels and for a large-scale three-dimensional fluid dynamics example.

ABSTRACT. The Brinkman penalization method is a technique that simplifies the numerical calculation of wave phenomena in complex domains. It avoids the need for complex meshing by adding a damping term to the equations in subdomains where obstacles exist. Because of this term, the resulting equations are conservative in some regions and dissipative in others. In this talk, we will explain a geometric property of these equations and discuss their applications in scientific machine learning.

ABSTRACT. Neural operators provide fast surrogate solvers for parametric partial differential equations (PDEs), but purely data-driven training often requires large datasets and can be brittle in small-data and out-of-distribution test settings. In this talk, I will present the Physics-Informed Laplace Neural Operator (PILNO), a physics-based extension of the Laplace Neural Operator (LNO). PILNO incorporates governing physics directly into training by enforcing PDE, boundary condition, and initial condition losses, improving data efficiency and robustness for time-dependent problems. To further strengthen generalization, we introduce an additional set of unlabeled “virtual” inputs that broaden the diversity of physics-only supervision, and we apply a temporal-causality weighting that prioritizes early-time dynamics to stabilize optimization. Across representative benchmarks including Burgers’ equation, Darcy flow, reaction-diffusion system, and a forced KdV equation, PILNO consistently improves accuracy in small-data regimes and shows better generalization beyond the training distribution.

ABSTRACT. In this talk we focus on Bayesian inverse problems in which the forward parameter-to-observable map is approximated in a stochastic way, for instance by Monte Carlo (MC) simulations. Such problems arise for example in uncertainty quantification of particle transport, where the parameter-to-observable map is defined through the solution of the Boltzmann equation. Approximating the likelihood with high accuracy requires MC simulations with many samples, which makes sampling from the posterior distribution expensive. We present an efficient method for sampling from the posterior based on pseudo-marginal sequential Monte Carlo (SMC) using likelihood tempering. To accelerate sampling from the posterior the method adopts a multilevel approach. A significant speedup is achieved compared to a single-level SMC using high-fidelity MC simulations while targeting the same posterior, which is demonstrated by numerical experiments.

ABSTRACT. Real-world traffic systems are inherently uncertain, due to heterogeneous driver behavior, environmental conditions, and infrastructure variability. Kinetic traffic models provide a flexible mesoscopic framework to describe vehicle interactions and incorporate uncertain parameters while capturing complex collective dynamics. To propagate these uncertainties, non-intrusive uncertainty quantification techniques are particularly suitable, as they allow the use of existing deterministic solvers for the kinetic equations. In this talk, we adopt a Monte Carlo framework to study the propagation of uncertainty, performing sampling in the stochastic space while solving the kinetic equation in the physical space. However, the resulting high-dimensional problem leads to significant computational costs when standard Monte Carlo estimators are employed. To mitigate this issue, we investigate control variate strategies, specifically multi-level and multi-fidelity Monte Carlo methods, which proved to be particularly effective. The latter exploits a hierarchy of models, where high-fidelity simulations provide accurate but computationally demanding solutions, while low-fidelity approximations offer computational efficiency at the cost of reduced precision. Combining these models, we achieve a significant reduction in computational cost while maintaining high accuracy. Numerical simulations indicate that these approaches provide substantial accuracy improvements over standard Monte Carlo methods. Moreover, by using appropriate low-fidelity surrogates, multi-fidelity methods can outperform multilevel Monte Carlo methods.

ABSTRACT. We present a numerical framework combining Telescopic Projective Integration (TPI) with Stochastic Galerkin polynomial chaos (SG-gPC) for the spatially inhomogeneous Boltzmann equation with uncertainty. The equation has a two-scale structure governed by the Knudsen number ε: fast collisional relaxation and slow hydrodynamic transport. Explicit time integration of this system requires resolving the fast scale, making simulations in the transitional regime computationally expensive. TPI exploits the spectral gap by performing a few inner time steps to damp the fast transient, then extrapolating the slow dynamics via projective Runge–Kutta methods at a larger time step. The stochastic uncertainty is handled through an SG-gPC expansion. The collision integral is evaluated via a fast spectral method, while spatial transport uses a MUSCL finite-volume scheme. We discuss the construction of the coupled TPI–SG hierarchy for the nonlinear Boltzmann equation and investigate accuracy, efficiency, and convergence properties through numerical experiments on benchmark problems.

ABSTRACT. We consider the modeling of neutral particles in the scrap-off-layer, i.e., outer region of tokamak fusion reactors. Here, one generally models plasma as a fluid, coupled with a kinetic equation modeling neutral particles. For reactor design, these equations are part of a set of PDE constraints for the optimization problem to be solved. Solving this problem requires both the ability to simulate the PDE, as well as to compute the corresponding gradient of the objective function to the design parameters.

These kinetic equations are typically simulated using particle-based Monte Carlo methods, which introduce noise to the simulation outputs and, consequently, the objective function being minimized. Care must then be taken when computing the gradient of the objective function to its design parameters, to ensure convergence. The discrete-adjoint method computes the gradient through a simulation that retraces the original particle trajectories backwards in time. For complex simulations, however, storing these trajectories requires a large amount of memory.

In this work, we tackle this memory issue by recomputing reversed trajectories from scratch. To this end, we use a reversible random number generator, i.e., one that can step both forwards and backwards through its sequence at identical cost. We demonstrate that this approach significantly reduces the memory required in a straightforward implementation of a simplified 1D coupled plasma-neutral model.

ABSTRACT. Neural networks are well known for their remarkable flexibility in approximating continuous functions, yet their capabilities reach far beyond this classical setting. In this talk, we turn to functional neural networks—an emerging and promising approach for approximating nonlinear smooth functionals or operators. By investigating the convergence rates of both approximation and generalization errors, we uncover key theoretical properties of these networks within the empirical risk minimization framework. This analysis not only deepens our understanding of functional neural networks but also sets the stage for their effective applications in functional data analysis and operator learning.

ABSTRACT. This talk addresses unbounded density ratio estimation--an underexplored yet critical challenge in statistical learning--and its application to covariate shift adaptation. We propose a three-step procedure utilizing unlabeled data from both source and target distributions: (1) estimating a relative density ratio; (2) truncating to control unboundedness; and (3) transforming the estimate back to the standard density ratio, which is then used as importance weights for regression. We establish non-asymptotic convergence guarantees for both the density ratio estimator and the resulting regression estimator, showing that under mild conditions, both achieve optimal or near-optimal rates. This work provides new theoretical insights into density ratio estimation and learning under covariate shift, extending classical theory to more practical settings.

ABSTRACT. We studied the problem of learning the Green's functions of partial differential equations from data, through reproducing kernel methods. With the help of a novel kernel design, we derived an algorithm of time complexity $O(m^3+m^2N)$ only, where $N$ was the size of training sample, and $m$ was the number of grid points. Minimax lower bound and upper bound of learning rates were derived. Numerical examples on elliptic equations demonstrated accurate approximation of Green's functions.

ABSTRACT. Numerical simulations of evolving surfaces, such as mean curvature flow and shape optimization, often suffer from mesh deterioration under large deformations. In this talk, I will present a numerical framework based on the Minimal Deformation Rate (MDR) strategy, whose main purpose is to preserve mesh quality throughout the evolution without remeshing.

First, for curvature-driven flows, including mean curvature flow and surface diffusion on both closed and open surfaces with moving contact lines, we design a new BGN--MDR scheme by exploiting the key observation that the Barrett–Garcke–Nürnberg (BGN) and MDR formulations differ only in one degree of freedom, thereby combining the discrete energy stability of the classical BGN approach with the mesh-quality preservation of MDR.

Second, for PDE-constrained shape optimization and surface hole filling, we couple a second-order inertial flow with MDR, which both accelerates convergence in flat energy landscapes and avoids frequent remeshing during large geometric deformations through a suitable MDR-based mesh-motion strategy.

ABSTRACT. A proof of optimal-order H^1-norm error estimates is presented for a backward difference full discretization of order 1 to 5 of Willmore flow for closed two-dimensional surfaces. The analysed numerical method discretizes a coupled system of evolution equations by evolving surface finite elements and backward difference method of order 1 to 5 in time. The presentation focuses on the novel techniques developed for the stability part of the convergence analysis, which is based on energy estimates exploiting the anti-symmetric structure of the second-order system, in combination with Dahlquist's G-stability and the multiplier technique of Nevanlinna and Odeh, with a new upper bound in the spirit of Dahlquist. Numerical experiments illustrate and complement the theoretical results.

ABSTRACT. The lattice-Boltzmann method (LBM) has become a popular tool for computational modeling of fluid flows in recent decades due primarily to its explicit time stepping and parallel scalability. However, in its standard form, the lattice-Boltzmann method is confined to regular, Cartesian grids. This geometrical limitation can prove to be a large impediment to accurate flow simulations. Indeed, when simulating contact-line dynamics, the restriction to Cartesian grids often leads to non-physical pinning of contact lines.

Conversely, while the finite-element method (FEM) is a versatile, geometry-agnostic tool for computational modeling, performing 3D simulations of fluid flows in FEM has been difficult due to the nonlinearity of the governing equations and the prohibitive cost of solving large linear systems.

Here we propose a new method wherein we simulate contact-line dynamics by solving the lattice-Boltzmann equations using the finite-element method. Our method preserves some of the parallel scalability and the explicit time stepping of LBM while inheriting the geometrical flexibility of FEM. We showcase our new method by simulating the spreading of a droplet on a patterned surface, and simulating a droplet bouncing on a spherical particle.

ABSTRACT. Numerical simulations play a central role in understanding complex dynamical systems, from fluid motion to large-scale physical models. In many such systems, important quantities—such as energy or spectral invariants—are preserved over time, and capturing these features accurately is essential for reliable long-term simulations.

In this talk, we consider isospectral Lie–Poisson systems on matrix Lie algebras, which arise naturally in areas such as mechanics and fluid dynamics. We present a new framework for analyzing structure-preserving numerical methods for these systems, focusing on isospectral symplectic Runge–Kutta methods.

Our approach builds on classical backward error analysis and is formulated in terms of Butcher series, an algebraic framework for analyzing numerical integrators. A key ingredient is a Lie–Poisson reduction of Butcher series. More precisely, we show that, in the isospectral Lie–Poisson setting, the modified equations arising from backward error analysis retain a Lie–Poisson structure. This provides a geometric interpretation of how the numerical method reflects the underlying dynamics. In particular, this perspective explains the good long-time preservation of invariants: the numerical method can be interpreted as exactly solving a nearby modified equation that remains Lie–Poisson and isospectral, and therefore preserves the same class of invariants.

Beyond its practical implications, this framework highlights connections between numerical analysis, geometry, and algebra, and offers new tools for the systematic study of structure-preserving algorithms. We illustrate these ideas with examples from fluid dynamics, emphasizing how geometric insight can guide the design of robust numerical methods.

This is joint work with Klas Modin.

ABSTRACT. In many applications, such as stochastic optimization, molecular dynamics, or statistical learning, one is interested in the long-term behavior of ergodic stochastic systems evolving under geometric constraints. Standard approaches often rely on penalization methods, embeddings in high-dimensional spaces, and extrinsic discretizations, which leads to high computational costs and severe restrictions on the timestep. We consider a general class of ergodic stochastic differential equations evolving on Riemannian manifolds, including in particular the Riemannian Langevin dynamics, and we develop numerical methods aimed at sampling their invariant measure with high order of accuracy. As the analysis of the order of convergence involves technical Taylor expansions, we introduce a framework based on Butcher series that provides a systematic way to describe the integrations by parts in the calculations. Numerical experiments confirm the theoretical results. This work is in collaboration with Adrien Busnot Laurent.

ABSTRACT. Aromatic Butcher series are the natural extension of the standard B-series that allows to compute the divergence of a Taylor expansion over the jet space of a given ODE. First used as a tool for the study of volume-preserving numerical integrators, they were soon studied for their numerous algebraic properties. In this talk, we will present an overview of recent results, both numerical and algebraic, concerning aromatic trees. After an intuitive definition of aromas with an invariant tensor theorem, we introduce the aromatic bicomplex, in the spirit of the variational bicomplex in variational calculus. This new object proves crucial for characterising Ker(Div) on aromatic trees, derive negative results on volume-preserving methods, and formulate the notion of symmetries and Noether theorems for tree-like structures. If time allows, we will give an overview of related algebraic structures, that are aromatic multi-indices, Hopf algebroids, and planar aromatic trees, and discuss how they will be used to tackle the open problems of geometric integration concerning volume-preservation.

ABSTRACT. Wave kinetic theory is an actively growing area of research concerned with deriving effective macroscopic wave equations from microscopic dynamics. For the nonlinear Schrödinger equation in dimensions d≥2, a rigorous mathematical framework for wave turbulence has been successfully established. In contrast, the one-dimensional setting remains much less understood, both physically and mathematically.

As in the higher-dimensional case, one may expand the solution into iterates, encoded by paired decorated trees, and reduce the proof to a counting problem on these pairings. However, in one dimension, this counting problem is significantly more difficult, leading to the divergence of some of these Feynman diagrams on time scales much shorter than the expected kinetic time.

In this talk, we show how this combinatorial difficulty can be overcome for the Benjamin–Ono (BO) equation by exploiting its integrable structure, and in particular a recent explicit solution formula. By suitably rescaling and iterating this explicit formula, combined with probabilistic arguments, we obtain insight on the Wave Kinetic theory for BO, up to the physically relevant kinetic time scale. As a consequence, we also obtain a numerical scheme for efficiently simulating these effective dynamics.

ABSTRACT. Langevin diffusions on Euclidean space are widely used to generate samples from given Gibbs distributions. Exploiting their underlying Hamiltonian structure, these diffusions can be extended to Lie groups of compact type, which are reductive Lie groups whose semisimple part is compact, since these allow for a compatible Riemannian structure via the existence of bi-invariant metrics. This leads to the momentum Langevin diffusion, which serves as a natural analogue of the kinetic Langevin diffusion on Euclidean space. Splitting methods provide a natural class of structure-preserving integrators for these equations. In this talk, we discuss this momentum Langevin equation and several splitting schemes for its numerical integration, and we show how the performance of these schemes can be analyzed in a natural way using the Baker–Campbell–Hausdorff (BCH) formula.

ABSTRACT. Weak convergence of the Euler approximation of SDEs on Riemannian manifolds is well understood, and has been studied in both intrinsic coordinates and on embedded submanifolds of Euclidean space. Despite this, the mean-square convergence rate of the Euler scheme is still unknown. Moreover, a rigorous derivation of the scheme has yet to be derived.

In this talk, I will derive the manifold Euler method and show that it has (on non-flat manifolds) global convergence rate of order 1/2 when the manifold is Cartan-Hadamard. I will then go further and derive a Milstein method whose mean-square global convergence rate is of order 1.

This is ongoing work with Karthik Bharath and Michael Tretyakov.

ABSTRACT. Lie-Poisson systems appear in many areas of physics. They exhibit a strong geometric structure, as they evolve on coadjoint orbits and have conserved quantities. By adding noise, it is possible to account for uncertainties and unresolved smaller scales. An important question is how to numerically integrate these systems in a structure-preserving manner. In this talk, we describe an integrator for a class of stochastic Lie-Poisson systems driven by Stratonovich noise. We describe how its derivation follows from discrete Lie-Poisson reduction of the symplectic midpoint scheme for stochastic Hamiltonian systems and discuss some of the properties of the scheme.

ABSTRACT. In molecular dynamics, nonabelian lattice gauge theories, and the constrained sampling problems ubiquitous in machine learning, one often needs to solve stochastic differential equations on curved spaces. Specifically, the efficient numerical integration of Langevin equations is required. By combining symplectic geometry and the theory of homogeneous spaces, we obtain an elegant theoretical framework that guides the numerical discretisation. For separable problems, we employ splitting integrators, and for nonseparable problems, we utilize variational integrators on groupoids.

ABSTRACT. When the nonconvex problem is complicated by stochasticity, the sample complexity of stochastic first-order methods may depend linearly on the problem dimension, which is undesirable for large-scale problems. To alleviate this linear dependence, we adopt non-Euclidean settings and propose two choices of non-smooth proximal terms when taking the stochastic gradient steps. This approach leads to stronger convergence metric, incremental computational overhead, and potentially dimension-insensitive sample complexity. We also consider further acceleration through variance reduction which achieves near optimal sample complexity and, to our best knowledge, is the first such result in the $\ell_1$/$\ell_\infty$ setting. Since the use of non-smooth proximal terms is unconventional, the convergence analysis deviates much from algorithms in Euclidean settings or employing Bregman divergence, providing tools for analyzing other non-Euclidean choices of distance functions. Efficient resolution of the subproblems in various scenarios is also discussed and simulated. We illustrate the dimension-insensitive property of the proposed methods via preliminary numerical experiments.

ABSTRACT. The Physics-Informed Neural Network (PINN) method has demonstrated significant effectiveness in solving partial differential equations (PDEs). However, a limited amount of rigorous convergence analysis supports this approach. Consequently, conducting analytical research on the convergence properties of the PINN method has become increasingly important. In this work, we focus on elliptic equations and analyze the L-infinity norm and L-p norm of the loss terms. We then introduce the concept of stability for PDEs in a continuous sense. Finally, numerical experiments are conducted to validate the theoretical estimations.

ABSTRACT. This talk will present a brain-inspired, spiking, neuromorphic algorithm for solving sparse linear systems (Ax = b) such as those arising in finite element methods for partial differential equations (PDEs), one of the most important techniques in modern numerical science and engineering. The algorithm embeds the sparse matrix into the synaptic connections between subpopulations of neurons directly without training or learning. Neural dynamics are defined such that the collective spiking activity of the whole network flows to an efficient spiking representation of the solution vector x, with comparable numerical accuracy to traditional algorithms. We demonstrate this algorithm on real neuromorphic hardware (Intel’s Loihi 2) and show close to ideal strong and weak scaling. We demonstrate the generality of the algorithm through several examples of PDEs in 2 and 3 dimensions, with nontrivial mesh topologies, and different boundary conditions. Our work establishes a direct connection between established numerical methods for PDEs and brain-like spiking neural networks, demonstrating the value of brain inspiration, and expanding the neuromorphic footprint in scientific computing.

ABSTRACT. Modeling three-dimensional magnetohydrostatic equilibria in stellarators is of paramount importance in the design and conduction of fusion experiments. Several workhorse methods of the community simplify the search for equilibria to configurations with nested flux surfaces. It is known, however, that relevant configurations in experiments feature magnetic islands and chaotic regions, violating this assumption. Magnetic relaxation codes provide an opportunity to solve for more general equilibrium solutions. We introduce a structure-preserving magnetic-relaxation solver (MRX) that does not require a nested flux surface assumption. The novelty of MRX is that it combines several crucial features for the first time. (i) Through structure-preserving mixed finite-elements, we retain central features of the relaxation problem in the discrete setting, such as div B = 0 to machine precision, guaranteed energy-dissipation and helicity preservation. (ii) The use of B-Spline mixed-finite elements allows us to use high-order basis functions that exhibit rapid convergence for regular solutions as well as non-uniform meshes. This flexibility also allows us to use the flux surface mapping from codes like VMEC in our computations (iii) The code is pure Python using the JAX ecosystem, making it very easy to install, run, and extend as well as performant on GPUs. It is also fully differentiable for future inverse design and optimization applications. The proposed magnetic relaxation solver is tested in several stellarator geometries at low and high values of β. Future work will address the integration of this code for 3D equilibrium optimization in stellarator fusion devices.

ABSTRACT. In this paper, we propose a high-order Lagrangian finite element method for ideal magnetic relaxation that preserves both topological structure and energy dissipation. In the Lagrangian setting, the electromagnetic fields are computed via pullbacks induced by the discrete flow map, thereby exactly preserving the topology of magnetic field lines. To achieve energy dissipation, we reformulate the discrete relaxation model as a Sobolev gradient flow that minimizes the discrete magnetic energy. The key insight is to identify the first variation of the discrete magnetic energy as a weak Lorentz force term derived from the Maxwell stress tensor. Combined with average vector field time integration, we construct a fully discrete scheme that preserves energy dissipation. The mesh and the electromagnetic fields are discretized by high-order Lagrange finite element spaces and finite element exterior calculus (FEEC) spaces, respectively, achieving arbitrarily high-order accuracy and improving tolerance under strong mesh distortion. With matrix-free techniques and sum factorization, the scheme can be implemented efficiently on graphics processing units (GPUs) for long-time, large-scale simulations. Numerical experiments demonstrate the accuracy, energy dissipation, topology preservation, and robustness to severe mesh distortion of the proposed numerical scheme.

ABSTRACT. We investigate the long-term dynamics of the advection-diffusion of finite element differential forms. Specifically, we prove that under certain conditions, the solutions remain bounded, thereby excluding dynamo action. These stationary solutions can be characterized as the "harmonic forms" associated with the advection-diffusion operators. In the absence of advection, classical Hodge theory establishes that the dimension of the space of harmonic forms coincides with the Betti numbers of the underlying manifold. In this work, we establish a discrete analogue of a result by V.I. Arnold: the number of linearly independent stationary solutions is no less than the corresponding Betti number. Furthermore, for the case of potential flows, we employ Witten transforms to provide a deeper characterization of these discrete forms. This is joint work with Jindong Wang.

ABSTRACT. Magnetic relaxation drives plasma toward lower-energy equilibria under helicity constraints. In ideal MHD, helicity is locally conserved; in resistive settings like Taylor relaxation, only global helicity is preserved, permitting local reconnection. We present two structure-preserving finite element methods based on finite element exterior calculus and Lagrange multipliers, which enforce helicity conservation at these two different levels. Applied to the magneto-frictional system (gradient flow of MHD), our schemes preserve a discrete Arnold inequality and maintain nontrivial magnetic topology over long time. Numerical experiments on braided and knotted fields confirm that local helicity preservation prevents spurious reconnection, while global-only conservation allows further relaxation. We then apply the relaxation idea to MHD turbulence model to preserve ideal invariants and capture the selective decay mechanism, which is believed to be a mechanism leading to the formation of large magnetic eddies.

ABSTRACT. Oscillatory Graph Neural Networks (OGNNs) are an emerging class of physics-inspired architectures designed to mitigate oversmoothing and vanishing gradient problems in deep GNNs. In this work, we introduce the Complex-Valued Stuart-Landau Graph Neural Network (SLGNN), a novel architecture grounded in Stuart-Landau oscillator dynamics. Stuart-Landau oscillators are canonical models of limit-cycle behavior near Hopf bifurcations, which are fundamental to synchronization theory and are widely used in e.g. neuroscience for mesoscopic brain modeling. Unlike harmonic oscillators and phase-only Kuramoto models, Stuart-Landau oscillators retain both amplitude and phase dynamics, enabling rich phenomena such as amplitude regulation and multistable synchronization. The proposed SLGNN generalizes existing phase-centric Kuramoto-based OGNNs by allowing node feature amplitudes to evolve dynamically according to Stuart-Landau dynamics, with explicit tunable hyperparameters (such as the Hopf-parameter and the coupling strength) providing additional control over the interplay between feature amplitudes and network structure. We conduct extensive experiments across node classification, graph classification, and graph regression tasks, demonstrating that SLGNN outperforms existing OGNNs and establishes a novel, expressive, and theoretically grounded framework for deep oscillatory architectures on graphs.

ABSTRACT. Operator learning is a rapidly growing field that aims to approximate nonlinear operators related to partial differential equations (PDEs) using neural operators. These rely on discretization of input and output functions and are, usually, expensive to train for large-scale problems at high-resolution. Motivated by this, we present a Multi-Level Monte Carlo (MLMC) approach to train neural operators by leveraging a hierarchy of resolutions of function dicretization. Our framework relies on using gradient corrections from fewer samples of fine-resolution data to decrease the computational cost of training while maintaining a high-level accuracy. The proposed MLMC training procedure can be applied to any architecture accepting multi-resolution data. Our numerical experiments, on a range of state-of-the-art models and test-cases, demonstrate improved computational efficiency compared to traditional single-resolution training approaches, and highlight the existence of a Pareto curve between accuracy and computational time, related to the number of samples per resolution.

ABSTRACT. Neural networks (NNs) may be viewed as discretizations of underlying continuous dynamical systems, both at the level of the network architecture and in the gradient flow used for parameter optimization. In both contexts, stability properties of the discretization methods can be relevant, for instance in enhancing adversarial robustness. To address challenges such as exploding or vanishing gradients, neural network feature and/or parameter spaces are often modeled as Riemannian manifolds. In this work, we propose a general framework for analyzing the stability of numerical integrators on Riemannian manifolds. As a concrete application, we analyze the explicit Euler method in the Riemannian setting—referred to as the Geodesic Explicit Euler (GEE) method—and establish precise stepsize stability bounds when the manifold has constant sectional curvature. Finally, we focus on Hyberolic spaces and provide non-expansive architectures defined on such spaces.

ABSTRACT. This talk introduces PDE-based Group Equivariant Convolutional Neural Networks (PDE-G-CNNs), offering a continuous, dynamical systems perspective on deep learning. We replace traditional discrete convolution, pooling, and activation layers with a geometric Hamilton-Jacobi evolution PDE solver.

Key highlights of this framework include: * Intrinsic Symmetries: Formulating the PDEs on homogeneous spaces guarantees built-in symmetries (such as roto-translation equivariance), completely eliminating the need for data augmentation. * Interpretable Parameters: Instead of learning opaque kernel weights, the network learns geometrically meaningful PDE coefficients. * Operator Splitting: We map PDE terms to specific network functions: convection handles data transport, fractional diffusion provides regularization, and morphological convolutions naturally subsume non-linearities like ReLUs and max/min-pooling.

Ultimately, framing deep learning as a continuous geometric evolution process allows PDE-G-CNNs to achieve highly competitive performance on imaging tasks while requiring only a fraction of the parameters used by traditional neural networks.

ABSTRACT. Variable surface-tension two-phase flows couple the evolution of an interface, the transport of surface-active species, geometric quantities such as curvature, and the incompressible flow. This coupling creates numerical challenges related to conservation, accuracy, efficiency, and stability. In this talk, I will discuss a space–time cut finite element approach for such problems. The starting point is a conservative high-order unfitted method for convection–diffusion equations in evolving domains, in which Reynolds’ transport theorem is used to ensure global mass conservation at the discrete level. The method combines a space-time formulation with quadrature in time and ghost-penalty stabilization. I will also discuss how the transport discretization can be coupled to unfitted approximations of curvature, surface-tension forces, and two-phase flow, with examples motivated by surfactant and active-surface dynamics.

ABSTRACT. Additive manufacturing techniques like 3D printing can produce complex structures and topologies. We introduce a shape and topology optimization problem in a phase field setting to obtain rigid structures capable of withstanding forces acting in the parts’ intended environment while ensuring constructability with additive manufacturing methods. To solve this problem numerically, the VMPT (Variable Metric Projection Type) method is applied in function space. A convergence result is stated. Convergence in function space indicates that the iteration numbers of the VMPT method are bounded, independent of the discretization level and the construction layers used. Numerical evidence of this is illustrated. Moreover, multiple solver specifications are discussed, which speed up the computation. The impact of the model parameters on the shape and topologies is shown. Furthermore, we present results in 2D and 3D, as well as with three materials, and for various manufacturing problem settings.

ABSTRACT. Today the classical homogenised models for simulating excitable tissue are challenged by new mathematical frameworks that explicitly represent and resolve the geometry of extracellular and intracellular spaces and cellular membranes (EMI models). These mixed-dimensional models crucially enable an abstract representation of heterogeneous distribution of membrane ion channels and realistic cellular morphologies. EMI models typically predict the electrical properties of the tissue, with the underlying assumption that the ion concentrations are constant. Although this assumption is only an approximation, the resulting models still give accurate predictions of neuronal electrodynamics in many scenarios. They do however fail in capturing the numerous phenomena related to shifts in the extracellular ion concentrations.

Here, we discuss an alternative approach to geometrically detailed modelling of excitable tissue that accounts for ion concentrations dynamics and electrodiffusion (KNP-EMI models). We introduce and numerically evaluate two new, finite element-based schemes for the KNP-EMI model, capable of efficiently and flexibly handling geometries of arbitrary dimension and arbitrary polynomial degree. The first scheme is based on a multi-dimensional formulation discretized with mortar finite elements, and the second is based on a single-dimensional formulation discretized with DG finite elements. We discuss stability and convergence properties of the two numerical schemes.

ABSTRACT. The Stokes mobility and resistance problems for rigid particles in viscous fluids arise in a wide range of natural and industrial settings. A persistent difficulty is the simulation of particles in near-contact, where discretizations become severely ill-conditioned and lubrication scales demand excessive global resolution. To address both issues, we introduce a two-body preconditioning strategy in a general framework suitable for any boundary value problem solver, and applicable more broadly to elliptic PDEs. The solution is expressed in a coarse one-body basis, obtained by solving one small boundary value problem per particle. For each close pair, this basis is locally corrected by a fine two-body solve that resolves the near-contact interaction. These corrections are precomputed, compressed to equivalent coarse representations, and applied within a global iterative solve that retains only coarse degrees of freedom. We demonstrate the approach using the method of fundamental solutions for two-dimensional Stokes flow, and illustrate it on suspensions of up to 10,000 particles. The resulting preconditioner delivers high accuracy in close-to-touching configurations, achieves linear scaling in the number of particles, and yields fast GMRES convergence.