Papers

We address the numerical approximation of a scalar conservation law with a highly oscillatory flux, in the linear setting. Taking advantage of the intrinsic numerical viscosity of the scheme employed that yields a parabolic regularization of the original hyperbolic equation, we use homogenization theory to derive an effective equation which in turn we solve numerically. Our numerical experiments show the performance of the approach on the example of the Lax-Friedrichs scheme.

This paper introduces a hybrid joint-selective optimization (HJSO) framework for large-scale numerical problems in which a small subset of trainable quantities is of primary interest. We partition the full parameter vector into a high-dimensional remaining block and a low-dimensional block of parameters of interest (POIs), perform joint first-order optimization over the full parameter set, and then freeze the remaining variables while applying a reduced-space Levenberg-Marquardt (LM) refinement to the POIs. The method is designed for settings in which the POIs are low-dimensional but strongly influence the quality of the computed solution, while the full parameter space remains too large for full-space second-order methods. The framework is evaluated on three representative problems: a matrix eigenvalue problem, an inverse Bratu problem solved with a physics-informed neural network, and a 100-dimensional nonlinear Black-Scholes problem solved with the DeepBSDE method. In each test, HJSO reaches prescribed POI-error thresholds faster than the corresponding joint first-order baseline and improves the final POI accuracy for the reported solver configurations. The contribution is therefore not a universal optimizer, but a practical reduced-space strategy for problems with known low-dimensional parameters of interest and expensive high-dimensional training variables.

The conjugate gradient (CG) method is a classical iterative solver for sparse symmetric positive definite (SPD) linear systems, but its convergence strongly depends on the spectral properties of the system matrix. Preconditioning can improve these properties; however, designing preconditioners that generalize across diverse and irregular sparsity patterns remains challenging. We propose BandPC (Band Preconditioner Combinations, where "Band" denotes residual-norm bands), a data-driven framework that leverages graph neural networks (GNNs) to predict multi-stage preconditioning strategies for the flexible conjugate gradient (FCG) method. BandPC partitions the FCG iteration into three residual-norm-based bands and defines a structured search space of 125 candidate combinations over five classical preconditioners. By exploiting FCG's ability to switch preconditioners across iterations, BandPC learns to map matrix structure directly to a promising preconditioner sequence. To improve training and generalization, we introduce a soft-labeling mechanism that retains near-optimal combinations and normalizes their scores into a label distribution, together with a hierarchical feature representation that captures node-level attributes, edge-level algebraic coupling strengths, and global matrix statistics. Experiments on a hybrid benchmark of synthetic SPD problems and diverse matrices from the SuiteSparse Collection show that BandPC achieves the best iteration count and solution time on 41.3% and 33.1% of test matrices, respectively, even when compared with each matrix's individually best traditional preconditioner. These results demonstrate that learned residual-band preconditioner scheduling can effectively accelerate flexible CG solvers.

Reduced and mixed precision have moved from a niche optimization to a central design axis in scientific computing and engineering, driven by energy constraints, heterogeneous accelerators, and the convergence of simulation and machine learning. This paper organizes the landscape around seven coupled themes—number formats, floating-point emulation, emerging architectures, hardware/software co-design, relation to other approximations, software design, and precision as a multilevel resource —and, for each theme, synthesizes the state of the art, future directions, and open questions. We emphasize energy per trusted solution as the core objective, and we frame \say{recklessly responsible} computing as a pragmatic doctrine: exploit low precision aggressively, but with systematic detection, escalation, and certification pathways.

This work is concerned with the numerical computation of Gross-Pitaevskii ground states using the tensor train (TT) format. We regard the problem as the minimization of the Gross-Pitaevskii energy functional under a unit mass constraint, and propose a first-order Riemannian optimization scheme that preserves both the mass constraint and the low-rank representation throughout the optimization, by operating on the manifold of unit-norm functions of fixed TT-rank. On this manifold, we derive the energy-adaptive Riemannian gradient and show how to compute it efficiently in the TT format. Spatial discretization is performed using a spectral method with numerical integration, which allows the required computations to be carried out efficiently while preserving the low-rank format. The method is shown to substantially reduce computational time compared to full-rank computations, while maintaining accuracy for both single- and multicomponent Gross-Pitaevskii equations.

For a controllable single-input linear system in dimension $n$, the controllability Gramian eigenvalues, ordered decreasingly, are known to have the successive small-time orders $T, T^3, \ldots, T^{2n-1}$. We refine this leading-order hierarchy by explicitly computing the next-order coefficient of every eigenvalue and the first-order variation of its eigendirection and of the nested spectral subspaces. The resulting formulas have intrinsic expressions in the orthonormal Krylov basis: eigenvalue corrections describe the action of the dynamics along each Krylov direction, while variations of the spectral subspaces describe coupling to the next direction. We also establish local joint real-analytic dependence on $(T, A, b)$ near $T=0$ for the eigenvalues divided by their leading powers of $T$, consistently oriented eigenvectors, and spectral projectors. This provides convergent expansions with remainders uniform on compact families of controllable pairs. These results yield refined asymptotics for the worst-case minimum energy required to reach a unit target from the origin, the Gramian determinant and condition number, and the Ornstein–Uhlenbeck Gaussian profile in moving principal coordinates. The energy equals the reciprocal of the smallest Gramian eigenvalue, whose eigendirection identifies the most energy-demanding targets for sufficiently small times. We further derive an exact Gramian-weighted energy identity for the forward Fokker–Planck equation and sharp anisotropic short-time smoothing asymptotics using, respectively, the Lyapunov equation and an exact Fourier norm formula combined with the graded Gramian factorization. For symmetric dynamics, the spectral data reconstruct $A$ and determine $b$ up to its global sign.

We investigate numerically fluid flow through a thin porous layer, which is modelled as a lower-dimensional interface separating two bulk domains. At this interface, we impose effective transmission conditions involving homogenized coefficients rigorously derived via two-scale methods in previous works. These coefficients are determined by cell problems encoding the microscale properties of the porous layer. In a first step, we analyze qualitatively how the effective coefficients depend on the pore geometry. Based on these computations, we then present numerical simulations of the macroscopic models incorporating the effective interface conditions. Although replacing the thin porous layer by an effective interface substantially reduces the computational cost, the resulting nonstandard transmission conditions pose additional challenges for the numerical implementation. We employ a Taylor–Hood-type mixed finite element discretization with discontinuous pressure and tangential velocity across the interface and prove stability as well as error estimates for the resulting scheme. Finally, for microscale geometries that are still computationally resolvable, we compare the results of the macroscopic models with those obtained from direct simulations of the microscopic models.

We study the backward-Euler scheme for a class of stochastic functional differential equations (SFDEs) with memory driven by a fractional Brownian motion with Hurst parameter H > 1/2. We first establish the local consistency of the method, with a local truncation error of order H -$\rho$- $\beta$ + 1, and then, as the main result, prove the uniform global pathwise convergence of the numerical approximation, on a set of probability one, with order H -$\rho$ -$\beta$, for $\beta$ $\in$ (1-H, 1/2) and 0 < $\rho$ < H -$\beta$. A uniform high-probability formulation with deterministic constants follows as a corollary. In particular, the convergence order can be taken arbitrarily close to 2H -1, matching the order obtained for the explicit Euler scheme. The convergence proof relies on an exact integral representation of the scheme at mesh points, an a priori bound for the numerical solution, fractional calculus estimates, and a Gronwall inequality for weakly singular kernels.

State-dependent Riccati equation (SDRE) control requires repeated online solution of continuous-time algebraic Riccati equations (CAREs). A projected-channel $H_\infty$ formulation preserves the physical actuator and sensor dimensions of under-actuated, partially observed systems. We derive a computable sufficient attenuation bound for positive-semidefinite stabilizing Riccati solutions and construct an update that enforces their spectral-radius coupling condition. A discrete Kalman recursion supplies the scheduling state estimates. The two CAREs arising in controller synthesis are solved by a structure-preserving doubling algorithm (SDA) or a warm-started Newton–Kleinman iteration equipped with a doubling Lyapunov solver. Numerical studies use a twelve-state F-16 model with complete aerodynamic forces, moments, and trim, together with quadrotor spiral tracking. With instantaneous command updates, all three CARE backends give essentially the same closed-loop response. In the tracking simulations, measured computation times determine when new control inputs are applied, and the previous inputs are held during computation. SDA and Newton reduce airspeed and pitch errors in all six paired aircraft tracking runs. In the faster quadrotor task with limited computational resources, SDA and Newton complete all paired runs, while runs using MATLAB icare terminate early.

Stokes-Brinkman-Darcy (SBD) models connect free flow and porous flow through a Brinkman layer of finite thickness. The proposed method for solving these coupled equations consists of a mixed first-order formulation, shared interface traces and hard traction constraints, and a pressure correction based on boundary mass balance. Independent stresses, Darcy flux and pressure-gradient auxiliaries avoid second derivatives of network outputs. Boundary liftings and shared traces impose exterior and interface kinematic conditions exactly, while a stress map enforces native Brinkman-Darcy force balance. An integral momentum relation determines the pressure correction after training, without reference-pressure labels. Across two manufactured solutions and three permeabilities, mean corrected upper-pressure relative $L^2$ errors range from 0.87% to 2.58%; correction reduces both upper-pressure $L^2$ errors and the full Brinkman-pressure $H^1$ error in every run. Paired controls assess the contributions of interface sharing and traction enforcement. A non-manufactured filtration problem tests mixed-boundary adaptation against an independently refined finite-element reference, with all four runs satisfying the prescribed field and physical-consistency criteria.

Phase-field simulations of solid-state dewetting based on polynomial potentials may exhibit spurious bulk-diffusion coarsening, inconsistent with the surface-diffusion-dominated physics. To address this issue, we formulate a degenerate Cahn–Hilliard model with the logarithmic Flory–Huggins potential and dynamic contact line boundary conditions for the film–substrate–vapor triple junction. The logarithmic barrier confines the phase variable to its physical range and suppresses spurious coarsening at the continuum level. We develop a fully discrete, structure-preserving scheme that combines a Lagrange multiplier approach with an upwind finite-volume discretization and rigorously guarantees pointwise boundedness, mass conservation, and energy dissipation without artificial cut-offs or projections. A dimensional-splitting strategy further reduces computational cost while preserving these properties in each one-dimensional sweep. We also derive an explicit estimate of the equilibrium radius contraction caused by spontaneous film shrinking, showing that the logarithmic potential produces weaker spurious shrinkage than the polynomial potential. Numerical experiments confirm the analysis and demonstrate that, for small temperature parameters, the proposed scheme suppresses spurious coarsening and pinch-off while accurately reproducing surface-diffusion-dominated dynamics and relaxation toward equilibrium island structures.

Out-of-distribution (OOD) generalization is a central challenge in scientific machine learning. We study regression problems in which the test distribution differs from the training distribution and ask: under what assumptions on the target function or operator is stable extrapolation possible, and how far beyond the training domain can one extrapolate? Existing theory controls the test error through additive penalties measuring the discrepancy between the training and test distributions. Such guarantees show robustness to small distribution shifts, but can be very pessimistic in comparison to OOD performance observed empirically. We identify classes of holomorphic functions and operators for which the OOD generalization error converges at algebraic rates even in the presence of large distribution shifts. This phenomenon stems from the increasing smoothness of higher-index coordinates, leading to what we term a `blessing of high dimensionality'. For learning with either polynomials, deep neural networks or deep neural operators, we derive explicit rates for arbitrary test measures supported on suitable domains and quantify how the admissible domain depends on the underlying regularity of the function or operator. Our extrapolation guarantees are independent of the test distribution, depending only on its support. We also present a series of numerical experiments across a range of functions and operators that support the main theoretical findings.

Physics-informed neural networks (PINNs) solve partial differential equations (PDEs) by incorporating governing physical laws into the training loss. For evolution equations, however, their conventional pointwise space–time representation does not explicitly encode temporal dependence, which can hinder accurate prediction. To mitigate this limitation, this work proposes a novel neural architecture termed a causal integral neural network (CinNet). The core module of CinNet is a Volterra-type causal integral term, which aggregates historical features to encode temporal dependence, thereby incorporating temporal causality at the architectural level rather than through training-level modifications as in many existing methods. Building on CinNet, we further develop a causal integral physics-informed neural network (CI-PINN) for solving evolution equations. Extensive numerical experiments on benchmark evolution equations demonstrate that the presented method outperforms various baseline PINN variants in terms of solution accuracy, with pronounced superiority under sparse-collocation scenarios. Additional empirical analyses show that CI-PINN exhibits low sensitivity to hyperparameter choices, while ablation studies confirm the effectiveness of the proposed network components.

The storage and recovery of a vector of length \(2^d\) can become prohibitively expensive as \(d\) grows, a manifestation of the curse of dimensionality. For certain structured functions and multiscale quantities, the resulting discretized vectors admit compact quantized tensor train (QTT) representations after binary tensorization: a length-\(2^d\) vector is reshaped into a \(d\)-way binary tensor with low tensor-train (TT) rank. Motivated by this construction, we study compressed sensing of tensorized signals that combine bounded TT rank with sparse low-order interaction structure, using randomly sampled Walsh–Hadamard measurements. For this structured model class, we establish a uniform restricted isometry property and derive a measurement condition guaranteeing uniqueness of noiseless recovery. We also propose a structure-aware initializer for the alternating linear scheme (ALS), obtained by projecting the adjoint backprojection onto the low-order interaction space before applying TT-SVD. We prove a uniform initialization error bound whose measurement requirement, for fixed structural parameters, grows polynomially with the tensor order \(d\) rather than with the ambient dimension \(2^d\). Numerical experiments support the predicted RIP and initialization scaling and demonstrate more reliable ALS recovery.

We develop a structure-preserving parametric finite element framework for geometric flows whose energy density couples the unit normal and the curvature. This class includes bending energies with orientation-dependent rigidity or spontaneous curvature. A common fully discrete formulation treats closed curves in two dimensions and closed surfaces in three dimensions, and accommodates the $L^2$ flow, curve or surface diffusion, and the area- or volume-constrained $L^2$ flow. The formulation couples the geometric and curvature updates so that their contributions satisfy a discrete energy inequality. It uses continuous piecewise linear elements, a surface energy matrix, and a mass-lumped projection of the curvature derivative of the density. For positive densities satisfying a directional condition and convexity in curvature, we prove energy dissipation without a time-step restriction. The diffusion and constrained flows also preserve the enclosed area or volume exactly. The analysis allows non-even anisotropies and nonseparable dependence on orientation and curvature. Numerical experiments exhibit approximately second-order convergence in the manifold distance and confirm the discrete structural properties. Shape relaxation under an anisotropic Helfrich-type energy illustrates the use of the framework for coupled directional and bending effects.

This study investigates generalized boundary representations and transfer learning for physics-informed neural-network (PINN) estimation of spatially varying convective heat transfer coefficients (HTCs) in jet impingement cooling. Piecewise-linear (PWL) and neural-network boundary estimators are inferred from transient interior-temperature measurements and evaluated against CFD-CHT reference HTC profiles at Reynolds numbers of 5000, 10000 and 15000. A discrete Adjoint-FEM formulation with gradient regularization based on the $H^1$-seminorm provides an independent inverse-method comparison. Transfer learning is investigated under two source-target conditions: transfer between noise levels at the same Reynolds number, and transfer from a controlled Gaussian-profile problem to CFD-CHT-derived target problems at different Reynolds numbers. For the representative Re = 5000 target with 5% measurement noise, quadratic, PWL, and neural-network PINNs achieve relative HTC-profile errors of 5.08, 5.31 and 5.19%, respectively, while generalized representations reduce the error for the Re = 15000 target from 16.76% for the quadratic PINN to 14.41% and 14.46%. Transfer learning reduces target-training time by 57.86-81.71% for fixed-Reynolds-number transfer and by 50.59-84.94% for Gaussian-to-CHT transfer. Its effect on reconstruction accuracy is target dependent: fixed-profile PWL transfer remains within 0.30 percentage points of scratch training, whereas Gaussian-to-CHT transfer produces differences of 0.10-1.52 percentage points across targets and representations. These results demonstrate that generalized PINN representations can accommodate spatially varying HTC profiles without a prescribed quadratic form, while transfer learning can substantially reduce target-training cost for repeated inverse estimation, with the resulting reconstruction accuracy depending on source-target compatibility.

We develop a fully relativistic extension of the energy-conserving, unstaggered electromagnetic-potential particle-in-cell method introduced in our earlier nonrelativistic work. As in the nonrelativistic generalized-momentum formulation, the fields are advanced by a Crank–Nicolson discretization of the Lorenz-gauge potential equations, the current is deposited from particle orbits, and charge is advanced through the discrete continuity equation. The new ingredient is a generalized-momentum adaptation of the Higuera–Cary relativistic particle pusher. The vector-potential force is decomposed into an orbit-discrete-gradient part, which enforces the same finite-difference chain rule used in our earlier nonrelativistic formulation, and a skew part, which generates a Higuera–Cary magnetic rotation in mechanical momentum. For prescribed fields and frozen orbit data, this generalized-momentum Higuera–Cary map is a unit-Jacobian conjugacy of the mechanical Higuera–Cary map and preserves particle phase volume in $(\xx,\PP)$. For the fully coupled implicit step, the same orbit data define the current deposit, mesh-to-particle gather, vector-potential chain rule, and relativistic kinetic-energy secant velocity. The resulting particle work cancels the Crank–Nicolson field-energy balance, yielding exact relativistic total-energy conservation up to nonlinear solver tolerance, orbit-quadrature error, and roundoff. A self-adjoint Fourier projection removes modes on even-grid Nyquist planes while preserving the particle–grid work identity. Numerical tests are organized around cold relativistic two-stream and Weibel/filamentation instabilities, with diagnostics for energy drift, Gauss's law, the Lorenz gauge, particle–grid work, and the orbit chain rule.

This article is aimed at introducing a new procedure to define reference solutions to which the numerical approximations of anisotropic problems arising in plasma physics may be compared. In this context, the main issue is related to the construction of functions with no gradients along the vector field defining the anisotropy direction. The so-called solution manufacturing, implemented so far in the literature, only applies to frame- work where all the computations can be carried out analytically. In the present paper, we propose a new method in which this requirement is not mandatory. Indeed, via a Riemannian metric, we use a propagation of the boundary value by the flows related to the vector field defining the anisotropy direction. This permits to extend the application of the solution manufacturing to a broader range of contexts.

Neural operators are typically trained in a supervised fashion, which requires a dataset to be generated with a classical solver. Training them physics-informed, i.e., purely from the governing equations, removes this large offline cost and allows fresh samples to be drawn at every optimization step, but has so far been limited to simplified problems and trails supervised training in accuracy. The obstacle is the ill-conditioning of physics-informed losses, which differential operators induce and which worsens as the discretization is refined. We therefore propose a preconditioned residual loss function and show mesh-independent conditioning for elliptic problems and greatly improved conditioning for saddle point problems. Realized through geometric and algebraic multigrid, the construction applies to linear and nonlinear equations, steady or time-dependent, on structured and unstructured meshes, is agnostic to the neural operator architecture, and adds no cost at inference. On the Poisson, Allen-Cahn and stationary Stokes equations, the resulting label-free training matches supervised training and is four to twenty-five times more accurate than previous physics-informed operator learning methods.

In this work, we develop a second-order optimization framework for physics-informed neural networks (PINNs) applied to high-dimensional parametric partial differential equations (PDEs). The framework is built on the Gauss–Newton pullback metric, which provides an operator-informed notion of curvature in parameter space and connects the method to the broader family of natural gradient schemes. We show that, for coordinate-separable neural architectures and linear differential operators (or linearized operators in the nonlinear case) admitting a finite separable representation, the residual Jacobian inherits a structured separable factorization. This yields an exact compressed formulation of the Gauss–Newton step in a reduced space, without assembling the full residual Jacobian on the exponentially large tensor-product collocation grid. The dimension of the reduced space (the effective compressed dimension) is determined by the local collocation grid sizes, the separable operator structure, and the contraction pattern of the architecture, thereby avoiding dependence on the full tensor-product grid size and replacing dense linear algebra in parameter space by a substantially smaller structured problem. Within our framework, we investigate canonical polyadic and tensor-train parametrizations and derive their full algebraic characterization relevant to the Gauss–Newton method, including the structure of the residual Jacobian, the resulting compressed system, and its effective compressed dimension. Numerical experiments on high-dimensional PDEs, including parametric problems, demonstrate the high efficiency of the proposed compressed Gauss–Newton method, which achieves substantially lower errors than tensor-compressed first-order baselines with orders of magnitude fewer iterations and only a fraction of the computing time.

We propose a localized reduced basis method for parametrized parabolic partial differential equations based on an adaptive partitioning of the parameter domain using hierarchical binary trees. A POD-Greedy procedure is employed to construct reduced spaces associated with each parameter subdomain. A certified a posteriori error estimator is derived from an inf-sup stability analysis and shown to be equivalent to the underlying approximation error. This equivalence, together with suitable regularity assumptions on the parameter-to-solution map, enables a rigorous convergence analysis and is used to derive a strategy for selecting the number of POD modes to retain in each local reduced space. Numerical experiments for a parametrized convection–diffusion problem demonstrate the effectiveness of the proposed approach and confirm the predicted convergence behavior.

We present a high-order finite element arbitrary Lagrangian-Eulerian (ALE) method for hydrodynamics with slip-wall boundary conditions on analytically curved domains. The method combines distinct Lagrangian, remesh, and remap phases while preserving the geometric features of curved boundaries and allowing tangential motion of boundary nodes. In the Lagrangian phase, we use weak boundary condition enforcement to impose the slip-wall condition without freezing tangential motion. In the mesh optimization phase, we adapt TMOP (target-matrix optimization paradigm) to analytic walls by replacing selected boundary coordinate degrees of freedom with curve or surface parameters. The Newton solve then moves boundary nodes tangentially with respect to the prescribed geometry while carrying the drift introduced by the Lagrangian phase. In the remap phase, we transfer the solution from the old mesh to the optimized mesh by a conservative advection-based remap over pseudo-time. We present numerical results in two and three dimensions demonstrating the method on high-order curved meshes and showing that it preserves boundary fidelity while producing the expected hydrodynamic solution features.

We introduce Neural Harmonic Measure Operator (NHMO), a neural solver for elliptic PDE problems on variable-shape domains. The harmonic measure of a domain is the boundary probability distribution that, integrated against any boundary data, returns the Dirichlet Laplace solution. It depends only on the geometry, not on the boundary data. NHMO parameterizes the density of this measure as a transformer-based boundary kernel supervised by Walk-on-Spheres exit samples, so one trained kernel handles different boundary values on a shape with no retraining. We extend it to Poisson via a classical decomposition, with an auxiliary network amortizing the source-induced correction and avoiding the singular volume quadrature that breaks direct evaluation. At inference, new boundary values and new sources both yield PDE solutions by re-integration against the fitted kernel and lift, with no retraining. NHMO improves over four prior baselines on the MCB-B 3D variable-shape Poisson benchmark across all five categories, and is competitive with major neural-operator baselines on a controlled 2D testbed.

Least-squares residual minimization approximates a differential equation through an optimization problem. This paper investigates the conditions required to transition from stationary points and minimizing sequences of the residual functional to the exact solution. Starting from the classical orthogonality relation between the residual and the range of the linearized operator, we study how boundary conditions and parameterized variations determine what can be inferred about the residual at a stationary point. Smooth spurious critical points can occur even for well-posed problems; for linear equations, a compatibility condition yields an exact projection characterization. For a scalar conservation law with transonic boundary data, every minimizing sequence converges to a continuous function that does not solve the equation, even when a stationary entropy shock exists. We give conditions under which approaching this limit forces the parameter norms to diverge. A Hamilton–Jacobi example shows that a function can have zero residual without being the viscosity solution. Favorable results for a viscous model and uniformly convex Monge–Ampère critical points identify conditions that restore the corresponding implications. For piecewise-smooth trial families, the regularity needed to define the residual functional does not by itself justify the sampled gradient used in training. Moving residual jumps can contribute terms absent from that gradient; an explicit example has identically vanishing sampled gradients and a nonzero continuous derivative. Together, these results clarify what stationarity and sampled gradients imply about the residual, and under which conditions minimizing sequences lead to the intended solution.

We study numerical time integration for the motion of a charged particle in a strong magnetic field, focusing on the regime in which the fast gyration is not resolved by the time step. Standard Boris-type methods are structure preserving and accurate when the gyration is resolved, but their error constants deteriorate with increasing field strength; filtered Boris methods improve this behavior but may suffer from singularities and resonance-induced error blow-up. In this work we consider the case of a constant strong magnetic field and derive a broad family of one-step exponential integrators from the variation-of-constants formula, formulated in terms of matrix-valued filter functions. This framework includes existing filtered Boris variants and permits the construction of schemes with uniformly bounded filters. Using a discrete variation-of-constants representation together with summation-by-parts arguments, we establish error bounds whose constants are independent of the magnetic field strength and characterize the filter conditions required for first- and second-order accuracy. In particular, we show that bounded filters avoid resonance blow-up and allow arbitrary time steps, at the cost of a controlled order reduction near resonant frequencies. Numerical experiments confirm the predicted convergence behavior, the different accuracy of position and velocity components, and the role of the leading error terms across a wide range of step-size-field-strength products.

While structure-preserving methods for the forward simulation of Hamiltonian systems are well established, the mathematical foundations of their data assimilation and filtering counterparts remain comparatively less developed, particularly regarding rigorous accuracy guarantees. To address this gap, we propose an online filtering algorithm to reconstruct an unknown function $u^\dagger$ from finitely many measurements, assuming that $u^\dagger$ solves a parametric Hamiltonian PDE with unknown inputs. The method uses low-dimensional symplectic approximation spaces that evolve in time according to the PDE model and informed by the measurement data. The method can also be coupled with a dynamical sensor-placement strategy designed to optimize reconstruction stability. We derive error bounds in terms of stability and best approximation, and show that the reconstruction remains in a symplectic space and preserves energy up to the approximation error. Numerical experiments demonstrate the effectiveness of the method and the benefits of dynamic over static sensor placement.

We build upon recent developments in high-order spline-based geometric methods for the three-dimensional initial boundary value problem of Maxwell's equations. Previous schemes were limited to single-patch geometries. To address the computational challenges associated with inverting mass matrices in multi-patch settings, we extend existing techniques for 0-forms mass and stiffness matrices to 1-forms mass matrices. Our approach leverages the tensor-product structure of spline spaces to construct a spectrally equivalent preconditioner for the 1-forms mass matrices. Specifically, we design an efficient preconditioner for single-patch domains that is robust with respect to both mesh size and spline degree. In the multi-patch case, this methodology is extended by integrating the single-patch preconditioner with an additive Schwarz approach, ensuring robustness and scalability with respect to the number of patches. We analyze the computational complexity of the proposed preconditioner and validate its robustness and efficiency through comprehensive numerical experiments.

We develop a high-order fully discrete ALE-SAV finite element method for sharp-interface two-phase Navier–Stokes flow that satisfies a discrete energy dissipation law. We introduce an SAV reformulation to ensure the stable treatment of the surface and gravitational potential energies, together with an energy-stable ALE transfer operator that controls the density-weighted kinetic energy across moving meshes. The resulting scheme is linear and satisfies a modified energy-dissipation inequality. Numerical experiments, including rising-bubble benchmarks, demonstrate the robustness of the method and confirm second-order convergence in time and arbitrarily high-order convergence in space.

For probability measures on $\mathbb{R}^d$ supported in fixed compact sets, we prove that quadratic optimal transport plans are quantitatively stable in Wasserstein distance under perturbation of both marginal measures, assuming that one of the initial measures satisfies an upper Ahlfors regularity condition with exponent strictly greater than $d-1$. Under the same assumption, we also prove novel quantitative stability results for optimal transport maps. Furthermore, we prove a coercivity theorem which states that any probability measure on the product space must be quantitatively close to the set of optimal plans, if it has similar marginals and a similar transport cost to the optimal value. Finally, we prove a quantitative uniqueness theorem which acts as a quantitative counterpart to Brenier's theorem. The Wasserstein diameter of the set of optimal plans is controlled by the distance of one marginal measure to a regular measure for which uniqueness holds. In this way, "almost uniqueness" of optimal plans is quantified by the source measure being "almost regular". Examples are provided which prove that the exponents of all estimates are sharp.

In this paper we study the numerical approximation and convergence analysis of a steady-state heat conduction problem with mixed boundary conditions. The physical model is governed by a hemivariational inequality depending on a heat transfer parameter $\alpha > 0$. We consider the finite element approximation of this penalized problem, as well as its corresponding limit problem with a prescribed constant temperature on a part of the boundary. The main theoretical result establishes the strong convergence of the discrete solutions. Specifically, we prove the double-limit convergence of the finite element approximations to the limit solution as the mesh size $h$ tends to zero and the penalization parameter $\alpha$ tends to infinity, independently and simultaneously. We further complement the convergence analysis with an error estimate of optimal order for the discrete limit problem, and we explain why an estimate uniform in the penalization parameter cannot be expected. The theoretical results are illustrated by numerical simulations on four examples verified by the method of manufactured solutions.