Papers

This paper proposes an efficient and scalable preconditioning strategy for eddy current problems involving coupled external circuits. The approach, named Electro-Magnetic Decoupling (EMD) preconditioner, decomposes the discrete system into vector and scalar potential components and applies tailored preconditioners to each. In particular, strong preconditioning is applied to the scalar component to address the spectral degradation induced by the discrete Laplacian. The method was evaluated on four eddy current models with varying frequencies, conductor topologies, and excitation types. Compared to the conventional incomplete Cholesky preconditioner, the EMD approach achieved up to 20 times fewer iteration counts and up to 20 times faster iterative solver time. Moreover, the method supports physics-level parallelization, allowing efficient treatment of independently excited conductor domains. The EMD framework is compatible with algebraic multigrid, domain decomposition, and direct solvers, offering flexibility and robustness for large-scale electromagnetic simulations.

In this work, we present a novel dynamic proximal point algorithm for unconstrained optimization. The method generates a sequence of proximal subproblems, where the quadratic regularization term is weighted by a diagonal matrix that is updated adaptively at each iteration. Each subproblem is solved using an inner Newton's method combined with a line search, which provides a global convergence mechanism for the nonlinear solver. At the outer level, the algorithm updates the reference point and adjusts the regularization parameter based on the performance of the inner Newton solver. We derive the reduced linear system used to compute the Newton step, define the corresponding merit function, and discuss practical approaches for constructing the diagonal scaling matrix from derivative information. The paper also provides implementation-oriented pseudocode and stopping criteria that are consistent with the proposed method.

A fully implicit two-stage fourth-order two-derivative time discretization was introduced previously as a temporal method. This paper closes that sequential integrator for stiff transport–relaxation equations by pairing a conservative finite-volume residual $\mathcal L_h$ with its discrete trajectory derivative $\mathcal G_h^{\rm tr}=D\mathcal L_h\,\mathcal L_h$. An ADER/Cauchy–Kowalevski predictor provides interface states and physical time derivatives; differentiating the same numerical flux and taking shared face differences yields a conservative approximation $\widetilde{\mathcal G}_h$. For linear constant-coefficient balance laws, $\widetilde{\mathcal G}_h=\mathcal G_h^{\rm tr}=\mathcal L_h^2$ exactly, although the derivative operator is assembled independently rather than by squaring the residual matrix. For nonlinear discretizations, the fourth-order temporal theory applies to $\mathcal G_h^{\rm tr}$, while a trajectory-closure consistency estimate controls the ADER approximation. The two unknown stage vectors are solved successively through two $N$-unknown systems. The completed step is fourth order and L-stable; the parameter $C_q=5/183$ cancels the leading inverse-power term and changes the deep-stiff amplification from $O(|z|^{-1})$ to $O(|z|^{-2})$. For fixed compatible spatial spaces, a slow–fast decomposition proves a full-step asymptotic-preserving operator limit with an $O(\delta)$ estimate and gives a preparation-dependent uniform-accuracy classification. Linear finite-volume, nonlinear relaxation, one- and two-dimensional damping, diffusion-limit, and modal experiments verify the corresponding closure, accuracy, stability, and singular-limit claims within their stated scopes.

High-dimensional sensitivity analysis and uncertainty quantification for multiscale gas dynamics, spanning the continuum to rarefied regimes, require computationally efficient and mathematically consistent gradient evaluation. This paper develops a discrete adjoint method for the unified gas-kinetic scheme (UGKS) based on a dual-consistent formulation. The adjoint system is derived directly from the discrete microscopic velocity-distribution equation coupled with the macroscopic-moment compatibility conditions. To resolve the stiff cross-scale coupling, we propose an asymptotic-preserving (AP) adjoint formulation constructed via macroscopic-moment projection and microscopic lifting. Under this framework, the AP adjoint formulation eliminates the stiff collision coupling and removes the collision-time step restriction in the continuum regime. Numerically, a memory-efficient residual-evaluation algorithm that mirrors the forward UGKS cell-vertex data structure is implemented to bypass the memory bottleneck in velocity space. Furthermore, a macroscopic–microscopic predictor–corrector implicit marching scheme is designed to accelerate convergence without solving a globally coupled system. The accuracy, consistency, and robustness of the proposed AP-adjoint scheme are rigorously verified against an independent linearized UGKS solver across a wide range of Knudsen numbers, including lid-driven cavity heat conduction, microchannel thermal creep flow, and hypersonic flow past a circular cylinder.

We study the numerical evaluation of an energy-based vector hysteresis model and its incorporation into finite element simulations based on a vector potential formulation. The inherent non-smoothness of the hysteresis model poses challenges for both numerical analysis and implementation. Using tools from convex analysis, we characterize the forward and inverse hysteresis operators in terms of energy densities. This characterization yields well-posedness of the resulting models and leads to robust algorithms for their evaluation based on generalized semi-smooth Newton methods. It furthermore enables a seamless integration into magnetic field simulations, leading to nonlinear and non-smooth optimization problems at every load step. We discuss the finite element discretization, present a semi-smooth Newton method for the iterative solution, and establish global linear convergence with mesh-independent convergence rates. The theoretical results are illustrated by numerical experiments.

We propose a gravity-consistent arbitrary Lagrangian–Eulerian finite element method for incompressible Navier–Stokes free-boundary problems with moving contact lines. A direct body-force discretization of gravity may fail to ensure consistency between the discrete gravitational work and the variation of the gravitational potential energy on the evolving domain, resulting in an artificial consistency error and persistent spurious velocities near equilibrium. To remove this inconsistency, we reformulate the gravitational potential energy variation as a moving-boundary integral over intermediate ALE configurations and evaluate it exactly using Simpson's quadrature rule. This leads to a mildly nonlinear fully discrete scheme in which the gravitational contribution is exactly consistent with the discrete potential-energy variation. The proposed method preserves volume exactly, satisfies a discrete energy-dissipation law including gravitational potential energy, and under suitable assumptions, drives the discrete velocity to zero in the long-time regime, thereby excluding persistent gravity-induced spurious velocities. Together with the BGN treatment of the free surface and the MDR bulk mesh extension, the scheme maintains accurate interface tracking and good mesh quality near the moving contact line. Numerical experiments in two and three spatial dimensions confirm the theoretical properties.

This paper presents a fully discrete analysis of pressure-robust enriched Galerkin (EG) methods for the time-dependent Stokes equations. The spatial discretization employs a velocity reconstruction operator that preserves the continuous component of the EG space while mapping its enriched discontinuous component into the lowest-order Raviart-Thomas ($\mathcal{R}T_0$) subspace of $H(\mathrm{div})$, ensuring strict local mass conservation and pressure-robustness. Fully discrete schemes are constructed using the backward Euler and Crank-Nicolson methods, with the reconstruction incorporated into both the forcing and discrete time-derivative terms. We prove unconditional stability and derive optimal-order, parameter-explicit a priori error estimates for both the velocity and pressure. In particular, the velocity estimates are independent of the continuous pressure and the irrotational component of the forcing term, contain no inverse-viscosity factors, and explicitly track the spatial mesh size, time-step size, and viscosity. Numerical experiments in two and three dimensions verify the predicted convergence rates, strict local mass conservation, and robustness with respect to the viscosity, while comparisons with non-pressure-robust variants demonstrate the importance of incorporating the reconstruction into the discrete time-derivative term.

We numerically investigate selection criteria for dissipative weak solutions of the compressible Euler equations. Using the Kelvin-Helmholtz problem as a benchmark, we compare four standard numerical schemes - the viscous finite volume (VFV) method, the discrete velocity Boltzmann finite volume (DVBFV) method, the Runge-Kutta discontinuous Galerkin (RKDG) method, and the discontinuous Galerkin spectral element method (DGSEM), and demonstrate that each may converge to a different dissipative weak solution. We evaluate these solutions with respect to several selection criteria based on entropy production, total energy, and energy defect, and examine the sensitivity of selection functionals to the numerical diffusion inherent in each scheme.

Forsythe conjectured in 1968 that the normalized residual directions of restarted conjugate gradients approach a two-cycle. We prove the conjecture for restart length two: unless the method terminates, its even and odd normalized residuals converge separately.

We introduce a new SVD-based tensor decomposition method for tensor networks with arbitrary graph topologies, extending classical hierarchical SVD-based techniques to networks with cycles and general connectivity. We also introduce addition and rounding procedures for traceable tensor graphs, enabling step-rounding time integration of high-dimensional PDEs directly in graph format, with rank truncation controlled to a prescribed tolerance at every time step. We demonstrate the new method on the decomposition of multivariate functions and on the numerical solution of the Fokker-Planck equation, and find that the graph-format representation attains comparable or better accuracy than the classical tensor train and hierarchical Tucker tensor formats, while using substantially fewer degrees of freedom at lower computational cost.

Multiphysics initial-value problems couple processes with distinct stability properties, such as advection, diffusion, and stiff local reactions. Standard implicit-explicit (ImEx) additive Runge–Kutta (ARK) methods can treat these processes accurately, but require globally coupled implicit solves when diffusion is grouped with reaction; operator splitting avoids such solves but typically provides weaker coupling and no inexpensive temporal error estimate; and PIROCK is tied to a specific Runge–Kutta–Chebyshev super-time-stepping (STS) construction. We introduce extended super-time-stepping (ExtSTS) methods, a family of time integration schemes that combine super-time-stepping methods for diffusive terms with explicit, implicit, or ImEx Runge–Kutta treatment of the remaining terms. The coupling is based on multirate infinitesimal techniques, yielding solve-decoupled methods that retain localized implicit solves, support embedded error estimation for adaptive time stepping, and allow flexible use of modern STS methods. We present the ExtSTS method family, provide a robust technique for ExtSTS method creation, formulate the corresponding linear stability theory, and construct embedded ImEx, explicit, and implicit ExtSTS methods. Numerical experiments on one- and two-dimensional advection-diffusion-reaction problems show that ExtSTS methods are robust across parameter regimes and boundary conditions, and are often more efficient than ARK, Strang splitting, and PIROCK methods, especially when strong coupling between operators is important.

Polymer quantization program applied in Loop Quantum Gravity/Cosmology leads to nonseparable Hilbert spaces. Commonly, one sidesteps this problem by singling out and working with a separable superselection sector. This is however often no longer accessible in more involved models beyond isotropic ones. In an alternative approach one builds a separable Hilbert space as an integral over all available sectors. Here we test the dynamics following from the latter on the example of a flat isotropic Universe admitting negative cosmological constant and a massless scalar field. There, numerical evolution of (initially) semiclassical states shows that there is no relevant difference in their long term semiclassicality properties in comparison to those in a single sector approach. Further, the older (problem-specific) numerical methods are compared against an application of more common and more efficient standard tools (eigen library).

The Poisson-Boltzmann (PB) theory is a cornerstone of implicit solvent models for electrostatic analysis, and has found a great success in various biomolecular applications. However, in calculating polar solvation energy, one should consider that the structure of the protein changes upon transition from vacuum to water phases. To address this, here we report for the first time a generalized PB framework capable of accommodating nonrigid conformational changes without suffering from self-energy artifacts. For regularized PB models, in which the charge singularities are captured by the Green's functions, self-energies in the water and vacuum states will be analytically canceled. For non-regularized PB solvers, such as APBS and DelPhi, a simple thermodynamic cycle is proposed for nonrigid proteins by adding a Coulombic correction in vacuum. The generalized PB theory is validated using a perturbed two-atom system and a diverse set of proteins with different structures in vacuum and water, demonstrating its accuracy and robustness, regardless of the choice of sharp-interface and diffuse-interface PB models and different numerical solvers.

We provide new estimates in Wasserstein distance for the asymptotic bias of the unadjusted Langevin algorithm, in the classical setting of log-smooth strongly log-concave measures. Our bound implies a Wasserstein mixing time of order $\kappa \sqrt{d}/\varepsilon$, where $\kappa$ is the condition number, $d$ is the dimension, and $\varepsilon$ is the target precision: this improves by a factor of $\sqrt{d}/\varepsilon$ over the previous state-of-the-art results.

When the spectrum of a system varies significantly over time, fixed choices of macro-, meso-, and micro-time steps, as well as burst lengths, become inadequate, necessitating adaptive and locally informed strategies. To address these challenges, this article proposes a novel and flexible generalised projective integration (GPI) scheme, designed to accommodate time-dependent spectral variation and dynamically evolving scale separation. The proposed framework unifies and extends several existing multiscale methodologies, thereby offering a more general and adaptable computational paradigm. A comprehensive stability analysis of the GPI scheme is carried out, including a detailed investigation of the splitting of the stability region, which forms a central component of this work. Furthermore, problem-dependent strategies for selecting the micro-, meso-, and macro-time steps, as well as the burst length, are developed and their impacts are systematically validated through numerical experiments. To assess the effectiveness of the proposed scheme, three representative problems with distinct types of spectral evolution are considered. The performance of the proposed GPI scheme is evaluated and compared with several existing projective integration methods as well as some widely used stiff solvers, based on (i) number of micro time steps, (ii) accuracy, (iii) computational time, (iv) memory usage and (v) proportion of micro-scale simulations. The results demonstrate that the GPI scheme consistently outperforms the existing methods in terms of performance.

We introduce a unified framework for the error analysis of generative models based on the entropy production rate of the forward-reverse diffusion process pair. For a pair of continuity equation flows, the rate admits a closed velocity form identity whose time integral decomposes the terminal Kullback–Leibler (KL) divergence into the sum of an initialization error, a score approximation error, and a time-discretization error. By analyzing the entropy production at the level of marginal distributions, rather than in path space, our framework yields a sharp convergence rate of $\mathcal{O}(h^2)$ for the Euler-Maruyama sampler, where $h$ is the step size. This improves upon the $\mathcal{O}(h)$ rates typically obtained from Girsanov's path-space analyses. Furthermore, our framework unifies the analysis of score-based SDEs, probability-flow ODEs, and stochastic interpolants by varying diffusion coefficients within a single inequality, revealing the trade-off between deterministic and stochastic sampling. Numerical experiments confirm the predicted scaling with step size and terminal time.

We develop a least-squares weak Galerkin (LS-WG) finite element method for the Cauchy problem of the biharmonic equation. The proposed approach reformulates the fourth-order equation as a coupled system of two second-order equations, which are discretized using discrete weak Laplacian operators on weak finite element spaces. The resulting least-squares formulation yields a symmetric positive definite linear system, thereby eliminating the discrete inf-sup condition required by mixed finite element methods while avoiding the construction of globally $C^1$-conforming finite element spaces. Furthermore, the weak Galerkin framework naturally accommodates general polygonal meshes, offering considerable flexibility in mesh generation and approximation. Under the assumption that the continuous biharmonic Cauchy problem admits a unique solution, we establish the uniqueness of the discrete LS-WG solution and derive optimal-order error estimates in a discrete energy norm. Numerical experiments confirm the theoretical convergence rates and demonstrate the accuracy, robustness, and effectiveness of the proposed method.

We introduce a mathematical framework for shape comparison based on mapping functions from the shape domain to a common reference domain. This Push-Forward Transform enables invariant and robust comparison of shapes, preserving intrinsic geometric information. Quantitatively comparing shapes and their temporal evolution is a fundamental challenge in image analysis. Meaningful shape comparison requires representations that are invariant to transformations that do not alter shape itself, such as translation, rotation, reflection, re-parametrization, and uniform scaling, while remaining sensitive to intrinsic geometric variation. Existing approaches often rely on sensitive parameterizations, landmark correspondence, or learned representations that are difficult to interpret and reproduce. We show that the Push-Forward Transform (PF-T) applied to Signed Distance Functions (SDFs) yields a continuous representation that captures both boundary and interior geometry. We derive an interpretable morphometric that quantifies shape similarity and reveals features such as skeletal topology and rotational symmetries. The push-forward transform applies consistently to two- and three-dimensional shapes, extends to time-evolving geometries, and supports the joint analysis of shape and additional scalar fields defined over shapes, such as intensity or molecular signals. We present the mathematical formulation, describe an efficient algorithm, and benchmark the approach on 2D, 3D, and temporal data sets.

Shock stabilization in compressible Euler flows remains a central challenge for high-order numerical methods. Existing shock-capturing approaches, including limiters, artificial viscosity, and reconstruction-based methods, involve tradeoffs between robustness, accuracy, preservation of fine-scale flow features, and computational complexity. In this work, we develop a discontinuous Galerkin (DG) discretization of the information geometric regularization (IGR) framework introduced by Cao and Schäfer for the compressible Euler equations. The method regularizes shocks at the PDE level by augmenting the Euler equations with the entropic pressure $\Sigma$, obtained from an auxiliary elliptic equation. Within the DG formulation, the regularization enters only through the augmented pressure $P+\Sigma$ in the Euler fluxes, preserving the conservative structure of the discretization while using a common approximation space for both the hyperbolic and elliptic equations. Numerical experiments spanning one and two-dimensional benchmark problems show the proposed formulation stabilizes shocks without shock-capturing limiters or artificial viscosity, although positivity-preserving methods may still be required when the density or pressure approaches zero. Compared with a characteristic TVB-limited DG formulation, the IGR-DG method resolves increasingly finer-scale flow features as the polynomial order is increased while maintaining stable shock resolution. The entropic pressure remains localized to regions of strong compression with minimal activation in smooth regions of the flow, providing selective PDE-level regularization while preserving the underlying solution elsewhere.

This paper establishes the first localized pointwise a posteriori error estimates for nonconforming finite element discretizations of the Poisson and biharmonic equations. For the Poisson problem, we derive localized estimates for the function-value and broken gradient errors of the Crouzeix–Raviart method. For the Morley discretization of the biharmonic equation, we derive a localized a posteriori estimate that controls the local Hessian error. This provides the first pointwise a posteriori error analysis for the biharmonic equation, for either conforming or nonconforming finite element methods. The key ingredient, absent from pointwise analysis of second order PDEs, is the design of two novel weight functions that facilitate sharp estimates of the $L^1$ norms of derivatives of a regularized Green's function for the biharmonic operator.

We propose a finite element method for Yang–Mills connections on nontrivial principal bundles with abelian structure group. Local connection forms satisfy internal jump conditions induced by the transition functions of the bundle. We discretize these forms in broken finite element exterior calculus spaces and enforce the jump conditions weakly by Lagrange multipliers. We prove well-posedness of the resulting saddle-point problem and derive an a priori convergence estimate. A numerical experiment for the Hopf fibration illustrates the method and exhibits linear convergence when the sphere is approximated by a piecewise linear mesh.

The tangent plane scheme (TPS) is a well-established discretization of the Landau-Lifshitz-Gilbert (LLG) equation. However, rigorous a posteriori error estimates have not been established. In this work, we derive a rigorous a posteriori error estimate for the TPS based on the second-order backward differentiation formula (BDF(2)) in time and finite elements of arbitrary polynomial degree in space. The proposed estimators provide computable upper bounds for the temporal and spatial discretization errors on adaptive meshes with variable time-step sizes. This result establishes the mathematical foundation for fully adaptive algorithms for the LLG equation.

We propose a novel UnCut finite element method (FEM) for the Poisson and Stokes equations on domains with curved boundaries represented by a level set function. Like the $\phi$-FEM, the method avoids numerical integration over cut subregions of boundary elements, but introduces a novel least-squares formulation that minimizes an $H^{-1}$ residual of the governing equations. This formulation ensures stability without requiring large stabilization parameters, thereby eliminating the need for user-tuned penalty parameters and improving the robustness of the computation. Optimal-order convergence of the UnCut FEM solutions is rigorously established in the $H^1$ norm for both the Poisson and Stokes equations, and numerical experiments are presented to support the theoretical analysis.

We study fractional diffusion equations on compact metric graphs, where the nonlocal dynamics is governed by fractional powers of the shifted Kirchhoff-Laplacian. Building on recent advances in the analysis of fractional operators on metric graphs, we establish a rigorous mathematical framework and propose a fully discrete scheme based on backward Euler time-stepping and finite element discretization. To approximate the action of the fractional operator, we employ rational approximations, reducing the problem to a sequence of sparse elliptic solves for efficient implementation. We derive error estimates for the temporal, spatial, and rational discretizations, and confirm convergence through numerical experiments.

Transcranial focused ultrasound (tFUS) requires accurate estimation of the intracranial acoustic field, which is distorted by skull-induced aberrations. Numerical solvers are accurate but computationally expensive for digital twins, where the field must be re-estimated repeatedly as treatment conditions change. Existing deep-learning surrogates are fast but typically use voxel-to-voxel regression on a fixed grid, with no mechanism reflecting how acoustic energy propagates through the skull. We instead cast tFUS simulation as an operator learning problem and propose tFUSOperator, a coordinate-aware neural operator that maps the free-field pressure, skull anatomy, and treatment parameters to the intracranial field within a shared physical coordinate frame. To our knowledge, this is the first operator-based formulation of tFUS field prediction. On both seen and unseen skulls, the model localizes the acoustic focus accurately-reaching about 90% and 72% Dice, respectively-and it performs nearly as well from magnetic resonance (MR) as from computed tomography (CT) input while running $5.6 \times 10^4$ times faster than numerical simulation. These results suggest a fast, radiation-free route to safe and practical digital twins for patient-specific tFUS treatment. The code is available at: https://github.com/CMME-Lab/tFUSOperator.git.

The regularity of solutions to elliptic partial differential equations degrades severely at mixed Dirichlet-Neumann boundary junctions, characterized classically by an $\mathcal{O}(r^{1/2})$ leading singular function. While this linear behavior is well documented, the introduction of gradient-dependent semilinear perturbations alters the local asymptotic profile. This article proves that gradient penalties scaling linearly with $|\nabla u|$ induce a highly localized $\mathcal{O}(r^{-1/2})$ source term that resonates with the half-integer spectrum of the principal homogeneous differential operator. This non-orthogonal resonance causes standard separable polynomial assumptions to fail at order $\mathcal{O}(r^{3/2})$. We establish a generalized resonance theorem that forces the emergence of a logarithmic anomaly, providing the exact analytical formulation of the resulting $r^{3/2} \ln(r)$ profile alongside rigorous local Sobolev regularity bounds for the remainder. By calculating the exact logarithmic coefficients and angular offsets for $\ell_1$, $\ell_2$, $\ell_\infty$, and arbitrary $\ell_q$-norm penalties, we demonstrate the universality of this obstruction. Finally, we formalize the corresponding enriched continuous Galerkin space (XFEM), establish its quasi-optimality, and present finite element experiments that confirm the predicted localized pollution, recover the predicted logarithmic coefficients across the $\ell_1$, $\ell_2$, and $\ell_\infty$ penalties on a curvature-free flat junction, and show that resolving the junction recovers optimal degree-of-freedom efficiency in standard finite element solvers; we close by outlining the targeted software architectures required for a fully enriched implementation.