ABSTRACT. Kernel methods, in the form of radial basis function and Gaussian process interpolation, have proved remarkably successful as a tool for various tasks in approximation and inference. This talk will focus on presenting recent advances in the design and numerical analysis of kernel methods in the context of uncertainty quantification in partial differential equations, which typically involves very high dimensional parameter spaces as well as physical constraints. To combat the challenge of high dimensions, we introduce anisotropic kernels and length-scale informed sparse grids that allow for efficient computations and accurate reconstructions also in this setting. To incorporate physical information, we present non-stationary kernels that are able to faithfully capture the structure and build this into the approximation.

ABSTRACT. Sessile liquid droplets on solid surfaces are ubiquitously found in nature and practical applications. They exhibit many intriguing properties and phenomena because their dynamics involve three phases interacting with each other, i.e., liquid, gas, and solid. Many physical processes participate in droplet dynamics, presenting a complex problem for fundamental understanding and practical applications.

Out of the many physical processes involved in droplet dynamics, two have been of continuous interest. The first is the moving contact line at which the evolving liquid-gas interface intersects the solid surface. As the three phases meet at the contact line, its motion involves the microscopic dynamics in the three-phase region coupled to the macroscopic hydrodynamics of the droplet. The second is the evaporation at the liquid-gas interface for evaporating sessile droplets. Coupled to both the moving contact line on the substrate and the liquid flow in the droplet, interfacial evaporation plays a critical role in controlling droplet dynamics.

Based on our earlier works on the modelling of moving contact line and thin-film dynamics, we derive a continuum model for evaporating droplets by applying Onsager’s variational principle. This approach ensures thermodynamical consistency in describing the coupling of many dissipative processes, including viscous momentum transport, contact line motion, evaporation, and vapor diffusion. Numerical results are presented to exhibit the diffusion-limited regime and the transition-limited regime, which are distinct from each other. The coffee-ring phenomena are also numerically investigated for droplets of colloidal suspension (coffee or tea).

This work is supported by the Hong Kong RGC General Research Fund (grant No. 16302525) and the Key Project of the National Natural Science Foundation of China (No. 12131010).

ABSTRACT. In this work, we develop a domain embedding model for Cahn-Hilliard-Navier-Stokes system with generalized Navier-slip boundary condition and dynamic boundary condition on complex domains based on Onsager's variational principle. The model enjoys the properties of thermodynamical consistency and asymptotic convergence to its limit. And the model is easily adapted to the numerical computation of two-phase flows with moving contact lines on an arbitrary domain. Numerical evidence on contact angle hysteresis on rough surfaces has been shown to demonstrate the model's computational efficacy.

ABSTRACT. Migrasomes are membrane-bound organelles whose mechanics remain poorly understood. Here, we develop a phase-field model to reveal how Tspan4 regulates migrasome behavior via membrane tension redistribution. We show that growth-induced membrane stretching generates tension gradients and stiffness heterogeneity, driving vertical elongation and a strong correlation between geometry and stiffness. By varying Tspan4 expression, we demonstrate that increased Tspan4 enhances bending rigidity and accelerates tension redistribution, stabilizing migrasomes, whereas Tspan4 depletion leads to central tension accumulation and higher rupture susceptibility. These findings identify membrane tension redistribution as a key mechanism linking Tspan4 expression to migrasome mechanics and stability.

ABSTRACT. In this talk, we present the dynamics of Hele-Shaw flow in a multi-connected region. We focus on our study on multiple interfaces in a radial Hele-Shaw cell. To compute the long-time interface morphologies and the motion of interfaces, we develop a spectral-accurate boundary integral method and a non-stiff time updating scheme. Numerical results reveal that the classic Saffman-Taylor instability is heavily influenced by the initial configuration of the interfaces and the physical parameters, e.g. the viscosity of the fluids. In particular, compared to the single-interface scenario, multiple interfaces can be used to promote or inhibit the development of fingering instability, thus facilitating the idea of morphological control in practice. Motivated by the self-similar theory and shape control strategy, one can preselect a limiting morphology, promote the growth of desired symmetry by a multi-interface setup, and realize it with a much smaller interface size by suppressing the unfavorable Saffman-Taylor instability.

ABSTRACT. In this talk I will review work on the analysis of motion capturing data and similar applications using techniques of shape analysis and deep learning. I will then consider a method for learning the Lagrangian and forces for mechanical systems using the discrete Lagrange d'Alembert principle. The case of manifold valued data and data on Lie groups will also be discussed if time permits. Applications to mechanical system will be considered.

ABSTRACT. We introduce Symplectic Neural Operators (SNOs) for learning the discrete-time flow of infinite-dimensional Hamiltonian systems arising from Hamiltonian PDEs. Our approach models the phase space as a Hilbert space endowed with a symplectic form induced by a skew-adjoint operator and constructs neural operators whose architectures are symplectic by design, ensuring structure preservation at the operator level. We provide a theoretical characterization of the symplecticity of the proposed operators and show how they approximate Hamiltonian flows in infinite-dimensional settings. Numerical experiments on canonical Hamiltonian PDEs demonstrate that SNOs achieve improved long-term stability and energy behavior compared to non-structure-preserving neural operators.

ABSTRACT. System inference for nonlinear ODE-driven dynamics becomes substantially harder when (i) only a subset of state components is observed, or (ii) measurements arrive as Repeated Cross-Sectional (RCS) snapshots (different individuals/units at each time) rather than full trajectories. To address this setting, we propose Simulation-based Generative Model for Imperfect Data (SiGMoID), an inference framework designed for partially observable system components and RCS data. SiGMoID combines two complementary ideas: (1) a physics-informed hyper-network (HyperPINN) that acts as a stable differentiable emulator of the system dynamics and (2) a Wasserstein GAN that matches the data distributions generated by the HyperPINN, allowing robust parameter estimation. SiGMoID quantifies observation noise, identifies ODE parameters, and reconstructs unobserved state components consistent with both the governing dynamics or the observed cross-sectional statistics. We validate the method on realistic experimental settings, demonstrating accurate recovery of full system dynamics under severe partial observability and cross-sectional sampling—supporting broad applications in biology, medicine, and engineered systems where longitudinal tracking is limited.

ABSTRACT. Stochastic gradient Langevin dynamics and its variants approximate the likelihood of an entire dataset, via random (and typically much smaller) subsets, in the setting of Bayesian sampling. Due to the (often substantial) improvement of the computational efficiency, they have been widely used in large-scale machine learning applications. It has been demonstrated that the so-called covariance-controlled adaptive Langevin (CCAdL) thermostat, which incorporates an additional term involving the covariance matrix of the noisy force, outperforms popular alternative methods. A moving average is used in CCAdL to estimate the covariance matrix of the noisy force, in which case the covariance matrix will converge to a constant matrix in long-time limit. Moreover, it appears in our numerical experiments that the use of a moving average could reduce the stability of the numerical integrators, thereby limiting the largest usable stepsize. In this article, we propose a modified CCAdL (i.e., mCCAdL) thermostat that uses the scaling part of the scaling and squaring method together with a truncated Taylor series approximation to the exponential to numerically approximate the exact solution to the subsystem involving the additional term proposed in CCAdL. We also propose a symmetric splitting method for mCCAdL, instead of an Euler-type discretisation used in the original CCAdL thermostat. We demonstrate in our numerical experiments that the newly proposed mCCAdL thermostat achieves a substantial improvement in the numerical stability over the original CCAdL thermostat, while significantly outperforming popular alternative stochastic gradient methods in terms of the numerical accuracy for large-scale machine learning applications.

ABSTRACT. Chromatin, the DNA-protein material which provides the building block of chromosomes, undergoes striking active reorganisations during the cell cycle. Replication constitutes one of the most striking such transformation. During S-phase, molecular motors such as helicases and DNA polymerases unwind the parental strands and synthesise new DNA, while chromatin is reshaped to accommodate the replication machinery. Textbooks often describe replication as a 1D process with the machinery tracking along the DNA template, but it actually occurs within a highly organised 3D nuclear environment, and the spatial dynamics remain incompletely understood. I will describe a 3D coarse-grained model for chromatin replication, and show how this naturally explains the spontaneous emergence of replication factories within a suitable range of parameters. These approaches provide mechanistic insight and lead to the predictions of a number of dynamical regimes which it would be interesting to test experimentally in the future. I will also outline a more general view of the genome as an active polymer maintained far from equilibrium, and some open problems within this context.

ABSTRACT. I will talk about our recent work on disordered and fluctuating active solids, including the active percolation transition described in Ref. [1]. This modelling effort involves a combination of deterministic and stochastic particle-based simulations, as well as their coarse-graining into continuum PDEs that are solved numerically. I will describe a variety of surprising phenomena that we uncover through our models, compare these models to experimental data, and emphasize the challenges of modelling elasticity far from equilibrium.

[1] More is Less in Unpercolated Active Solids. Jack Binysh, Guido Baardink, Jonas Veenstra, Corentin Coulais, and Anton Souslov. Phys. Rev. X 16, 021012 (2026)

ABSTRACT. For equilibrium hard spheres the statistical geometry of the insertion space, the room to accommodate another sphere, relates exactly to the equation of state. We begin to extend this idea to active matter, analyzing the insertion space for repulsive active particles in one and two dimensions using both on- and off-lattice models. In one-dimension we derive closed-form expressions for the mean insertion cavity size, cavity number, and total insertion volume, all in excellent agreement with numerical Langevin simulations. In two dimensions we develop a simple and efficient computational algorithm to track insertion cavities, which passes analytical tests. Strikingly, activity increases the total insertion volume and tends to keep the insertion space more connected. We also find that insertion space metrics contain signatures of collective phase behaviors occurring at previously predicted packing fractions. Taken together, our work provides the first quantitative foundation for the statistical geometry of active matter.

ABSTRACT. Inverse problems in differential equations are central to many scientific and engineering applications, requiring the estimation of model parameters based on noisy or incomplete observations. Traditional numerical methods for solving these problems are computationally expensive, especially when evaluating the likelihood in a Bayesian approach that involves high-dimensional parameter spaces and complex models. In this work, we investigate neural networks as surrogates to address this challenge. By incorporating a Laplace Approximation into the neural network in a Bayesian framework, our method efficiently approximates the forward model and provides calibrated uncertainty estimates. Compared to traditional methods, this approach significantly reduces computational costs while maintaining accurate posterior approximations. These findings underscore the potential of neural networks for scalable and reliable solutions to inverse problems in complex systems.

ABSTRACT. Many problems in uncertainty quantification require working with probability distributions that are often too complex to be treated directly. This creates a need for efficient numerical representations. An important class of approximation methods constructs mixture distributions, i.e. weighted sums of simple components such as Gaussian or Student’s t distributions. In this work, we present a unifying framework for greedy mixture methods that iteratively add new components to a mixture approximation in order to progressively improve the representation of the target distribution. This perspective brings together a range of seemingly distinct approaches—including iterated Laplace approximation, incremental and adaptive multiple importance sampling (IMIS/AMIS), variational boosting, and kernel-based methods such as kernel herding—within a common algorithmic framework. A general convergence result for the unified scheme is presented, and suitable assumptions so that it applies to existing methods are found; moreover, systematic numerical comparisons across different targets and performance criteria are performed, with particular emphasis on importance sampling and related sampling tasks, where mixture approximations yield efficient proposals.

ABSTRACT. In the geosciences, it is often necessary to interrogate the subsurface of the Earth for information about its structure and characteristics. A good way to do so can be to use recordings of either naturally or artificially induced seismic waves, since these have passed through the subsurface and are altered by subsurface (an)elastic properties. Spatially 3-dimensional, seismic full waveform inversion (FWI) is a computationally demanding inverse problem that constructs estimates of 3D subsurface seismic velocity structures using seismic waveform data. Bayesian FWI estimates the family of all such structures that are consistent with any available geological prior information, and with the observed data. This talk demonstrates that 3D Bayesian FWI can now be used to image the subsurface, analyse different prior geological hypotheses, and aid decision making about subsurface operations. We analyse a variety of prior hypotheses to reveal the sensitivity of the inversion process to different assumptions about the subsurface, using a method of variational prior replacement. This demonstrates that fully probabilistic, 3D Bayesian FWI can be performed, and can be used to test different prior hypotheses, at a cost that may be practical for certain applications.

ABSTRACT. The parareal algorithm combines a coarse and a fine time integrator to parallelize time-stepping schemes in time. When the material coefficient of a parabolic evolution problem depends on time and varies on a small time scale, the standard parareal algorithm produces substantially larger errors than in the classical setting. To ensure convergence with respect to the coarse time step size, the small time scale must be resolved. This talk addresses the use of averaged coarse integrators as a means of overcoming this difficulty. We present error bounds that extend the standard parareal error analysis, even for the classical heat equation, and that apply to the present setting. All bounds are explicit and require no regularity assumptions on the solution. Numerical experiments illustrate the theoretical results and demonstrate the effectiveness of the proposed approach.

ABSTRACT. This work studies the optimal choice of overlap in optimized Schwarz methods for the Poisson problem when the cost of the subdomain solves is taken into account. While increasing the overlap may accelerate convergence, it also enlarges the subdomain problems and thus can increase the cost of each iteration. We therefore introduce an efficiency function that balances error reduction and iteration cost under a general complexity model for the subdomain solver. We then analyze this function with respect to the overlap, establishing monotonicity and extremum properties and identifying regimes in which the optimal overlap is either the maximal admissible one or a strictly interior value. Numerical experiments support the theory and illustrate the influence of geometry and solver complexity on the optimal overlap.

ABSTRACT. We consider dynamic coupled problems such as heat transfer and fluid structure interaction, or more specifically PDEs that interact through a lower dimensional interface. Our general goal is a partitioned method that is high order, allows for different and adaptive time steps in the separate models, makes efficient use of hardware resources, is robust, and contains fast inner solvers. A prime candidate to fulfill this wishlist are waveform relaxation methods.

In order to design and effectively use such methods, it is important to have accurate error estimates. To this end, we analyze Dirichlet-Neumann waveform relaxation for two coupled 1D heterogeneous linear heat equations in both continuous and time discrete settings. Using a weighted Fourier transform, we derive new superlinear error bounds in the L2 norm in both settings. These bounds depend on the material parameters, domain sizes, the simulation time T as well as the time-step size. We further compare the continuous and time-discrete estimates and analyze how the time discretization influences the convergence behavior. Finally, we present numerical experiments that illustrate the accuracy of both the continuous and time-discrete error bounds.

ABSTRACT. The $p$-Laplace problem seeks to minimize a nonlinear functional of the form $J(u) = \int_{\Omega} \|\nabla u\|_2^p + fu \, dx$, where $\Omega \subset \mathbb{R}^d$ is a domain, $p \geq 1$ is a parameter, $f$ is a given forcing, and the unknown $u \in W^{1,p}(\Omega)$ is subject to some Dirichlet boundary conditions. When $p=2$, this coincides with the usual linear Laplace problem, expressed in energy minimization form (i.e. the Dirichlet principle). When $p \neq 2$, this is a nonlinear problem with applications, e.g. in compressed sensing. The objective functional $J(u)$ is convex in the unknown $u$ so the solution of the $p$-Laplace problem can be obtained by methods of convex programming. The advantage of this approach is that it produces an algorithm whose global convergence can be analyzed. In this talk, we shall discuss some versions of the barrier method for solving the $p$-Laplace problem. The basic barrier method converges in $\tilde{O}(\sqrt{n})$ Newton iterations, where the tilde indicates we neglect some polylogarithm. We shall also discuss newer multigrid-based algorithms that converge in $\tilde{O}(1)$ Newton iterations, for finite elements and for spectral elements, provided some regularity conditions are satisfied.

ABSTRACT. We develop a sparse spectral method for solving partial differential equations on a class of two-dimensional geometries bounded by algebraic curves. The numerical method uses generalised bivariate Koornwinder polynomials which form a complete orthogonal basis, but one which is not graded in terms of polynomial degree. The polynomials are built from new families of univariate semiclassical orthogonal polynomials whose associated operator matrices (Jacobi matrices, raising matrices and differentiation matrices) are computed with linear complexity in the number of basis functions. When used to discretise partial differential equations, the resulting matrices are sparse, enabling efficient numerical solutions. Moreover, we develop fast transforms from values on a grid to expansion coefficients. The efficiency and accuracy of the resulting spectral method are illustrated through a series of numerical experiments on geometries whose boundaries are smooth and piecewise smooth including non-convex geometries. We observe algebraic convergence for geometries with corners, which accelerates to exponentially fast (spectral) convergence when the boundary is smooth.

ABSTRACT. Full-potential Kohn-Sham density functional theory have been playing an increasingly important role in cutting-edge attosecond physics. In this talk, based on traditional self-consistent field iteration framework, a class of tetrahedral spectral element methods will be introduced for the spatial discretization, including the design and organization of basis functions, as well as solvers for the derived linear system. Notably, an h-adaptive mesh method is proposed for the framework to further resolve the singularity from the full-potential model. Numerical experiments show the spectral convergence of numerical solutions of model problems, and potential of the methods is clearly demonstrated by delivering ground state solutions of electronic structures towards the chemical accuracy.

ABSTRACT. We propose a stabilised auxiliary variable approach and construct high-order time-discrete schemes for a large class of gradient flows. By introducing an auxiliary variable a linearised time-discrete scheme can be obtained to achieve the energy stability, however, its discrete energy may deviate from the one in the original system. We define a specific auxiliary variable r=0 and construct an auxiliary equation for the auxiliary variable that stabilise the computation of the auxiliary variable and achieve a long-term good approximation between the discrete and original systems. We design second-order and higher-order time-discrete schemes and show the unconditional energy stability and in many cases the original energy law can be directly recovered. We also establish the optimal error estimate for the designed second-order scheme without any restriction on the time step. Ample numerical experiments including comparisons with some existing auxiliary variable methods will be presented as well. This is a joint work with Dr Zhaoyang Wang.

ABSTRACT. We consider the time-harmonic Maxwell equation in heterogeneous media with coefficients that are piecewise analytic and the interfaces across which the coefficients are allowed to jump are analytic. The equations are considered at large wavenumbers k and discretized by high order edge elements. We show that quasi-optimality is achieved if (a) the approximation order p is selected as p = O(log k) in conjunction with (b) a mesh size h such that kh/p is small. As in the related case of the Helmholtz equation, the analysis relies on a k-explicit regularity theory for Maxwell’s equations that decomposes solutions into two components: the first component is a piecewise analytic, but highly oscillatory function while the second one has finite regularity but features wavenumber-independent bounds.

ABSTRACT. The Serre-Green-Naghdi (SGN) equations are fluid models that incorporate dispersion (through higher order asymptotic expansions) into the traditional shallow water equation (SWE). While the presence of dispersion results in the modeling of important physical effects, the SGN equations contain a differential PDE constraint which poses a key numerical challenge.

This talk will focus on developing time integration and preconditioning strategies for handling the differential PDE constraint. First, we establish a constant coefficient preconditioner and rigorous bounds on the preconditioned conditioning number. The conditioning bounds incorporate the effects of bathymetry in two dimensions, are quasi-optimal within a class of constant coefficient operators, highlight fundamental scalings for a loss of conditioning, and ensure mesh independent performance for iterative Krylov methods.

Utilizing the conditioning bounds, we devise and test two time integration strategies for solving the full SGN equations. The first class combines classical explicit time integration schemes (4th order Runge-Kutta and 2nd--4th order Adams-Bashforth) with the new preconditioner. The second is a linearly implicit scheme where the differential constraint is split into a constant coefficient implicit part and remaining (stiff) explicit part. The linearly implicit methods require a single linear solve of a constant coefficient operator at each time step. We provide a host of computational experiments that validate the robustness of the preconditioners, as well as full solutions of the SGN equations including solitary waves traveling over an underwater shelf (in 1d) and a circular bump (in 2d).

ABSTRACT. Stiff systems of differential equations arise naturally in the modeling of complex physical phenomena, posing significant challenges for time integration. Exponential integrators are particularly well suited to this setting thanks to their favorable stability properties in the stiff regime. In this talk, we examine, from both theoretical and practical perspectives, the so-called directional split exponential integrators, a class of methods tailored to partitioned systems with d-dimensional Kronecker sum structure. On the computational side, we show in particular how these schemes can be efficiently implemented on modern computer architectures such as multi-core CPUs and GPUs, as they rely on high-performance level 3 BLAS. Numerical experiments demonstrate strong gains in efficiency and scalability compared to state-of-the-art methods in a variety of situations, for instance for diffusion--reaction systems exhibiting Turing pattern formation.

ABSTRACT. Large-scale systems of ordinary differential equations (ODEs) often exhibit multiple time scales, and many applications require understanding how their solutions depend on model parameters, for example through gradient computation. This talk presents forward and adjoint sensitivity analysis for MRI-GARK, a family of multirate methods that advances different ODE partitions with different time steps. This structure allows method design and adjoint checkpointing strategies to be tailored to each partition. Numerical experiments demonstrate improved performance and reduced checkpointing storage compared with single-rate methods. We also show how the framework can incorporate machine learning and lower-fidelity models into traditional sensitivity analysis workflows.

ABSTRACT. Modern gyrokinetic simulations of fusion plasmas, such as those performed with the GENE-X code, face significant computational bottlenecks due to the multi-scale nature of the governing equations. The system’s stiffness arises from the interplay between the advective terms of the Vlasov operator—including self-consistent electromagnetic effects—and the large diffusive terms associated with collisions and neutrals operators.

In this work, we present the first application of the Partitioned Runge–Kutta–Chebyshev (PIROCK) method to gyrokinetic simulations, implemented in the global, full-f GENE-X code. Unlike traditional operator-splitting techniques, PIROCK treats the full system within a unified framework, employing the explicit stabilized method ROCK2 for stiff diffusion while handling the advective Vlasov-electromagnetic terms with standard RK3. We briefly discuss the transition from a baseline PIROCK implementation using RK3 for the Vlasov part to one using RK4. Importantly, replacing RK3 with RK4 can, in some cases, provide up to 22% acceleration in total simulation time, a gain driven by the superior stability properties of RK4 for advective operators.

Furthermore, we explore the implementation of adaptive PIROCK schemes, which utilize a spectral radius estimator to dynamically adjust the timestep and the number of stages in response to transient physical phenomena. Preliminary tests indicate that this approach significantly improves efficiency compared to existing schemes (classical splitting and standard RK4). We will also discuss ongoing optimization efforts and outline possible future directions for enhancing large-scale gyrokinetic simulations.

ABSTRACT. On a two-dimensional Riemannian disk, there is a unique differential one-form whose exterior derivative equals the curvature and whose exterior coderivative vanishes, subject to appropriate boundary conditions. This one-form encodes the Levi-Civita connection with respect to a certain orthonormal frame. This talk will discuss a finite element method for approximating this one-form in settings where the Riemannian metric is only known approximately. We assume that the approximate metric belongs to the Regge finite element space, and we prove a convergence result for the computed connection one-form under the assumption that the approximate metric approaches the smooth metric under refinement. Our analysis relies on techniques that have been developed in recent years for studying distributional curvature, variational crimes, and Riemannian generalizations of trace inequalities and inverse inequalities.

ABSTRACT. Runge-Kutta methods are affine equivariant: applying a method before or after an affine change of variables yields the same numerical trajectory. However, for some applications, one would like to perform numerical integration after a quadratic change of variables. For example, in Lie-Poisson reduction, a quadratic transformation reduces the number of variables in a Hamiltonian system, yielding a more efficient representation of the dynamics. Unfortunately, directly applying a symplectic Runge-Kutta method to the reduced system generally does not preserve its Hamiltonian structure, so many proposed techniques require computing numerical trajectories of the original, unreduced system.

In this talk, we characterize when a Runge-Kutta method in the original variables descends to a numerical integrator expressible entirely in terms of the quadratically transformed variables. In particular, we show that symplectic diagonally implicit Runge-Kutta (SyDIRK) methods, applied to a quadratic projectable vector field, are precisely the Runge-Kutta methods that descend to a method (generally not of Runge-Kutta type) in the projected variables. We illustrate our results with several examples in both conservative and non-conservative dynamics.

ABSTRACT. We study a class of Hamiltonian partial differential equations introduced by Bridges (2006) in space and space-time. These systems satisfy a multisymplectic conservation law, which can be viewed as a multidimensional generalization of the symplectic properties inherent in Hamiltonian ODEs. A central question in structure-preserving discretization is determining which finite element methods maintain this property at the discrete level. To address this, we introduce a unified hybridization framework that allows for the simultaneous analysis of local conservation properties across diverse discretization families. We demonstrate that while many methods, such as certain hybridizable discontinuous Galerkin (HDG) schemes, satisfy a strong version of multisymplecticity, the standard conforming finite element exterior calculus (FEEC) methods of Arnold, Falk, and Winther (2006) satisfy this property only in a weak sense. Finally, we establish the bridge between this local multisymplectic perspective and the classical treatment of PDEs as infinite-dimensional Hamiltonian ODEs.

ABSTRACT. The typical energy estimate for the Navier–Stokes equations provides a bound for the gradient of the velocity; energy-stable numerical methods that preserve this estimate preserve this bound. However, the bound scales with the Reynolds number (Re) causing numerical solutions to become unstable (i.e. exhibit spurious oscillations) for large Re.

We propose a mixed finite-element discretisation for the incompressible Navier–Stokes equations, making use of a discrete Stokes complex, that exactly preserves the evolution of both energy and enstrophy (the H1 norm of the velocity). In two dimensions, this includes the strict dissipation of enstrophy, implying a Re-robust bound on the velocity gradient that naturally stabilises the scheme for under-resolved flows. In three dimensions, while the preserved enstrophy evolution equation is not a dissipation inequality (the convective term does not in general vanish), we observe numerically that preserving the behaviour of the enstrophy still has a stabilising effect on the numerical solution.

ABSTRACT. The Koopman operator plays a crucial role in analyzing the global behavior of dynamical systems. Existing data-driven methods for approximating the Koopman operator or discovering the governing equations of the underlying system typically require a fixed set of basis functions, also called dictionary. The optimal choice of basis functions is highly problem-dependent and often requires domain knowledge. We present a novel gradient descent-based optimization framework for learning suitable and interpretable basis functions from data and show how it can be used in combination with EDMD, SINDy, and PDE-FIND. We illustrate the efficacy of the proposed approach with the aid of various benchmark problems such as the Ornstein-Uhlenbeck process, Chua's circuit, a nonlinear heat equation, as well as protein-folding data.

ABSTRACT. The orthogonal group synchronization problem, which aims to recover a set of $d \times d$ orthogonal matrices from their pairwise noisy products, plays a fundamental role in signal processing, computer vision, and network analysis. In recent years, numerous optimization techniques, such as semidefinite relaxation (SDR) and low-rank (Burer-Monteiro) factorization, have been proposed to address this problem and their theoretical guarantees have been extensively studied. While SDR is provably powerful and exact in recovering the least-squares estimator under certain mild conditions, it is not scalable. In contrast, the low-rank factorization is highly efficient but nonconvex, meaning its iterates may get trapped in local minima. To close the gap, we analyze the low-rank approach and focus on understanding when the associated nonconvex optimization landscape is benign, i.e., free of spurious local minima. Recent works suggest that the benignness depends on the condition number of the Hessian at the global minimizer, but it remains unclear whether sharp guarantees can be achieved. In this work, we consider the low-rank approach which corresponds to an optimization problem over the Stiefel manifold ${\rm St}(p,d)^{\otimes n}$. By formulating the landscape analysis into another convex optimization problem, we provide a unified characterization of the optimization landscape for all parameter pairs $(p,d)$ with $p \geq d+2$ for $d\geq 1$ and $p = d+1$ for $1\leq d\leq 3$ which gives a much improved dependence on the condition number of the Hessian. Our results recover the known sharp state-of-the-art bound for $d=1$ which is extremely useful for characterizing the Kuramoto synchronization, and significantly improved the guarantees for the general case $d \geq 2$ with $p \geq d+2$ over the existing results. The theoretical results are versatile and applicable to a wide range of examples.

ABSTRACT. Two-timescale learning algorithms are often applied in game theory and bi-level optimization, using distinct update rates for two interdependent processes. In this talk, I will focus on the two-timescale gradient descent-ascent (GDA) algorithm, which is designed to find Nash equilibria in min-max games with improved convergence properties. Through a PDE-inspired approach, we analyze the convergence of this algorithm for both finite- and infinite-dimensional cases. In finite-dimensional quadratic min-max games, we revisit long-time convergence in near quasi-static regimes through a hypocoercivity perspective. For mean-field GDA dynamics, we investigate convergence under a finite-scale ratio in the whole space by considering a weighted Poincare inequality setup.

ABSTRACT. Solving time-dependent partial differential equations (PDEs) that exhibit sharp gradients or local singularities is computationally demanding, as traditional physics-informed neural networks (PINNs) often suffer from inefficient point allocation that wastes resources on regions already well-resolved. This paper presents an adaptive sampling framework for PINNs aimed at efficiently solving time-dependent partial differential equations with pronounced local singularities. The method employs a residual-driven strategy, where the spatial–temporal distribution of training points is iteratively updated according to the error field from the previous iteration. This targeted allocation enables the network to concentrate computational effort on regions with significant residuals, achieving higher accuracy with fewer sampling points compared to uniform sampling. Numerical experiments on representative PDE benchmarks demonstrate that the proposed approach improves solution quality.

ABSTRACT. We provide new convergence guarantees in Wasserstein distance for diffusion-based generative models, covering both stochastic (DDPM-like) and deterministic (DDIM-like) sampling methods. We introduce a simple framework to analyze discretization, initialization, and score estimation errors. Notably, we derive the first Wasserstein convergence bound for the Heun sampler and improve existing results for the Euler sampler of the probability flow ODE. Our analysis emphasizes the importance of spatial regularity of the learned score function and argues for controlling the score error with respect to the true reverse process, in line with denoising score matching. We also incorporate recent results on smoothed Wasserstein distances to sharpen initialization error bounds.

ABSTRACT. Permutation-structured functions arise in a variety of applications, including symmetric models and antisymmetric wavefunctions in quantum many-body physics. In this talk, I will present a tensor-train (TT) sketching framework for efficiently estimating high-dimensional functions with permutation structure. The proposed approach combines randomized sketching with low-rank tensor representations to mitigate the curse of dimensionality. The framework naturally accommodates both symmetric and antisymmetric settings, with particular emphasis on the latter due to its relevance in Fermionic systems. Numerical experiments illustrate the effectiveness of the approach in capturing high-dimensional structure from limited samples.

ABSTRACT. We present an estimate of the Wasserstein distance between the data distribution and the generation of score-based generative models. The sampling complexity with respect to dimension is $\mathcal{O}(\sqrt{d})$, with a logarithmic constant. In the analysis, we assume a Gaussian-type tail behavior of the data distribution and an $\epsilon$-accurate approximation of the score. Such a Gaussian tail assumption is general, as it accommodates practical target distributions derived from early stopping techniques with bounded support.

The crux of the analysis lies in the global Lipschitz bound of the score, which is shown from the Gaussian tail assumption by a dimension-independent estimate of the heat kernel. Consequently, our complexity bound scales linearly (up to a logarithmic constant) with the square root of the trace of the covariance operator, which relates to the invariant distribution of the forward process.

ABSTRACT. Important nonlinear dynamics, such as those found in plasma and fluid systems, are typically hard to simulate on classical computers. Thus, if fault-tolerant quantum computers could efficiently solve such nonlinear problems, it would be a transformative change for many industries. In a recent breakthrough [Liu et al., PNAS 2021], the first efficient quantum algorithm for solving nonlinear differential equations was constructed, based on a single condition R<1, where R characterizes the ratio of nonlinearity to dissipation. This result, however, is limited to the class of purely dissipative systems with negative log-norm, which excludes application to many important problems. In this work, we correct technical issues with this and other prior analysis, and substantially extend the scope of nonlinear dynamical systems that can be efficiently simulated on a quantum computer in a number of ways. Firstly, we extend the existing results from purely dissipative systems to a much broader class of stable systems, and show that every quadratic Lyapunov function for the linearized system corresponds to an independent R-number criterion for the convergence of the Carlemen scheme. Secondly, we extend our stable system results to physically relevant settings where conserved polynomial quantities exist. Finally, we provide extensive results for the class of non-resonant systems. With this, we are able to show that efficient quantum algorithms exist for a much wider class of nonlinear systems than previously known, and prove the BQP-completeness of nonlinear oscillator problems of exponential size. In our analysis, we also obtain several results related to the Poincaré-Dulac theorem and diagonalization of the Carleman matrix, which could be of independent interest.

ABSTRACT. Computational fluid dynamics (CFD) is a cornerstone of classical scientific computing, and there is growing interest in whether quantum computers can accelerate CFD. In this work, we first show that the existing proposals for fault-tolerant quantum algorithms for CFD suffer from a range of severe bottlenecks that negate conjectured quantum advantages. We then develop a novel algorithm for the incompressible lattice Boltzmann equation that circumvents these obstacles, and provide a detailed analysis of our algorithm, including all potential sources of algorithmic complexity, as well as gate count estimates. We find that for an end-to-end problem, a modest quantum advantage may be preserved for selected observables in the high-error tolerance regime. We lower bound the Reynolds number scaling of our quantum algorithm in dimension D at Kolmogorov microscale resolution with O(Re^{3/4(1+ 𝐷/2)} x q_M), where q_M is a multiplicative overhead for data extraction with q_M = O(Re^{3/8}) for the drag force. This upper bounds the scaling improvement over classical algorithms by O(Re^{3D/8}). However, our numerical investigations suggest a lower speedup, with a scaling estimate of O(Re^{1.936} x q_M) for D=2. Finally, to support our theoretical analysis, we provide a classical numerical study illustrating the accuracy, complexity, and convergence of the algorithm for representative incompressible-flow cases, including the driven Taylor–Green vortex, the lid-driven cavity flow, and the flow past a cylinder. Our results give robust evidence that small, but nontrivial, quantum advantages can be achieved in the context of CFD, and motivate the need for additional rigorous end-to-end quantum algorithm development.

ABSTRACT. Molecular dynamics simulations are a central computational methodology in materials design for relating atomic composition to mechanical properties. However, simulating materials with atomic-level resolution on a macroscopic scale is infeasible on current classical hardware, even when using the simplest elastic network models (ENMs) that represent molecular vibrations as a network of coupled oscillators. To address this issue, we introduce Quantum Elastic Network Models (QENMs) and utilize the quantum algorithm of Babbush et al. (PRX, 2023), which offers an exponential advantage when simulating systems of coupled oscillators under some specific conditions and assumptions. Here, we extend their algorithm in 2D systems and demonstrate how our method enables the efficient simulation of planar materials. As an example, we apply our algorithm to the task of simulating a 2D graphene sheet. We analyze the exact complexity for initial-state preparation, Hamiltonian simulation, and measurement of this material, and provide two real-world applications: heat transfer and the out-of-plane rippling effect. We estimate that an atomistic simulation of a graphene sheet on the centimeter scale, classically requiring hundreds of petabytes of memory and prohibitive runtimes, could be encoded and simulated with as few as $\sim 160$ logical qubits.

ABSTRACT. Quantum computing is often framed in terms of asymptotic speedups under highly idealised assumptions, whereas in practical scientific computing the more immediate question is how quantum methods can be integrated into real simulation workflows to address meaningful problems. This talk will present two frameworks in this direction that seek both to exploit the potential benefits of quantum computing and to confront the subtleties of practical implementation. The first method [1] considers the use of quantum annealing for load balancing in high-performance computing, with applications to adaptive mesh refinement and smoothed particle hydrodynamics, a problem of central importance for scientific codes running on today’s massively parallel HPC systems. The second [2] presents a predictor–corrector reformulation of the HHL algorithm for time dependent simulations, designed to mitigate the classical–quantum data transfer bottlenecks that often limit practical advantage. Although these approaches target different computational tasks, both point toward the same broader conclusion that quantum algorithms may be most effective when deployed selectively to accelerate particularly costly subproblems within larger classical workflows. Achieving this in practice requires careful integration of quantum and classical resources, so that any gains are not offset by surrounding overheads. The talk will discuss the performance of these methods, the constraints imposed by current hardware, and the prospects for incorporating quantum resources into future heterogeneous computing environments.

ABSTRACT. The numerical simulation of systems in which solids immersed in a fluid can come into contact is a challenging problem that raises significant modeling, mathematical, and numerical issues. Even in the absence of contact, fluid-structure interaction (FSI) problems are challenging due to moving geometries and potentially strong coupling between fluid and solid subsystems. When contact is included, additional difficulties arise: (i) in some configurations, FSI models with no-slip boundary conditions fail to predict contact; (ii) the straightforward introduction of contact constraints may lead to mechanically inconsistent models; and (iii) contact induces nonlinearly evolving interface conditions and disparate space and time scales. To address these issues, modified interface conditions and poroelastic roughness modeling have been proposed in the literature. In this talk, we will discuss several of these mechanically consistent modeling approaches and their unfitted mesh approximation.

ABSTRACT. In recent years we have proposed a fictitious domain approach with distributed Lagrange multiplier (FD-DLM) suited for the finite element approximation of fluid-structure interaction problems. The main advantage of this technique is that it allows for the treatment of moving interfaces avoiding the need of body-fitted mesh. We present here the a priori and a posteriori error estimates for the approximation with the FD-DLM method for a simplified stationary coupled problem consisting of Stokes-Stokes interface and Stokes-elliptic problems. In particular, we present the stability properties of our approach leading to optimal error estimates and derive a reliable and efficient error estimator. This represents a first step towards the error analysis of the FD-DLM method applied to FSI problem.

ABSTRACT. In recent years, the framework of Unfitted Finite Element Methods was developed in order to facilitate the handling of complex geometries in Finite Element simulations. We build on these results and present the application of space-time methods to solve moving domain problems. In particular, an isoparametric mapping is generalised in space and time to yield computationally feasible discrete geometries. We start by presenting numerical results of convergence of an arbitrary high order of the proposed method, first at the example of a moving domain convection-diffusion bulk problem. Next, we review the rigorous mathematical proof of this property, involving both a geometry and a discretisation error analysis. Finally, we present the application of these computational tools to a coupled surface-bulk convection-diffusion problem and the transport equation.

ABSTRACT. We construct finite-dimensional de Rham subcomplexes generated by fixed-neuron shallow ReLUk neural networks, a class of spaces known to provide optimal approximation rates. For neurons of the form s_i(x)=wi · x+bi, the space of neural p-forms is built from the ridge functions ReLU^k-p(s_i) together with constant-coefficient p-forms. These spaces are compatible with the exterior derivative because differentiating a ReLU power lowers its order by one, and for each fixed neuron, differentiation amounts to exterior multiplication by the fixed one-form ds_i. Under a linear independence assumption on the lowest-order family {ReLU^k-p(s_i)}, the global complex decomposes into independent neuron-wise Koszul complexes. We prove exactness in arbitrary dimension and provide a geometric sufficient condition for the required linear independence.

ABSTRACT. Full-waveform inversion (FWI) and least-squares reverse time migration (LSRTM) are powerful tools for high-resolution subsurface reconstruction and imaging. However, FWI is highly sensitive to the initial model because of local minima, while LSRTM often suffers from slow convergence and high computational cost due to repeated wavefield simulations in iterative updates. In this work, we investigate how source wavelets influence the optimization landscape of FWI and the conditioning of the linear inverse problem underlying LSRTM, and develop a unified source-manipulation framework to improve both tasks. The proposed approach decomposes the original problem into two stages. In the first stage, the recorded data are transformed into equivalent data associated with a desired source wavelet. For the known-source case, this source transformation is formulated as a regularized deconvolution problem; for the unknown-source case, it is handled by a convolutional neural network that directly maps the observed data to the target-source data. In the second stage, the transformed data are used in a conventional FWI or LSRTM procedure. For LSRTM, our analysis further shows that properly designed source wavelets lead to a more favorable eigenvalue distribution of the normal operator and thus improved conditioning and faster iterative convergence. Numerical experiments on benchmark models demonstrate that the proposed source-manipulation strategy improves gradient quality for FWI, accelerates LSRTM convergence, reduces the computational burden, and yields higher-resolution reconstructions and images in both known- and unknown-source scenarios.

ABSTRACT. We develop the $L^p$ space and approximation theory for shallow neural networks with $\mathrm{ReLU}^k$ activations. The central object is the Radon-domain $L^p$ space $\mathcal{R}L^p_k(\Omega)$ containing all functions on a bounded domain $\Omega$ that admit a ridge integral representation whose coefficient density belongs to $L^p$ in the Radon domain. In the Hilbert case $p=2$, we prove by elementary Fourier analysis that this space recovers the critical Sobolev space $H^{k+(d+1)/2}(\Omega)$. For general $1<p<\infty$, the identity becomes a Sobolev sandwich. The sharp gap of each side is exactly the Seeger--Sogge--Stein loss for the Radon transform as a Fourier integral operator. This also clarifies how the activation regularity and Radon back-projection jointly produce the regularity. As an application, we discretize the integral representation using a deterministic interpolation skeleton plus uniform sampling. This yields high-probability $L^p$ approximation rates and the optimal Hilbert rate $O\!\big(n^{-\frac12-\frac{2k+1 {2d}}\big)$ at $p=2$ for linearized neural networks.

ABSTRACT. The identification of Partial Differential Equations (PDEs) has emerged as a prominent data-driven approach for mathematical modeling and has attracted considerable attention in recent years. The stability and precision in identifying PDE from heavily noisy spatiotemporal data present significant difficulties. This problem becomes even more complex when the coefficients of the PDEs are subject to spatial variation. In this paper, we propose a Weak formulation of Groupsparsity-based framework for IDENTifying PDEs with varying coefficients, called WG-IDENT, to tackle this challenge. Our approach utilizes the weak formulation of PDEs to reduce the impact of noise. We represent test functions and unknown PDE coefficients using B-splines, where the knot vectors of test functions are optimally selected based on spectral analysis of the noisy data. To facilitate feature selection, we propose to integrate group sparse regression with a newly designed group feature trimming technique, called GF-Trim, to eliminate unimportant features. Extensive and comparative ablation studies are conducted to validate our proposed method. The proposed method not only demonstrates greater robustness to high noise levels compared to state-of-theart algorithms but also achieves superior performance while exhibiting reduced sensitivity to hyperparameter selection.

ABSTRACT. Over the past two decades, the mathematical modelling of Fractional Differential Equations (FDEs) has attracted growing attention and expanded significantly across various scientific fields. In this context, finding analytical solutions is often more challenging than for classical ordinary differential equations, while accurate and reliable numerical methods can be hindered by the potential nonsmoothness of the solution and/or vector field at the initial time. Moreover, the nonlocality of the differential operator and the persistence of the intrinsic memory term can render long-time simulations computationally demanding. To mitigate these issues, the class of Runge-Kutta type methods, known as Fractional Hamiltonian Boundary Value Methods (FHBVMs), is presented, covering its design, development and analysis. In particular, a novel extension is discussed, allowing for a mixed graded/uniform mesh for time step selection and resulting in an updated version of pre-existing Matlab codes. Numerical experiments show that this approach is especially effective for problems having nonsmooth vector field/solution at the initial time, with solution of oscillatory type, achieving higher accuracy in reproducing the initial nonsmooth behavior, while maintaining efficiency over long time periods. Finally, a generalisation to fractional multi-order problems is introduced and applied to model predator-prey dynamics with intraguild predation, effectively accounting for potentially different rates of change of the populations with respect to their own time history.

ABSTRACT. We consider a stochastic time-fractional Burgers equation driven by additive fractional Gaussian noise with Hurst parameter $H>\frac{1}{2}$. The interplay between the Caputo fractional derivative, the non-globally Lipschitz convection term, and the low regularity of the noise presents substantial analytical challenges. We propose a fully discrete numerical scheme that combines a spectral Galerkin method in space with a backward Euler convolution quadrature in time. Due to the lack of semimartingale structure in fractional Brownian motion and the temporal singularity induced by the fractional derivative, classical stochastic calculus techniques are inapplicable. Instead, we establish pathwise estimates of $L^2$-norm using a generalized Henry-type Gronwall inequality and derive probabilistic error bounds via the Markov inequality. Under suitable regularity assumptions, we prove the convergence in probability of the numerical solution. The analysis highlights the delicate balance between stochastic regularity and fractional temporal dynamics. Numerical experiments are provided to illustrate the convergence behavior and verify the theoretical results.

ABSTRACT. We study a class of fractional delay differential equations involving the Caputo–Fabrizio derivative with non-singular exponential kernel. The problem is reformulated as an equivalent integral equation, which allows the application of progressive contraction techniques to obtain existence and uniqueness results under Lipschitz-type assumptions. Moreover, we establish continuous dependence on the initial condition, prove Ulam–Hyers stability, and analyze the asymptotic behavior of solutions as the delay parameter tends to infinity.

ABSTRACT. We propose a second-order numerical framework for variable-order (VO) time-fractional reaction--diffusion equations of the form \[ _SD_t^{\alpha(t)}u = \mathcal{L}u + F(u), \] where the fractional order evolves according to \[ \alpha(t)=\alpha_2+(\alpha_1-\alpha_2)e^{-ct}, \qquad 0<\alpha_1,\alpha_2<1, \] and \(\mathcal{L}\) denotes either the classical diffusion operator \(\Delta\) or the nonlocal fractional diffusion operator \(-(-\Delta)^s\). The nonlinear reaction term \(F(u)\) is chosen from representative models including Fisher--KPP and Gray--Scott equations.

The main difficulty is that the VO operator $_SD_t^{\alpha(t)}$ has a time-dependent memory structure, which makes its numerical treatment substantially more challenging than in the constant-order case, especially when both accuracy and efficiency for long-time simulations are required. This challenge is further compounded by the need to handle nonlinear reactions and both local and nonlocal diffusion within a unified framework. To this end, we develop a second-order IMEX time-stepping scheme based on BDF2 convolution quadrature. To make the method practical for long-time simulations, an adaptive compressed-history technique is introduced to reduce the cost of evaluating the memory term. In space, periodic boundary conditions are imposed, and second-order spatial discretizations are used for diffusion operators.

At each time step, the resulting linear systems are solved by GMRES with a specially designed preconditioner, which significantly improves the computational efficiency of the fully discrete method. Numerical experiments demonstrate second-order convergence and illustrate the significant effects of VO memory and nonlocal diffusion on front propagation and pattern formation.

ABSTRACT. We construct a finite element method for the numerical solution of a fractional porous medium equation on a bounded open Lipschitz polytopal domain, where d = 2 or 3. The pressure in the model is defined as the solution of a fractional Poisson equation, involving the fractional Neumann Laplacian in terms of its spectral definition. We perform a rigorous passage to the limit as the spatial and temporal discretization parameters tend to zero and show that a subsequence of the sequence of finite element approximations defined by the proposed numerical method converges to a bounded and nonnegative weak solution of the initial-boundary-value problem under consideration. This result can be therefore viewed as a constructive proof of the existence of a nonnegative, energy-dissipative, weak solution to the initial-boundary-value problem for the fractional porous medium equation under consideration, based on the Neumann Laplacian. The convergence proof relies on results concerning the finite element approximation of the spectral fractional Laplacian and compactness techniques for nonlinear partial differential equations, together with properties of the equation, which are shown to be inherited by the numerical method. We also prove that the total energy associated with the problem under consideration exhibits exponential decay in time.

ABSTRACT. We develop a hypercomplex operator calculus to treat fully analytically boundary value problems for the homogeneous and inhomogeneous fractional Helmholtz equation, where fractional derivatives in the sense of Caputo and Riemann–Liouville are applied. Our method extends the recently proposed fractional reduced differential transform method (FRDTM) by using fractional derivatives in all directions. For the special separable case in three dimensions, we obtain completely explicit representations for the fundamental solution.

ABSTRACT. Equilibrium properties in statistical physics are obtained by computing averages with respect to Boltzmann-Gibbs measures, sampled in practice using ergodic dynamics such as the Langevin dynamics. Some quantities however cannot be computed by simply sampling the Boltzmann-Gibbs measure, in particular transport coefficients, which relate the current of some physical quantity of interest with the forcing needed to induce it. For instance, a temperature difference induces an energy current, the proportionality factor between these two quantities being the thermal conductivity. From an abstract point of view, transport coefficients can also be considered as some form of sensitivity analysis with respect to an added forcing to the baseline dynamics.

There are various numerical techniques to estimate transport coefficients, which all suffer from large errors, in particular large statistical errors. I will review the most popular methods, namely the Green--Kubo approach where the transport coefficient is rewrittten as some time-integrated correlation function, and the approach based on longtime averages of the stochastic dynamics perturbed by an external driving (so-called nonequilibrium molecular dynamics). I will make precise in each case the various sources of errors, in particular the bias related to the time discretization of the underlying continuous dynamics, and the variance of the associated Monte Carlo estimators. I will also briefly present some recent alternative techniques to estimate transport coefficients.

ABSTRACT. In this talk, we introduce a new approach that integrates Kolmogorov–Arnold Networks (KANs) with linear multistep methods (LMMs) for discovering and approximating dynamical systems. We first provide approximation and error bounds for two-layer B-spline KANs when representing the vector fields of dynamical systems. Based on these results, we show that for certain classes of LMMs, the total error can be described in terms of both the step size of the method and the approximation accuracy of the network. We also study the difference between the numerical solutions generated from the learned vector fields and the true trajectories of the systems.

ABSTRACT. Coarse-graining and model reduction are essential for extending the accessible time scales of molecular and stochastic simulations, yet classical reduced models often fail to simultaneously preserve equilibrium statistics and reproduce dynamical observables. This limitation is particularly important for scientific machine learning.

In this talk, I will present a recently developed Markovian coarse-graining framework for overdamped Langevin dynamics. Instead of relying solely on the classical potential-of-mean-force approximation, we construct improved reduced models through a hierarchical learning procedure that combines force corrections with locally inferred spatiotemporal rescaling rules. This approach incorporates energy-landscape information into data-driven effective dynamics. We prove that the resulting model preserves equilibrium statistics while eliminating systematic errors in dynamical quantities such as the mean-squared displacement. These results establish a rigorous scientific machine learning framework for interpretable, physics-consistent, and dynamically accurate coarse-grained modeling of stochastic multiscale systems.

This is a joint work with Dr Thomas Hudson from the University of Warwick.

ABSTRACT. Many Hamiltonian systems with energy-dissipative terms possess a conformal symplectic structure, wherein the symplectic form is proportionally preserved. In this study, we propose Conformal Symplectic Neural Networks that learn the solution maps for different energy-dissipative systems by treating the Hamiltonian as input information.

ABSTRACT. Sampling methods based on Langevin dynamics have applications in nonequilibrium molecular dynamics, Bayesian inference and machine learning. We focus on the noise and the deviation from the invariant distribution caused by the stochastic gradient methods, and then propose the idea of using variable stepsize to control gradient noise. We design an adaptive stepsize adaptive Langevin (Ad2L) algorithm. In combination with the idea of thermostats, the dynamics are changed by introducing the monitor function and an invariant measure-preserving transformation. We apply the splitting method to implement a discretised numerical method to achieve the purpose of dynamically adjusting the stepsize to cope with the problem of steep gradient changes and variable noise. Numerical experiments show that the Ad2L algorithm not only performs well in toy models with steep gradient and large noise, but also effectively improves the convergence speed and achieves higher sampling accuracy when applied to higher-dimensional problems.

ABSTRACT. We consider a generic class of stochastic particle-based models whose state at an instant in time is described by a set of continuous degrees of freedom (e.g., positions), and the length of this set changes stochastically in time due to birth-death processes. Using a master equation formalism, we write the dynamics of the corresponding (infinite) set of probability distributions: this takes the form of coupled Fokker-Planck equations with model dependent source and sink terms. We derive the general expression of the entropy production rate for this class of models in terms of path irreversibility. To demonstrate the practical use of this framework, we analyze a biologically motivated model incorporating division, death, and diffusion, where spatial correlations arise through the division process. By systematically integrating out excess degrees of freedom, we obtain the marginal probability distribution, enabling exact calculations of key statistical properties such as average density and correlation functions. We validate our analytical results through numerical Brownian dynamics simulations, finding excellent agreement between theory and simulation. Our method thus provides a powerful tool for tackling previously unsolved problems in stochastic birth-death dynamics.

ABSTRACT. The clustering of self-motile and repulsive particles, so-called motility-induced phase separation (MIPS), is one of the clearest signatures of active physics. Typically, increasing the amplitude of self-motility increases the degree of clustering, however for high enough self-motility the homogeneous phase is reentered. Here, we numerically report that such reentrance naturally emerges in an overdamped Hamiltonian model, with such a model recapitulating properties of (active) bird flocks, and exhibits clustering behaviour reminiscent of MIPS. We demonstrate the reentrance of the homogeneous phase and identify the underlying mechanism as a competition between the amplitude of a spin-velocity coupled drive and mobility-limited kinetic frustration. We reveal, both numerically and analytically, that strong spin-velocity coupling suppresses transverse diffusion, which causes a dimensional reduction for the particle dynamics, thereby leading the system into an arrest that closes the window for phase separation. Overall, our work offers a Hamiltonian bridge between reentrant physics across equilibrium and non-equilibrium materials.

ABSTRACT. Proliferation is a distinctive source of activity in soft and biological matter: It creates new degrees of freedom, generates stresses and flows, changes intrinsic length scales, and couples local mechanical events to collective dynamics. These features pose challenges both at the level of numerical simulation, where growth and division alter the structure of the system itself, and at the level of statistical-physics analysis, where changing particle numbers, finite cell lifetimes, and expanding or remodeling geometries complicate the definition of trajectories, transport coefficients, and effective parameters. In this talk, I will discuss how particle-based simulations and effective models inspired by active Brownian particles can be combined to identify transition mechanisms in proliferating cellular active matter. Numerical consistency is essential in this setting: Aberrant force discontinuities, discontinuous trajectories, or loss of cell identity can directly affect mechanical observables, stresses, and displacement statistics. I will discuss modeling strategies for dividing cells based on smooth, force-based dynamics, which provide a controlled way of following mechanical observables and trajectories across division events. Beyond division itself, additional cellular activities require similarly careful modeling, such as cell motility based on crawling through or across neighboring cells, which might not be adequately represented by traditional self-propulsion and instead requires reciprocal forces. The broader goal then is to use large-scale simulations to expose the macroscopic consequences of proliferation and extract the physical mechanisms underlying them. Following this strategy, our studies have revealed that, in expanding cell aggregates, local rearrangements can coexist with lineage confinement, so that motility must overcome steric suppression and expansion by growth before global mixing is possible. In homeostatic mixtures of proliferating and motile particles, on the other hand, the proliferating component can act as a dynamic medium that renormalizes motility and persistence, enhances diffusion, and mediates effective interactions, producing a new type of phase separation. On the technical side, these examples demonstrate that effective models, such as active-Brownian-particle-like descriptions with parameters shaped by proliferation, provide a bridge from detailed simulations to coarse-grained physical mechanisms. Conceptually, our results show how transitions in proliferating active matter may emerge from the competition between active motion, mechanically generated constraints, and population-level growth or turnover.

ABSTRACT. Four-dimensional data assimilation (4D-Var) aims to find the best state of a dynamical system by minimizing an objective function that measures the distance of a model trajectory to observations over a given time window, while remaining close to a prior forecast. In practice, the problem is usually solved using an approximate Gauss-Newton iterative method to minimize the function. However, the time window over which the minimization is performed is limited by the error growth rate of the model. A difficulty arises when applying 4D-Var to a coupled dynamical system with different timescales, as the time window must be chosen based on the fastest error growth. We consider this problem in the context of coupled atmosphere-ocean prediction, now used at operational weather forecasting centres. For this problem the chosen time window must reflect the fast error growth in the atmosphere, which limits the information that can be extracted from observations of the slower ocean. Here, a new method is proposed to address this issue. After running the standard coupled assimilation for several short windows, a new 'smoother' procedure is performed to find an ocean trajectory that best fits all ocean observations over these windows. This allows observations at later assimilation windows to update the state at earlier windows. Furthermore, it permits the use of late-arriving observations that were not available to the short-window assimilations. Using an idealized model with two timescales we illustrate how the new method is able to improve the estimate of the state and the subsequent coupled model forecast.

ABSTRACT. For partially observed systems depending on unknown parameters, applying the Ensemble Kalman Filter (EnKF) equations to an augmented state–parameter formulation induces a dynamic on the parameter that empirically steers the ensemble towards the true parameter, thereby providing parameter identification in addition to state estimation. When the state arises from a PDE discretisation, the joint estimation problem becomes computationally challenging due to expensive, typically nonlinear, forward solves required for each ensemble member. However, filtered systems frequently display a low-rank structure: correcting along a few dominant, time-varying directions often suffices to ensure adequate signal tracking.

We propose a Dynamical Low-Rank (DLR) filtering framework for state/parameter estimation that evolves the filtering density in a low-dimensional, dynamically adapting subspace. Building on existing DLR Kalman filter theory for linear systems, we extend the approach to the nonlinear setting. For efficient time-stepping, we combine a forecast/analysis scheme with a Basis Update and Galerkin (BUG) integrator, enabling efficient in-subspace assimilation with rank-adaptive variants. We furthermore introduce a hyperreduction technique to cheaply evaluate the nonlinear forward maps.

We demonstrate the method on a Fisher–KPP system with unknown diffusion coefficient and on a reduced human arterial network with 55 arteries observed at 3 locations. Both test cases confirm computational gains while preserving filtering accuracy comparable to the full order model EnKF.

ABSTRACT. Solving Bayesian inverse problems governed by partial differential equations (PDEs) is notoriously expensive, as each likelihood evaluation demands a PDE solve. This talk presents Gaussian process (GP) surrogates tailored for linear PDEs that deliver accurate posterior approximations in the small-data regime. Building on Raissi et al. (2017), we construct PDE-informed GP priors that couple the unknown parameters and spatial eld through the operators in the PDE, yielding a joint prior over the solution, source, and boundary data. This enables training with cheap, auxiliary physics data and admits closed-form gradients for e cient MCMC via the Metropolis adjusted Langevin algorithm.

ABSTRACT. We consider the common setting in which samples from a target distribution are obtained via a poorly converged sampling method, such as an MCMC chain with limited mixing or constrained by expensive likelihood evaluations. This frequently arises in inverse problems where each density evaluation requires the solution of a differential equation. Although the resulting sample set may be insufficient for reliable Monte Carlo estimation, exact evaluations of the target density are typically available at the sampled points.

We introduce Posterior Regression ImpOrtance Reweighting (PRIOR), a method that exploits this additional information by learning a surrogate of the posterior, for example using Gaussian processes or neural networks. The parameter space is partitioned into cells, and sample weights are adjusted so that the total contribution of samples within each cell matches the posterior mass of that region under the surrogate. This provides a mechanism to correct for poor sample coverage and imbalance.

Our analysis yields error decompositions that inform the adaptive construction of the partition, allowing the method to concentrate resolution where it most improves accuracy. We present initial numerical experiments on simple test problems that illustrate the behaviour of the method.

ABSTRACT. Optimal control problems arise in a wide range of applications, including aerodynamics, mathematical finance, bioprocesses, and epidemiology. Using the Lagrange multiplier technique, optimal solutions can be characterized by the first-order optimality system. When the governing PDEs are time-dependent, this system typically exhibits a forward-backward structure, and classical time-stepping methods cannot be applied to solve this system. Solving the entire system at once can become computationally expensive, especially in higher spatial dimensions. To address this challenge, parallelization techniques are crucial. In this talk, I will present some recent developments using domain decomposition methods, supported by both theoretical results and numerical examples.

ABSTRACT. In this talk, we consider the theoretical convergence of flexible GMRES. While convergence of standard GMRES is well studied, there exist few results of similar nature for flexible GMRES. The aim of this talk is to discuss and fill in this gap. In addition, we report on experiences of using these ideas in the context of randomized sketching.

ABSTRACT. Standard boundary integral equation techniques do not naturally result in 2nd-kind integral equations for magnetostatics in the presence of an open surface. We begin by reviewing the classical formulation for closed surfaces. Then we focus on numerical solvers for axisymmetric closed surfaces, where the axisymmetry yields substantial computational advantages. Next, we formulate a second-kind integral equation for the genus-zero open-surface case. Finally, we introduce a consistency condition and augmented representation of the fields that ensures well-posedness when the boundary genus exceeds zero. We will take a look at the numerical difficulties and future work that arise from these problems.

ABSTRACT. Moving interface problems arise across biology, fluid dynamics, and materials science, where evolving geometries must be captured accurately to model physical behaviour. This talk introduces a selection of Arbitrary Lagrangian–Eulerian (ALE) schemes tailored to such problems, focusing on clear formulation and practical performance. I will outline several moving mesh strategies—motivated by both geometric and analytical considerations—and show how they integrate interface tracking with bulk discretisation while maintaining mesh quality under large deformations. Computational examples will illustrate the accuracy, robustness, and efficiency of the proposed methods, as well as the trade-offs between different ALE formulations.

ABSTRACT. In this talk I would like to discuss error estimates for anisoptropic mean curvature flow of closed surfaces. We will first derive the evolution equations for the anisotropic flow, which formally resemble to those for the mean curvature flow (isotropic case), yet have substantially more complicated proofs. The discretisation in space uses evolving surface finite elements. Thanks to this formal resemblance, stability and consistency proofs are quite similar to previous works, and they will lead to optimal-order H^1-norm error estimates. Various numerical experiments and some regularisation ideas will also be presented.

The talk is based on joint work with Harald Garcke (Regensburg) and Klaus Deckelnick (Magdeburg).

ABSTRACT. We present a sliding-interface finite element formulation for fluid--structure interaction between an incompressible fluid and a rotating rigid body. The method combines a rotational arbitrary Lagrangian--Eulerian framework with a skew-symmetric Nitsche stabilization imposed on an artificial sliding interface. This design preserves an energy-dissipating property at the continuous level and leads to a first order fully discrete scheme that retains the same dissipation mechanism, providing a stable and accurate time-marching strategy. On the theoretical side, we address a key gap in the literature by establishing an inf--sup condition for the isoparametric FEM on non-matching and overlapping meshes across the sliding interface, extending classical results beyond fitted-interface settings. Based on this inf--sup stability and the energy estimate, we prove unique solvability of the fully discrete scheme. Numerical experiments in both two and three dimensions demonstrate convergence, efficiency, and consistent energy dissipation of the proposed approach.

ABSTRACT. In this talk, we present a numerical analysis of the Eyles-King-Styles tumor growth model, a free boundary problem coupling a Poisson equation in the bulk \Omega with a forced mean curvature flow on its boundary \Gamma. Unlike existing evolving surface analyses based on integer-order Sobolev spaces, this bulk-surface coupling requires H^{1/2}-order regularity on \Gamma. We establish a fractional Sobolev framework that admit a rigorous convergence analysis for continuous finite elements of polynomial degree at least three.

ABSTRACT. In centrifuged convection, the convection cell is placed at the end of a long arm that rotates rapidly about an axis orthogonal to the arm. The normals of the hot and cold walls are parallel to the arm, producing a centrifugal acceleration that can greatly exceed gravitational acceleration. This configuration differs from rotating convection with gravity parallel to the rotation axis, where centrifugal buoyancy is weak and often neglected. In this talk, we will present the governing equations, numerical scheme, and flow dynamics of centrifugal-acceleration-dominated flows inside a rotating cube across various Prandtl and Ekman numbers. We will also discuss how the heat flux varies with changes in these dimensionless parameters.

ABSTRACT. Minimum action methods provide a powerful framework for analyzing rare transitions in small-noise-driven dynamical systems, but their practical performance is often limited by time truncation and parameter sensitivity in infinite-horizon problems. In this paper, we develop an efficient Laguerre spectral minimum action method (LMAM) for computing quasi-potentials associated with fixed points of dynamical systems. Based on the large deviation framework, the method computes minimum action paths by formulating the problem on a semi-infinite time interval and discretize the temporal direction using Laguerre functions. An adaptive time rescaling strategy based on our recent theoretical framework [SINUM 64,125–147; arXiv.2602.03083.642] is proposed to enhance accuracy and convergence of the Laguerre spectral approximation. To efficiently handle nonlinear terms, we employ an improved procedure for evaluating Laguerre--Gauss--Radau quadrature [J. Sci. Comput. 101:72. doi:10.1007/s10915-024-02725-9], which enables stable and accurate double-precision computations with a large number of Laguerre modes. Precise numerical analysis for the linear problem and a local result for the nonlinear case are developed. Numerical experiments including both ordinary and partial differential equations (Allen-Cahn and Navier-Stokes) are presented to illustrate the accuracy and efficiency of the proposed method.

ABSTRACT. Long-time integration of high-dimensional polynomial dynamical systems requires both accurate error control and preservation of geometric structures. This work develops an algebraic framework for constructing structure-preserving splitting methods using generalized ridge representations. We study finite block representations of real polynomial functions and polynomial vector fields under the parameters n, d, and k, and connect these representations with locally solvable subflows and geometric compositions.

For homogeneous polynomial Hamiltonians with k = 1, we construct a fixed family of prime-Vandermonde directions and prove that every real homogeneous polynomial of degree d in n variables admits a unique representation as a sum of d-th power terms along these directions. For Hamiltonian systems, one explicit Euler step for each one-dimensional power block is shown to coincide with the exact flow map, and Lagrangian blocks of dimension n/2 are used to reduce the number of blocks. For source-free systems, a skew-symmetric matrix-valued potential is obtained through a homotopy formula, leading to explicit volume-preserving shear splittings and fixed-sequence partitioned Runge-Kutta compositions. For contact systems, we identify structural obstructions for finite-linear-function blocks and construct analytically solvable contact monomial blocks on the positive domain.

Numerical experiments with high-accuracy DOP853 reference solutions confirm the construction procedure and the preservation of symplectic, volume-preserving, and contact structures. Comparisons with standard fourth-order geometric integrators show that the proposed block compositions usually reduce CPU cost while maintaining comparable trajectory accuracy.

ABSTRACT. The spectral collocation method is easy to implement--one simply replaces derivatives by spectral differentiation matrices at Gaussian-type collocation points. This plug-and-play nature makes it well suited for variable-coefficient and nonlinear problems. However, the collocation method is often severely ill-conditioned, which limits its application, particularly in multiple dimensions for solutions that require high resolutions. In this talk, we introduce optimal preconditioners for multi-dimensional collocation methods by using the notion of preconditioning differentiation by integration via Birkhoff interpolation and stable matrix diagonalization techniques. With this at our disposal, we discuss the preconditioning of the Hierarchical Poincare–Steklov (HPS) spectral-element method, which is a fast and accurate direct spectral solver using the collocation methods as local solvers.

ABSTRACT. Power system dynamics is one of the most computationally challenging areas in energy system modeling and analysis. In the wake of a fault, relay trip, or equipment malfunction, electrical surges propagate and refract along high voltage transmission lines at nearly the speed of light. Power system dynamics is modeled by networks of connected differential algebraic equations, typically of index-1.

We propose efficient multirate infinitesimal integrators for large-scale dynamic power grid simulations. The new multimethod framework leverages the system structure and separates the algebraic solutions, which involve variables from all subsystems, form the time advancement of subsystem DAEs, which are integrated in parallel and with different discretization schemes and step sizes.

ABSTRACT. The real-time Boltzmann transport equation (rt-BTE) provides a first-principles framework for understanding nonequilibrium dynamics in materials to enable predictive modeling of charge and heat transport as well as ultrafast pump-probe experiments. The widely separated timescales of electron-phonon (e-ph) and phonon-phonon (ph-ph) interactions and the high cost of evaluating collision integrals has limited simulations to short times and small system sizes. In this presentation, we discuss the application of adaptive multirate time integration methods from the SUNDIALS library to accelerate rt-BTE simulations in PERTURBO. Multirate infinitesimal (MRI) methods with step size adaptivity at both the fast and slow time scales, allow for the relatively inexpensive, fast e-ph dynamics and the more expensive, slower ph-ph dynamics of the system to be advanced with different, dynamically selected time steps. The resulting increase in efficiency leads to 10x faster results in 2D graphene simulations and enables simulating bulk materials (silicon and gallium arsenide) that were previously intractable. These results demonstrate that adaptive multirate integration can significantly expand the scope of first-principles nonequilibrium simulations in materials.

ABSTRACT. Nonlinearly partitioned Runge–Kutta (NPRK) methods are a new family of time-integration schemes that generalize additive and component-partitioned Runge–Kutta methods. Specifically, NPRK methods allow different factors within nonlinear terms to be treated with differing levels of implicitness. In this talk, I will provide an overview of the NPRK framework, present new methods with varying types of implicitness (implicit–explicit and implicit–implicit), and discuss several recent developments regarding multirate NPRK methods.

ABSTRACT. Due their mechanism of generating high-order approximations from lower order information, Runge-Kutta methods can incur order reduction when applied to stiff problems. Implicit-Explicit Runge-Kutta (ImEx RK) methods fundamentally benefit from a lower triangular structure, and therefore high stage order is generally not a viable option to recover the full convergence order. This talk first characterizes how order reduction manifests in classical ImEx RK methods, and then presents how stiff order conditions can be enforced in new ImEx RK methods that mitigate order reduction while maintaining a lower triangular structure.

ABSTRACT. Time integration of nonlinear Dirac equations in the nonrelativistic limit regime is notoriously challenging, because the solution oscillates in time with high frequency. When traditional time integrators are applied to such problems, an acceptable accuracy can only be expected if the oscillations are resolved by many tiny time steps, which is hopelessly inefficient.

In this talk, we present a new splitting method which, in contrast to standard splitting methods, yields high accuracy even when the step size is larger than the wavelength of the oscillations. We discuss three different bounds for the local error, which correspond to three different parameter regimes. Numerical experiments indicate, however, that surprisingly the accuracy of the new integrator is even better than what can be expected from our analysis.

ABSTRACT. Runge-Kutta methods are widely known and used for the numerical integration of differential equations of the form x'=f(x). When advancing one step, these methods use several intermediate values of x, denoted by k, and in general there are as many k's as stages in the method. This can be a drawback for problems where the dimension of x is very large, as it may lead to excessive memory usage. For this reason, Runge-Kutta methods have been developed that can be implemented using only two registers per iteration. These are known as Low-storage Runge-Kutta methods. In this presentation, we introduce a new approach in which splitting methods can be used as Low-storage methods, and we discuss the advantages they offer compared to classical Low-storage Runge-Kutta methods. This is a joint work with Sergio Blanes and Abutalib Mohammed

ABSTRACT. This study introduces a new idea of splitting methods for solving singular and nonlinear reaction-diffusion partial differential equations involving mixed derivative terms. Intercardinal splitting finite difference approximations are built for fast and accurate simulation approximations of underlying solutions. The exploration will demonstrate that the implicit schemes accomplished are numerically stable, convergent, efficient, and preservative for key physical parameters such as positivity and monotonicity. Dynamic orders of accuracy of the numerical algorithms will be estimated through generalized Milne’s devices. Simulation experiments of the intercardinal splitting methods will also be presented. Further open problems will be outlined.

ABSTRACT. In this talk we provide a comprehensive overview of numerical techniques for integrating the time-dependent Schrödinger equation with an autonomous potential. These methods fall into two main categories: polynomial-based approximations of the exponential operator and methods that exploit the separable structure of the Hamiltonian. In the first category, we discuss Taylor and Chebyshev expansions, the Lanczos procedure, and a recent algorithm based on symplectic splitting schemes; all of these are subject to step size restrictions tied to spatial discretization. The second category comprises unitary splitting methods, which are unconditionally stable and inherently preserve the unitary evolution, although numerical resonances may arise for specific step sizes. We analyze the key properties of each approach, illustrate their practical performance through representative examples and provide some guidance about their use for practical applications.

ABSTRACT. Geometric flows, where an immersed manifold evolves in time according to its own geometry, exhibit important structural properties.For example, surface diffusion dissipates surface area while conserving volume; it is desirable to preserve these properties on discretization. This has motivated a substantial body of research on structure-preserving discretizations for these flows, albeit at low order in time. In this work, we present the first discretization of geometric curvature flows (curve shortening/mean curvature flow and curve/surface diffusion) that preserves the evolution of area and volume at arbitrary order in space and time. The key idea is to introduce auxiliary variables in a particular way so that the derivation of the area dissipation law can be replicated after discretization with continuous Petrov--Galerkin in time. These auxiliary variables are indicated by a general strategy for structure-preservation in time that applies to many other problems. The proposed scheme also preserves mesh quality in the same manner as the minimal deformation rate strategy. We demonstrate its structure-preserving properties and high-order convergence on several benchmark examples.

ABSTRACT. Structure-preserving discretisations, by construction, enable the conservation of many invariants of the dynamics of PDEs. For example, energy, momentum, helicity (or en- strophy in 2D), in the Navier-Stokes equations, also potential enstrophy in the shallow waters equations, magnetic and cross helicity in MHD, etc. Besides the formal relevance, there are also concrete practical benefits of employing numerical discretisations that preserve some of the structure of the original equations. Structure preserving discretisations fundamentally rely on a set of function spaces that constitute a discrete de Rham sequence and is extensively explored in Finite Element Exterior Calculus (FEEC). In this work we will consider the linear wave equation as a prototypical case to explore the extension to space-time FEEC discretisations.

ABSTRACT. The construction of numerical optimal control methods relies on a discretization of the underlying dynamical system. Much effort in designing numerical methods for the simulation and optimization of dynamical systems has been put into schemes which are structure preserving. One particular class are variational integrators which are momentum preserving and symplectic. Whereas for Lagrangian or Hamiltonian dynamics the canonical symplectic form is defined on the space of configuration and momenta, in optimal control theory, there exists a sympelctic structure defined on the space of state and adjoint variables stemming from the necessary optimality conditions, which is preserved for the optimal solution. In recent years there have been many works on analyzing and preserving symplecticity in the state adjoint system of the optimal control problem. However, typical approaches do not allow for a consistent and unified formulation of the optimal control problem and its discretization.

In this talk, we present a new framework that restates the optimal control problem using a Lagrangian formulation, whose Euler-Lagrange equations provide the necessary conditions for optimality. The new Lagrangian is nondegenerate which allows to naturally derive the Hamiltonian version of necessary optimality conditions based on the Legendre transform. Furthermore, the development of variational integrators for the new Lagrangian naturally leads to symplectic discretizations of the necessary conditions defined on the state-adjoint space. If the underlying dynamical system is itself described by forced Euler-Lagrange equations, the resulting numerical schemes are not only symplectic in the state-adjoint space, but automatically yield a symplectic discretization of the dynamic state equations, leading to doubly symplectic schemes. We will discuss geometric and numerical properties of our new formulation and illustrate the theoretical results by means of examples for controlled mechanical systems.

ABSTRACT. Adjoint systems are widely used to inform control, optimization, and design in systems described by ordinary differential equations or differential-algebraic equations. We explore the geometric properties and develop geometric discretization methods for such adjoint systems. In particular, we utilize symplectic and presymplectic geometry to investigate the properties of adjoint systems associated with ordinary differential equations and differential-algebraic equations, respectively. We show that the adjoint variational quadratic conservation laws, which are key to adjoint sensitivity analysis, arise from (pre)symplecticity of such adjoint systems. We discuss various additional geometric properties of adjoint systems, such as symmetries and variational characterizations.

As an application of this geometric framework, we discuss how the adjoint variational quadratic conservation laws can be used to compute sensitivities of terminal or running cost functions. Furthermore, we develop structure-preserving numerical methods for such systems using Galerkin Hamiltonian variational integrators which admit discrete analogues of these quadratic conservation laws. We also consider the generalization to flows on Lie groups by studying Lie group variational integrators with a novel Type II variational principle on the cotangent bundle of a Lie group which allows for Type II boundary conditions. Finally, we investigate the geometric structure of adjoint systems associated with evolutionary partial differential equations at the fully continuous, semi-discrete, and fully discrete levels and the relations between these levels.

ABSTRACT. In this talk, I will present some recent results on quantum algorithms for differential equations. First, I will present some lower bounds for quantum algorithms to solve certain fluid dynamical equations like Euler equation and the Korteweg-de Vries equation. For the Euler equation, one can show that certain instabilities such as the Kelvin-Helmholz instability will cause the solutions to move far from the initial states. This, in turn, means that the quantum algorithm will take time exponential in the simulation time to solve the differential equation. Next, I will present some recent work on quantum algorithms for differential equations with an algebraic component as well as applications of these equations.

ABSTRACT. Simulating the time evolution of quantum systems remains one of the most promising applications of quantum computing. In this talk, we will present an efficient quantum algorithm designed to simulate slowly varying time-dependent Hamiltonians. By leveraging Floquet theory alongside a smooth extension of the Hamiltonians to periodic systems, our approach achieves near-optimal scaling, specifically, an almost linear and additive dependence on evolution time and error parameters. We will also discuss how to extend this algorithm to general slow non-unitary dynamics using the linear combination of Hamiltonian simulation technique.

ABSTRACT. We present an efficient, nearly optimal quantum algorithm for solving linear matrix differential equations, with applications to the simulation of open quantum systems and beyond. For unitary or dissipative dynamics, the algorithm computes an entry of the solution matrix with query complexity \widetilde{O}(L^2 /\epsilon), where L involves the time integral of the norm of the evolution operator and epsilon is the error. This L is linear in the evolution time t for unitary dynamics and can be bounded by a constant for dissipative dynamics. Our result contrasts prior quantum approaches for differential equations that typically require exponential time for this problem due to the encoding in a quantum state, which can lead to exponentially small amplitudes. We demonstrate the utility of the algorithm through an end-to-end application, namely the simulation of quantum Langevin dynamics for non-interacting fermions, and comment on applications to other quantum and classical systems. We compare with classical algorithms and give evidence of significant polynomial quantum speedups for systems in a lattice, which become more pronounced for systems with long-range interactions, and are exponential in general. We also provide a lower bound of \Omega(L^2/\epsilon) for unitary or weakly dissipative dynamics that proves our algorithm is optimal up to logarithmic factors.

ABSTRACT. We introduce \emph{random-LCHS}, a circuit-efficient randomized-compilation framework for simulating linear non-unitary dynamics of the form $\partial_t u(t) = -A(t)\,u(t) + b(t)$ built on the linear combination of Hamiltonian simulation (LCHS). We propose three related settings: the \emph{general random-LCHS} for time-dependent inhomogeneous linear dynamics; the \emph{observable-driven random-LCHS}, which targets estimation of an observable’s expectation at the final time; and the \emph{symmetric random-LCHS}, a time-independent, homogeneous reduction that can exploit physical symmetries. Our contributions are threefold: first, by randomizing the outer linear-combination-of-unitaries (LCU) layer as well as the deterministic inner Hamiltonian simulation layer, random-LCHS attains favorable resource overheads in the circuit design for early fault-tolerant devices; second, the observable-driven variant employs an unbiased Monte-Carlo estimator to target expectation values directly, reducing sample complexity; and third, integrating the physical symmetry in the model with the sampling scheme yields further empirical improvements, demonstrating tighter error bounds in realistic numerics. We illustrate these techniques with theoretical guarantees as well as numerical verifications and discuss implementation trade-offs for near-term quantum hardware.Quantum algorithms & computation

ABSTRACT. Spatial self-organization is widely recognized as a mechanism that enhances ecosystem resilience, yet most studies have focused on stationary patterns arising from scale-dependent feedbacks. Here, we investigate a different route to resilience driven by excitable dynamics in systems with delayed negative feedbacks. Using a minimal reaction–diffusion model coupling plant biomass and toxin accumulation, we show that ecosystems can avoid collapse through the emergence of dynamic spatiotemporal structures, including defect turbulence, spiral waves, wave trains, and traveling pulses. As environmental stress increases, the system undergoes a predictable sequence of regimes, culminating in isolated traveling pulses that persist well beyond the tipping point of the corresponding homogeneous system. These moving structures continuously escape locally accumulated toxicity, effectively redistributing stress in space and allowing positive feedback to dominate at their leading edge. As a result, the ecosystem maintains higher biomass and survives under conditions where non-spatial models predict extinction. The framework captures observed patterns in seagrass meadows and applies broadly to plant–soil systems with mediated negative feedbacks. Our results highlight the importance of nonequilibrium dynamics and excitability in ecological resilience, providing new indicators of regime shifts and emphasizing the role of transient, mobile structures in sustaining ecosystems under global change.

ABSTRACT. In recent years, the development and analysis of advection-diffusion-reaction models have gained increasing attention due to their applicability in various sustainability-related contexts. Notable examples include the modeling of ecological dynamics, such as vegetation patterns in arid and semi-arid ecosystems [1]. A further relevant application, in the framework of sustainable innovation, arises in the mathematical modeling of supply chains and production networks, where PDE-based continuum descriptions can be used to represent the evolution of material and product flows across interconnected processors, buffers and distribution nodes [6]. These models are typically described by systems of two-dimensional Partial Differential Equations (PDEs). Solving such problems often requires numerical integration over long time intervals combined with very fine spatial discretizations. This becomes especially demanding in tasks such as parameters estimation, where the model is calibrated to reproduce observed data. Indeed, this generally involves solving the PDEs system multiple times within optimization algorithms. Therefore, the use of efficient numerical methods, capable of providing a stable and accurate solution in short computing times, is crucial. Building on recent advances in linearly implicit W-methods [2], in this talk we focus on the construction and analysis of schemes coupled with splitting strategies based on Approximate Matrix Factorization (AMF) [3], and with Matrix-Oriented (MO) techniques [4], specifically designed for multidimensional advection-diffusion-reaction problems. These approaches exploit the structure of the discretized operators, leading to a reduction in computational costs, while preserving the consistency of the underlying W-method. We analyze the properties of the new methods in terms of accuracy, stability and computational cost [5]. Numerical experiments show the effectiveness of the proposed approaches, including in parameters estimation tasks for vegetation PDEs models calibrated on real data.

Acknowledgments: this research falls within the activities of PRIN-MUR 2022 project 20229P2HEA Stochastic numerical modelling for sustainable innovation, CUP E53D23017940001, granted by the Italian Ministry of University and Research (call relating to scrolling of the final rankings of the PRIN 2022).

Bibliography [1] M. Abbas, F. Giannino, A. Iuorio, Z. Ahmad, F. Calabrò, PDE models for vegetation biomass and autotoxicity. Math. Comput. Simul., 228, 386–401 (2025). [2] D. Conte, J. Martin-Vaquero, G. Pagano, B. Paternoster, Stability theory of TASE-Runge-Kutta methods with inexact Jacobian. SIAM J. Sci. Comput. 46(6), A3638–A3657 (2024). [3] D. Conte, S. Gonzalez-Pinto, D. Hernandez-Abreu, G. Pagano, On Approximate Matrix Factorization and TASE W-methods for the time integration of parabolic PDEs. J. Sci. Comput., 100, 34 (2024). [4] M. C. D’Autilia, I. Sgura, V. Simoncini, Matrix-oriented discretization methods for reaction–diffusion PDEs: Comparisons and applications. Comput. Math. Appl., 79(7): 2067-2085 (2020). [5] D. Conte, S. Iscaro, G. Pagano, On Matrix-Oriented and AMF W-methods for advection-diffusion-reaction problems. In preparation. [6] M. Herty, A. Klar, B. Piccoli, Existence of solutions for supply chain models based on partial differential equations. SIAM J. Math. Anal., 39(1), 160–173 (2007).

ABSTRACT. We investigate the preservation of diffusion-driven (Turing) instability under first-order matrix-oriented splitting integrators for two-species reaction-diffusion systems on 2D domains. Starting from tensor-product spatial discretizations, the semi-discrete problem is reformulated as a differential matrix equation, yielding matrix-oriented splitting schemes in which diffusion and reaction processes are treated separately. Within this framework, we derive fully discrete modal amplification matrices and characterize the onset of instability through the associated Jury conditions. The analysis reveals a clear separation between the continuous Turing mechanism, inherited from the underlying reaction-diffusion model, and purely discrete effects introduced by the time discretization. Focusing on the Gierer-Meinhardt system as a benchmark problem, we show that different splitting strategies may substantially alter the continuous Turing instability region: some schemes generate stable spurious numerical patterns in parameter regimes where the continuous model is Turing stable, whereas others may suppress pattern-forming instabilities. For the IMEX scheme, we further derive a time-step condition ensuring preservation of the continuous Turing region. Numerical experiments confirm the theoretical predictions and illustrate how each integrator induces its own discrete Turing region, emphasizing the need to assess the preservation of pattern-forming properties when interpreting numerical solutions.

ABSTRACT. We talk about the development and analysis of a higher-order unfitted finite element method for the incompressible Stokes equations, which yields a strongly divergence-free velocity field up to the physical boundary. The method combines an isoparametric Scott--Vogelius velocity-pressure pair on a cut background mesh with a stabilized Nitsche/Lagrange multiplier formulation for imposing Dirichlet boundary conditions. We construct finite element spaces that admit robust numerical implementation using standard elementwise polynomial mappings and produce exactly divergence-free discrete velocities.

ABSTRACT. Reaction-diffusion equations on surfaces are widely used for modeling various phenomena in biology. In this talk, we will introduce recent progress on solving a surface reaction-diffusion system coupled with the evolution of the surface. In particular, we will present the application of a stabilized trace finite element method for the reaction-diffusion system on evolving surfaces. The surface motion is computed by a level-set method. Both the trace finite element space for the reaction-diffusion system and the finite element space for the level-set function are defined in a narrow band region near the surface on a bulk mesh. The method is fully decoupled and allows for easy handling of topology changes.

ABSTRACT. We propose an energy-stable inertial level-set flow for solving PDE-constrained topology optimization problems. This flow integrates the second-order time derivative of the level-set function into the normal velocity mimic an inertial dynamic, thus is capable of escaping from local minima. Unconditional energy stability is rigorously established for the continuous-in-time dynamical system, and its corresponding nonlinear temporal and fully discretized scheme. The existence and uniqueness of the level-set function update at each time step of the scheme are proven via Banach’s fixed point theorem. Furthermore, we present a convergence analysis for the long-time asymptotic behavior of the temporal discretization scheme. Numerical examples verify that the proposed level‑set flow features broad applicability and yields highly competitive performance in both structural and fluid topology optimization problems.

ABSTRACT. In this talk, we introduce a unified neural flow framework that provides an infinite-depth formulation for deep neural networks and operators. Two representative dynamical systems recover plain and ResNet-type architectures through time discretization. We establish well-posedness and develop approximation theory for both networks and operators. The framework also incorporates various spatial discretizations for inter-neuron linear operators, enabling coverage of existing neural operator architectures and yielding approximation results for finite-depth DNNs, CNNs, and neural operators within a single continuous perspective.

ABSTRACT. A lot of physical phenomea are modelled by Hamiltonian systems, which are entirely determinated by a real function called the Hamiltonian, which represents the energy of the system. One of the major property of the flow is that it is a symplectomorphism, i.e. an application which preserves a differential 2-form called the symplectic form or symplectic structure of the phase space over time. When numerically solving Hamiltonian equations, taking this property into account using methods that preserve the symplectic structure ensures better stability over long time and physical consistency. When using a neural network to solve a Hamiltonian equation while preserving the symplectic structure, a specific architecture must be considered, using symplectic layers. A well-known symplectic architecture used to learn Hamiltonian ODEs was proposed in [1]. Here, we propose a way to extend this framework to Hamiltonian PDEs and propose an architecture of symplectic neural operators. We compare the performance with non-symplectic neural operators for wave equation in a heterogeneous medium.

ABSTRACT. We propose a physics-informed neural network framework, termed scaled TW-PINN, for computing traveling wave solutions of n-dimensional reaction-diffusion equations with various reaction and diffusion coefficients. By applying a scaling transformation based on the traveling wave form, we reduce the original problem to a one-dimensional scaled reaction-diffusion equation with unit reaction and diffusion coefficients. This reduction enables the construction of a single PINN solver that can be reused across different coefficient regimes and spatial dimensions. We prove a universal approximation property of the proposed solver for traveling wave solutions. Numerical experiments in one and two spatial dimensions, together with comparisons to the wave-PINN method, demonstrate the accuracy and computational advantages of the proposed approach. Finally, we extend the proposed framework to the Fisher’s equation with general initial conditions and present corresponding numerical results.

ABSTRACT. Neural operators have gained much attention for accelerating physics simulations. However, they often suffer from capturing the laws of physics from data. This study proposes a structure-preserving operator learning method for non-canonical Hamiltonian partial differential equations (PDEs). Under periodic boundary conditions, the equations are formulated as finite-dimensional Hamiltonian systems via truncated Fourier spectral representations, and their time evolution is characterized as symplectic mappings. By modeling these mappings with neural networks, the proposed approach enables fast data-driven simulations while preserving the laws of physics such as energy conservation. We also provide experimental results for some non-canonical Hamiltonian PDEs.

ABSTRACT. In this talk, we present a model for solving mixed-dimensional problems with highly heterogeneous coefficients, a type of problem that commonly appears in e.g. modelling of fractured porous media but can be computationally challenging to solve numerically. Thin structures are often modelled as lower-dimensional interfaces embedded in a higher-dimensional bulk domain, leading to the mixed-dimensional model problem.

Our method is based on the Localised Orthogonal Decomposition (LOD) method and constructs locally supported basis functions on a coarse mesh that does not resolve the fine-scale variations of the coefficients. The basis functions are adapted to the problem at hand and thus carry the fine-scale information in order to ensure optimal convergence with respect to the coarse mesh, independent of the coefficient regularity. This method leads to an exponentially decaying localisation error. We present numerical experiments to validate the theoretical findings and demonstrate the computational viability of the method.

ABSTRACT. Low-rank approximation techniques and domain decomposition methods have both attracted considerable attention in recent years as promising approaches for tackling high-dimensional problems. In this work, we combine these two ideas and introduce a domain decomposition dynamical low-rank method for the efficient solution of high-dimensional radiative transfer and similar kinetic equations. The algorithm uses a separate low-rank approximation on each spatial subdomain, which means that, for a given accuracy, we can often use a smaller overall rank compared to classic dynamical low-rank methods. Additionally, our algorithm only transfers boundary data between subdomains and is thus very attractive for distributed memory parallelization, where classic dynamical low-rank algorithms suffer from global data dependency. In this talk, we will first briefly recap the classic dynamical low-rank approach and then demonstrate how it can be naturally combined with domain decomposition, presenting the resulting algorithm and its efficiency on challenging test cases featuring both very optically thin and thick regions.

ABSTRACT. The solution to a differential Lyapunov equation can be expressed in closed form as a matrix-valued integral, the so-called finite-horizon Gramian. Such Gramians also have applications in many other areas, such as optimal control and Gauss-Markov regression. The Gramian is positive semi-definite, and often it is more useful to have a Cholesky factorization of it rather than the Gramian itself. This work considers a new efficient numerical method for computing such Cholesky factors of finite-horizon Gramians without first computing the full Gramian. There is currently no alternative general-purpose method for this task; all plausible combinations of methods known to us either break down when applied to many interesting problems or are much less efficient. The proposed method is a generalization of the general-purpose scaling-and-squaring approach for approximating the matrix exponential. It exploits a similar doubling formula for the Gramian, and thereby keeps the required computational effort modest. Most importantly, we have performed a rigorous backward error analysis that guarantees that the approximation is accurate to the round-off error level in double precision if the method parameters are chosen appropriately. This is complemented by a forward analysis based on estimating the condition numbers of computing the Gramian. In the talk I will describe the method, outline the main backward result, and show some experimental results produced by our efficient and freely available Julia code.

ABSTRACT. Partial differential equations (PDEs) model a wide range of physical phenomena, including electromagnetic fields and fluid mechanics. Methods such as sparse identification of nonlinear dynamics (SINDy) and PDE-Net 2.0 identify and model PDEs from data via sparse optimization and deep neural networks, respectively. While PDE models are less commonly applied to fMRI data, they can uncover hidden connections and essential components of brain activity. Using the ADHD200 dataset, we applied canonical independent component analysis (CanICA) and uniform manifold approximation (UMAP) for dimensionality reduction of fMRI data. We then used sparse ridge regression to identify PDEs from the reduced data. We applied advanced PDE features for classification, achieving high accuracy in distinguishing individuals with attention-deficit/hyperactivity disorder (ADHD). This study demonstrates a novel approach to extracting meaningful features from fMRI data for the analysis of neurological disorders, to understand the role of oxygen transport (delivery & consumption) in the brain during neural activity, which is relevant for studying intracranial pathologies.

ABSTRACT. This paper investigates the consensus problem for multi-agent systems modeled by reaction-diffusion partial differential equations. Considering the practical constraints that complete state information is often unavailable and that control inputs can only be applied at boundary of spatial domain, a boundary control strategy based on partial spatial-domain measurements is proposed. First, to address the issue of unmeasurable states, an observer is designed using measurement information from piecewise spatial subdomains or discrete points. Based on the observed states, a boundary control protocol is developed to achieve consensus. By employing the Lyapunov direct method and Poincaré-Wirtinger inequality, sufficient conditions for the asymptotic stability of the closed-loop error system are derived and formulated as Linear Matrix Inequalities. Finally, the effectiveness of the proposed method is validated through numerical simulations.

ABSTRACT. This paper studies a class of generalized nonlinear functional integral equations that arise in models of infectious disease transmission. The problem is formulated in the Banach space C[0; 1]. Using fixed point theory, we establish sufficient conditions for the existence and uniqueness of solutions. These results ensure that the model is well posed and provide a unified setting for several related nonlinear integral equations.

We also develop a hybrid iterative method that combines a modified Homotopy Perturbation technique with the Adomian Decomposition Method to compute approximate solutions. The convergence of the method is proved rigorously, and explicit error estimates are derived to measure its accuracy and stability. Numerical examples support the theoretical results and show that the method is accurate and computationally efficient for solving nonlinear functional integral equations arising in epidemiological modeling.

Finally, possible directions for future research and further extension of the proposed approach are briefy discussed.