Papers

In this work, we consider a class of continuous-time autonomous dynamical systems that model various important phenomena and processes encountered in real-world situations. We construct a second-order nonstandard finite difference (NSFD) method that simultaneously preserves two essential properties of the dynamical systems for all finite step sizes, namely the positivity of the solutions, the set of equilibrium points and their asymptotic stability. This NSFD method is constructed based on an appropriate choice of nonstandard denominator functions and a weighted discretization of the right-hand side functions. Under easily-verified conditions, the denominator functions guarantee second-order convergence, whereas the weights ensure the dynamic consistency. By taking advantage of the specific structure of the right-hand side functions, a simple discretization is utilized instead of the nonlocal discretization approaches commonly used in previous works. This simplifies the construction of the proposed NSFD method and, in particular, makes its asymptotic stability analysis easier. As an illustration and an important application, we apply the constructed second-order NSFD method to a well-known two-stage structured species model with recruitment. Consequently, a simple second-order NSFD scheme for the considered two-stage structured species model is derived, improving upon a first-order NSFD scheme constructed in a previous work. Numerical experiments demonstrate the advantages of the second-order NSFD scheme over a standard second-order numerical method, namely, the explicit trapezoidal method. The proposed NSFD method is simple and can be applied to a broad class of dynamical system models arising in both theory and applications. Moreover, it can be readily combined with the Richardson extrapolation technique to improve its accuracy.

An efficient aggregation-based multigrid approach is developed in this work to solve linear algebraic systems generated by discretizing the Landau-Lifshitz equation in micromagnetics. The discretization couples the first-order backward differentiation formula with first-order extrapolation, and the underlying equation describes magnetization dynamics within micromagnetic simulations. For these linear systems, standard iterative methods and conventional multigrid solvers intended for symmetric problems suffer from deteriorated convergence under mesh refinement and poor computational performance. Existing algorithms struggle to handle the system's intrinsic features, namely its sparse spectral properties and non-symmetric matrix structure. To overcome these limitations, a multigrid framework is constructed, where smoothing and coarse-grid correction operators are customized for the target linear system. Numerical tests confirm the robustness and mesh-independent convergence of the resulting method. Compared with well-established Krylov-subspace solvers and conventional multigrid techniques, the aggregation-based multigrid solver cuts down iteration counts and overall computational cost considerably, while producing faithful representations of magnetization dynamics.

We develop a rigorous mathematical framework for the governance of systems of K self-adapting generative AI models under the principles of Model Risk Management (MRM). When multiple models share a meta-learning coupling through an interaction matrix, the per-agent Lyapunov analysis that underpins standard MRM is provably insufficient: individual agents can each satisfy their declared stability bounds while the joint system is in a regime of emergent ensemble-level drift. We formalize this gap through the Joint Lyapunov Proof (JLP)—a cryptographic and stochastic protocol that attests, without revealing proprietary weights, that the aggregate dynamics satisfy MRM Ongoing Monitoring standard at every validation epoch. Our main contributions are the following. We give a complete characterization of the infinitesimal generator of the joint quadratic Lyapunov function. We derive the exact critical coupling threshold above which the system loses mean-square stability. We prove a Noise-Floor Theorem and identify the correct target for zero-knowledge attestation. A per-epoch Succinct Non-Interactive Argument of Knowledge (SNARK) on the live weights is derived. All theoretical claims are validated against five numerical studies using a multi-agent softmax system.

Forward uncertainty propagation in complex physical systems can induce structured covariance across field-valued outputs. For a probabilistic surrogate, the total predictive covariance comprises the covariance of conditional means across input realizations and the average conditional predictive covariance. Probabilistic DeepONet (Prob-DeepONet) provides lightweight uncertainty quantification by predicting pointwise Gaussian means and variances in a single forward pass, but its conditional predictive covariance is restricted to a diagonal form. To represent cross-location conditional dependence without explicitly parameterizing a full high-dimensional covariance matrix, we develop a two-step mean-variance DeepONet (two-step MV-DeepONet) through two principal modifications. First, two-step training is used to decouple output-basis learning from the input-to-coefficient mapping, together with basis orthogonalization and subspace rotation. Second, Gaussian probabilistic modeling is transferred from the high-dimensional physical output space to the low-dimensional rotated coefficient space. Mapping these probabilistic coefficients through the shared basis induces a generally non-diagonal conditional predictive covariance in the physical output space while retaining single-pass inference. A Frobenius-norm error decomposition and corresponding upper bound identify low-rank covariance compressibility, trunk-subspace approximation, finite-sample statistical error, and coefficient-space covariance estimation as the principal factors governing covariance recovery. Numerical experiments on three representative problems governed by partial differential equations (PDEs) and a hypersonic blunt-body aerothermal problem show improved generalization, more structured uncertainty bands, and accurate recovery of off-diagonal correlation patterns compared with Prob-DeepONet.

Model risk management (MRM) guidance assumes a static model lifecycle, in which models are developed, independently validated, and implemented without further autonomous modification. Continually self-adapting generative AI systems — models that update their own weights during production deployment — fundamentally violate this assumption and render point-in-time validation inadequate. This paper addresses the resulting governance problem in two parts. Part I develops a rigorous telemetry architecture for such models, operating simultaneously in discrete and continuous time. We establish a Minimal Sufficient Statistic for audit purposes, construct a tamper-evident Merkle chain for discrete weight sequences, derive the appropriate continuous-time generalization via the Ito formula, and propose event-driven logging via KL divergence stopping times that is both computationally tractable and meaningful for validation. Part II asks what happens when the model provider is adversarial. A firm deploying such a model has strong incentives to conceal learning updates that would trigger mandatory validation review. We formalize this as the Model Hiding Problem and provide a systematic taxonomy of six distinct attack strategies against the Part I architecture, spanning discrete and continuous time, with a formal countermeasure for each. Together the two parts establish a dual-regime architecture in which continuous telemetry is necessary but not sufficient, narrowing—but never replacing—periodic invasive audit. The framework is model-architecture-agnostic and is designed to satisfy the three pillars of traditional MRM.

We consider the probability that a discrete random matrix $M_n(\xi)$ is strongly non-singular, meaning all its leading principal submatrices are non-singular. This property is equivalent to the existence of an LU factorization. We show that for any discrete random variable $\xi$ with finite support and $|\xi|_\infty < 1$, there is a constant probability that $M_n(\xi)$ is strongly non-singular with a growth factor bounded by $n^{5/2+\delta}$. Furthermore, we provide a tight asymptotic lower bound for this probability as $|\xi|_\infty \to 0$. Finally, we provide exact counts for strongly non-singular binary matrices up to $n=9$ and use these to derive improved upper bounds for the Bernoulli case.

Each quasi-static load step of phase-field fracture is a bound-constrained minimization of a nonconvex, coupled displacement-damage energy under an irreversibility bound on the damage. Monolithic Newton stalls once the nonlinearity localizes at the advancing crack front, and staggered (alternate-minimization) schemes converge slowly there. We present an on-demand nonlinear-elimination preconditioned trust-region Newton method: an energy Steihaug-Toint trust region, a primal-dual active set for irreversibility, and a bound-constrained field-split sweep that eliminates an algebraically-identified "hard set" spanning both fields before each step. The elimination is applied on demand – triggered by the coupled Newton's own stalling and otherwise skipped – so the method reduces to monolithic Newton at no surcharge where the step is already healthy. We find the robustness to come from the energy trust region: with that globalization fixed, monolithic Newton already completes every loading history without cutbacks, where residual-merit Newton death-spirals, alternate minimization stalls, and the full-field sweep loses robustness. Against that well-globalized baseline, the on-demand elimination cuts outer nonlinear iterations by 19-25% (brittle) and 17% (ductile), with always-on elimination reaching 26-28% and about $39\%$ at the ductile nucleation step. Measured machine-independently, as a full-mesh-equivalent assembly-work proxy rather than wall-clock, it is competitive with – not faster than – monolithic Newton (within about 10%), whereas an always-on sweep adds up to 30%. Nonlinear elimination is thus an iteration-reduction mechanism whose overhead the on-demand gate bounds, with no demonstrated total-work advantage over well-globalized monolithic Newton.

An annular liquid jet may enclose a volume of gas and, when surface tension dominates over inertia, close on the symmetry axis at a finite distance. In the one-dimensional membrane model, closure provides an algebraic relation among the state variables at the free end. Earlier domain-adaptive schemes differentiate this relation to obtain an evolution equation for the domain length. We instead absorb the endpoint closure identically in the dependent variable and discretise the reduced system by Chebyshev collocation. Convergence is geometric, published steady data are recovered, and the unsteady calculations are checked by componentwise off-grid residuals, temporal refinement, and an independent upwind finite-difference discretisation. For jets of finite thickness, the convergence length decreases with nozzle thickness when the Froude number, Weber number, nozzle angle, and pressure coefficient are held fixed and the terminal zero remains transversal; this is a parameter-specific model result rather than a universal experimental thickness law. The forced response has a pressure-coefficient maximum near a Strouhal number of 0.093, whereas the convergence-length gain increases monotonically over the globally admissible frequency interval studied. Most importantly, a global positivity check changes the interpretation of the first finite-time event. Under body-force forcing, the mean radius first reaches zero at an interior point while the tracked endpoint remains transversal. The interior contact is quadratic, and the independent finite-difference scheme converges to the same first-contact time. After this event, the tracked endpoint is no longer the first physical closure point, so the single-interval annular-jet representation is no longer globally admissible.

Many computational problems possess the following three features: (1) the problem could be solved efficiently if a linear operator could be diagonalized, (2) direct diagonalization is infeasible due to scale, but, (3) the operator is known to commute with a permutation since the problem is (a) symmetric with respect to a change of coordinates, or (b) the operator is block circulant. For example, high-resolution template matching demands repeated convolution of large, high resolution images with an exhaustive list of related templates. The associated linear operator could be compressed, via a low rank approximation, if diagonalized. The required decomposition expensive, but, the entire problem is symmetric to in-plane rotations. In these cases, the eigenfunctions are constrained by the symmetry. These constraints allow efficient decomposition. We illustrate a parallelized algorithm that allows fast diagonalization of any block-circulant matrix. When the index space can be partitioned into $l$ classes of $m$ interchangeable elements, the algorithm reduces storage costs by a factor of $m$ and, if $w$ workers are available, the algorithm demands $\mathcal{O}(l^2m\log(m)/w) + \mathcal{O}(ml^3/w)$ floating point operations per worker yielding a $m^2$ speedup. We use this procedure to decompose a high precision template matching matrix. We compare runtime on a reduced problem, where the fast procedure ran 205 times faster per recovered feature, for 22.5 times as many features. On a typical template matching matrix, the decomposition required 25 times less time than computing the matrix entries. We also demonstrate decomposition of a $\sim$30 times larger matrix that covers the complete set of possible projections represented in a cryo-EM image at 2~Å resolution in 14 minutes. This procedure is more stable and allows recovery of in-plane rotations to arbitrary precision.

Recent work of Colbrook and Stepaniants produced the first counterexamples to the planar Pompeiu and Schiffer conjectures and introduced the conformal fixed-disc, disk-polynomial, and validated-tail machinery used here. By adapting this framework to the complementary Dirichlet endpoint, we disprove the unrestricted planar Berenstein conjecture. Specifically, we construct a bounded simply connected domain $\Omega$ with real-analytic Jordan boundary, which is not a disc and for which there exist $k\in(27.4381178838,27.4381198839)$ and a nonzero real-valued function $u\in C^\omega(\overline\Omega)$ satisfying $(\Delta+k^2)u=0$ in $\Omega$, with $u=0, \partial_\nu u=\text{constant}\ne0$ on $\partial\Omega$. Thus the overdetermined Dirichlet–Neumann data do not characterize the disc without an additional sign assumption on $u$. The domain has dihedral symmetry of order $26$, but is neither a disc nor centrally symmetric, and the corresponding eigenfunction changes sign. Equivalently, its boundary arclength measure satisfies $\widehat{\sigma_{\partial\Omega}}(k\omega)=0$ for $\omega\in\mathbb S^1$. After conformally transferring to the unit disc, exact support identities and quantitative disk-polynomial estimates yield rigorous control of the infinite-dimensional tail. A Newton–Kantorovich argument then reduces existence to finitely many explicit inequalities, which are certified using interval arithmetic. The extension from the Pompeiu–Schiffer problem is not formal. The earlier construction absorbs both boundary conditions into a single inverse-Laplacian equation. At the Dirichlet endpoint considered here, the nonzero Neumann datum forces the harmonic source modes to remain, producing a coupled interior–boundary system involving the full zero-Dirichlet inverse and its Neumann trace, together with a separate sign-recovery problem.

We consider the problem of PDE discovery from possibly noisy observations under the hypothesis that the underlying dynamic is symmetric in all rigid motions. Rather than using a generic library of derivative monomials, we leverage this assumption to construct libraries whose candidate terms are themselves rigid-motion-equivariant, and combine them with sparse regression to benchmark such libraries over five different equations. The advantage is most pronounced when the noise itself breaks rigid-motion symmetry (e.g., radially or axially varying noise), and when the ambient spatial dimension increases, in which case they are also less resource intensive.

Late-time wave tails decay at different rates along future null infinity and along timelike worldlines at finite radius. A compactified numerical evolution must represent both the slower decay at null infinity and the faster interior decay, producing an increasingly sharp transition between the two regimes. We address this difficulty for semilinear wave equations in Minkowski spacetime using homothetic hyperboloidal coordinates adapted to the scaling structure of the tail. In these coordinates, the tail approaches a smooth radial profile with the same decay rate at every compactified radius. The formulation therefore avoids the steepening of the radial profile seen in stationary hyperboloidal evolutions, and it reaches late times in a number of steps that grows only logarithmically with retarded time. We demonstrate this approach using pseudospectral simulations in 3+1 dimensions and reproduce the generic decay rates conjectured by Rinne. We also provide numerical evidence consistent with a nongeneric codimension-one cancellation of the leading tail coefficient at null infinity, resulting in a faster decay rate.

Volumetric parameterization, the process of mapping a 3-manifold onto a simplified volumetric domain, is important for many tasks in computer graphics and imaging science. However, most prior volumetric parameterization approaches have only utilized standardized domains such as a solid ball regardless of the overall shape of the given 3-manifolds, which introduces significant geometric distortion and affects the subsequent shape processing and analysis tasks. To overcome this issue, in this work we propose a novel volumetric parameterization framework for simply connected 3-manifolds. Specifically, the proposed framework jointly controls local shape and mass distortions, while adapting the target domain during the optimization process. It enables three progressively more flexible target-domain settings for the parameterization: a prescribed solid ellipsoid, a volume-normalized adaptive ellipsoid with variable radii, and a sea-embedded free-boundary domain. For each setting, the parameterization algorithm consists of a 3D quasi-conformality shape update, a diffusion-based density-equalizing update, and a geometric correction procedure for removing element foldings, thereby allowing for volumetric parameterizations with different desired effects. Experimental results are presented to demonstrate the effectiveness of our proposed framework. Moreover, our framework can be easily applied to multiresolution and localized adaptive volumetric remeshing, volumetric registration, and volumetric morphing. Altogether, our work provides a new way for the representation, processing, and analysis of 3-manifolds.

Fractional elliptic stochastic partial differential equations (SPDEs) are widely used in statistics and machine learning for computationally efficient and flexible modeling of Gaussian random fields. The computational efficiency of the SPDE approach relies on finite element approximations combined with numerical quadrature and mass lumping, which enable sparse matrix methods during inference. Although many works have studied finite element approximations of fractional SPDEs, the effect of the mass lumping and quadrature approximations used in practice has not been fully analyzed. To fill this gap, we derive convergence rates for numerical approximations of fractional SPDEs based on finite element discretizations combined with numerical quadrature and mass lumping. Specifically, we obtain explicit convergence rates for the mean-squared error of the covariance function in a general framework that covers the main settings where mass lumping is used in the SPDE approach. We also analyze non-stationary variance-control factors of the form $L^\beta(\tau u)=\mathcal{W}$, where $\tau$ is spatially varying, and derive covariance error estimates showing how the regularity of $\tau$ affects the convergence rate. As specific examples, we provide results for random fields on bounded Euclidean domains, Riemannian manifolds, and metric graphs. Numerical experiments are presented that confirm the theoretical results.

Coarse finite-volume averages do not generally form a predictive state: realized interface fluxes close a conservative update, but distinct fine-grid states with identical parent averages can generate different future coarse histories. We analyze this failure for periodic scalar advection discretized by a first-order upwind finite-volume method with forward Euler time integration. We derive a finite-horizon observation-rank law: each additional observation exposes one new child-cell layer and contributes one fewer independent direction than the number of parent cells, until the unresolved layers are exhausted. Centralized prediction therefore requires one fewer additional coordinate than the number of parent cells per exposed layer, whereas product-local prediction requires one coordinate per layer in each parent. An anchored flux-divergence queue attains the centralized bound, identifies the periodic flux gauge, and, at saturation, forms a minimal autonomous predictive state with the parent averages. We then distinguish exact observability from stable recoverability. The collar-to-queue map becomes rapidly ill-conditioned as the Courant number decreases, so algebraically visible delayed information may fall below a prescribed numerical tolerance. For Courant numbers strictly between zero and one, we prove contraction of the nonconstant component under bounded arithmetic perturbations, with separate control of mean drift, an explicit perturbation neighborhood, and a grid-dependent decay time. Numerical experiments illustrate the rank ladder, effective-rank loss, queue conditioning, step-function flattening, and delayed coarse separation followed by dissipative decay. The results provide a solvable benchmark for assessing state sufficiency in coarse, reduced, multiscale, and learned scientific models.

Fourier neural operators (FNOs) provide efficient nonlocal spectral learning, but varying geometries and independently chosen discretizations remain difficult to accommodate. We introduce the ambient-domain extension Fourier neural operator (ADEx-FNO), a deterministic framework that incorporates geometry without modifying the defining Fourier-operator layers. Each physical domain is embedded in a fixed ambient hypercube and represented by a signed distance function. Inputs and solution fields are deterministically extended to the ambient domain, transferred to a common, potentially nonuniform rectilinear latent grid, processed by the FNO, then interpolated to an independently chosen target discretization and restricted to the physical domain. All geometry-transfer operations lie outside the optimization procedure and require no trainable graph, point-cloud, deformation, or geometry-decoding modules. ADEx-FNO achieves relative l2 errors of 0.32%-0.77% on held-out smooth-domain nonlinear Poisson and advection-reaction-diffusion problems in 2D and 3D, and is also evaluated on unseen nonsmooth geometries. A single ADEx-FNO inference is then used to initialize conventional CFD solvers. For all 29 converged 2D and 3D RANS cases, pseudo-time iterations decrease, with mean reductions of 44.17% and 43.03%, respectively, with comparable gains across three mesh resolutions. URANS cases reduce post-window physical-time advances by 18.52%-27.51%. In transfer from 2D URANS training data to DNS at different Mach and Reynolds numbers, the bootstrap interval decreases by 23.47%-48.21%, depending on the target statistic. In all CFD tests, ADEx-FNO provides only the initial field; the governing-equation solver controls the subsequent solution, while physical or statistical consistency is assessed separately from computational savings.

Differentiable simulation is a key component in learning, control, and inverse problems, where gradients through nonlinear implicit solvers are required. Existing approaches either rely on unrolled automatic differentiation, whose memory grows with solver depth, or on equation-level implicit differentiation, which assembles global Jacobians and solves large sparse adjoint systems, discarding the locality of the forward solver – and differentiating the converged equation rather than the finite computation that actually ran. We propose solver-level differentiation, which differentiates the executed solver itself. When a solver is composed of block implicit updates, its discrete adjoint is obtained by applying the corresponding adjoint updates in reverse order, yielding a reverse-sweep formulation whose backward pass mirrors the forward solver. From an operator perspective, the forward pass realizes an approximate inverse through ordered local solves, and the backward applies its transpose through reverse local adjoint solves, constructing no global system. We instantiate this idea on Vertex Block Descent, yielding a differentiable solver whose reverse colored Gauss-Seidel sweeps are composed entirely of local $3\times 3$ adjoint solves. The backward matches automatic differentiation through the identical executed forward to machine precision at every solver depth, where the equation-level adjoint is off by 37% after one sweep; in a controlled same-codebase, same-GPU comparison it is 33x faster and uses 71x less memory than unrolled automatic differentiation; and the same construction is exact on projective dynamics and extended position-based dynamics. We scale differentiable elastodynamics to $10^6$ contact-coupled soft bodies (8M vertices) on one GPU. Overall, this work highlights solver structure as a practical organizing principle for efficient differentiable simulation.

The Generalized Optimized Schwarz Method (GOSM) originally proposed in [Claeys, 2021] is a variant of Deprés algorithm, a domain decomposition strategy for the solution of harmonic wave propagation problems. It imposes transmission conditions through interfaces by means of a non-local exchange operator. Conducting our analysis at the continuous level, in an infinite dimensional setting, we propose a coarse space construction for the preconditioning of the GOSM formulation, and provide estimates for the convergence of GMRes applied to the preconditioned equation.

We assemble numerical algorithms to search for $J$-holomorphic curves in symplectic manifolds. Each algorithm employs several different numerical techniques, each technique addressing a different aspect of the geometric problem. We separately consider both classical Fourier expansion and deep neural networks in our algorithms and compare their performance. Our algorithms take as input a smooth curve in a given homology class and search for a $J$-holomorphic curve in the same homology class. We first verify we can produce explicitly known holomorphic curves in complex manifolds, for example the Weierstrass $\wp$ function on the torus and curves in $S^2\times S^2$ with the standard complex structure. Then we search for $J$-holomorphic curves in $S^2\times S^2$ with non-integrable almost complex structures: essentially we start with a known holomorphic curve in an integrable almost complex structure $J_0$, deform $J_0$ to a nearby nonintegrable almost complex structure $J_\epsilon$, and use our methods to find the nearby $J_\epsilon$-holomorphic curve.

Linear multistep methods (LMMs) combined with neural-network approximation provide a high-order framework for learning governing vector fields of dynamical systems from discrete trajectory data. This paper studies two issues in LMM-based discovery that are not resolved by existing grid-level convergence theory. First, in non-auxiliary Adams–Bashforth (A-B) and Adams–Moulton (A-M) discovery systems, we observe that zero-residual grid solutions are nonunique but consistent along the trajectory, with differences limited to the boundary layer. We explain this phenomenon through a kernel analysis of the non-auxiliary discovery matrices. Under the corresponding discovery-stability conditions, the differences between zero-residual grid solutions are exponentially localized near the initial indices for A-B schemes, whereas they form two-sided boundary layers near the initial and terminal indices for A-M schemes. Second, we derive whole-trajectory generalization estimates for both auxiliary and non-auxiliary formulations. Once the learned vector field is restricted to a fixed observed trajectory, each component of the error becomes a scalar function of time. For the auxiliary formulation, the trajectory error is $O(h^p)$ under the corresponding grid accuracy, approximation, and trace regularity assumptions. For non-auxiliary formulations, the global estimates contain additional boundary-layer terms. On fixed interior subintervals, these terms are exponentially damped. Numerical experiments illustrate the convergence behavior.

We establish the global convergence of the (non-homogeneous) Fitness Centrality algorithm for general graphs, deriving an explicit convergence bound for the corresponding fixed-point iteration. Furthermore, we show how the convergence can be dramatically improved by Anderson acceleration and by switching to Newton's method once a sufficiently good approximation to the fixed point has been found. The efficacy of this strategy is illustrated by numerical experiments on different types of graphs.

Tumor progression is an inherently multiphysical phenomenon in which interstitial fluid dynamics, biochemical transport, and cellular mechanics interact across multiple spatiotemporal scales. Classical mesh-based solvers, although accurate, impose prohibitive computational costs for the repeated evaluations demanded by inverse parameter identification and future patient-specific predictive pipelines. In this work we introduce a Physics-Informed Neural Network (PINN) framework for a tractable chemo-fluidic continuum model of tumor growth that couples an advection-diffusion-reaction (ADR) equation for the tumor volume fraction with a quasi-static Darcy pressure equation for the interstitial fluid pressure. By intentionally decoupling the solid-mechanical equilibrium, we obtain a three-equation system whose gradient structure is stable under automatic differentiation, enabling robust deep-learning optimization. The network simultaneously learns both state variables from physics constraints alone (forward problem) and recovers hidden transport parameters from sparse, noisy synthetic measurements (Data-Assimilation PINN, DA-PINN, inverse problem). We verify the forward solver against a high-resolution finite-difference (FD) reference, achieving a mean absolute error below 0.002. For the inverse problem, starting from an initial permeability estimate of 0.08 (a factor of 4x above the true value of 0.02) with only 5% spatially sparse observations corrupted by 5% Gaussian noise, the DA-PINN recovers the permeability with a relative error below 5%. These results demonstrate that physics-informed deep learning constitutes a viable, computationally efficient route to multiphysics oncology modeling and lays the mathematical groundwork for future integration into clinical data assimilation pipelines.

This paper presents robust numerical formulations for suppressing low-frequency cancellation errors in the magnetic field integral equation (MFIE) for perfect electric conductors and the Müller integral equations for homogeneous and composite dielectric objects. The proposed approach combines quasi-Helmholtz projectors, which separate the Helmholtz components of the surface currents, with a tailored rescaling procedure that preserves the components essential for accurate far-field evaluation from finite-precision cancellation. The resulting formulations retain the conditioning properties of the mixed MFIE and Müller discretizations in the dense-mesh and low-frequency regimes. Numerical experiments involving simply and multiply connected scatterers, composite dielectric structures, multiscale and high-contrast configurations demonstrate that the proposed methods yield accurate far-field predictions across a broad frequency range, from the full-wave regime to extremely low frequencies.

We propose and analyze a frequency-stable boundary integral equation method for scalar wave scattering by periodic diffraction gratings in two dimensions. The standard quasi-periodic Green function is singular as a function of the frequency at Wood anomaly frequencies, which prevents a direct frequency-uniform error analysis of conventional integral equation discretizations. Our method is based on a modified quasi-periodic Green function obtained by subtracting the singular Rayleigh modes and introducing a finite-dimensional correction that restores the radiation condition. The resulting combined-field boundary integral equation is shown to be uniquely solvable whenever the underlying diffraction problem is uniquely solvable, and the corresponding scattered field depends continuously on the frequency, including at Wood anomalies. We then construct a fully discrete Nyström scheme using an Ewald representation of the Green function and a quadrature rule adapted to the logarithmic singularity of the kernel. For compact frequency intervals that may contain Wood anomaly frequencies, we prove unique solvability of the discrete system for sufficiently large discretization and truncation parameters and establish a frequency-uniform error estimate and bound on the condition number. Numerical examples illustrate the convergence of the scheme and its stability across Wood anomaly frequencies.

We present a rigorous error analysis for a semidiscrete immersed CutFEM applied to fluid–structure interaction (FSI) problems involving an incompressible viscous fluid and an elastic solid. The method employs a Lagrange multiplier to enforce the fluid-solid coupling conditions while fully accounting for the effect of the structural deformation on the fluid domain. For sufficiently regular FSI solutions and finite elements of degree $k\geq3$, we prove convergence of order $k$ for the fluid velocity and solid displacement and convergence of order $k-1$ for the fluid pressure in the energy norms associated with the FSI problem. To the best of our knowledge, this is the first rigorous convergence analysis of an unfitted finite element method for a fully coupled FSI system with a deforming interface.

We present a neural network-based method for numerical analytic continuation of solutions of nonlinear partial differential equations (PDEs) and detection of their complex singularities. The proposed framework employs "unsafe" Padé activation units (PAUs) as activation functions, together with a novel backpropagation-free method for computing the weights and biases of the hidden layers. The training procedure is designed specifically for learning meromorphic functions with pole-type singularities. Unlike existing methods with fixed hidden-layer parameters, the proposed approach treats the weights and biases as time-dependent functions, allowing the network to adapt to the evolution of the complex singularities. Using this method, we can efficiently locate the complex singularities of the extended solution, track their trajectories over time, and, based on this dynamics, infer the formation of characteristic phenomena in the solutions of nonlinear PDEs. To demonstrate the performance of our method, we analyze three well-known cases: (1) the nonlinear heat equation, which exhibits finite-time blow up; (2) the nonlinear Burgers equation, which develops a shock; and (3) the nonlinear Schrödinger equation, in which rogue waves form.

This paper presents a family of high-resolution weighted essentially non-oscillatory compact least-squares schemes with implicit time integration for the compressible Navier-Stokes equations on curvilinear grids. Compared with the original compact least-squares schemes, the proposed method introduces two main improvements. First, instead of enforcing the accuracy constraints over the entire computational domain, including discontinuous regions, it constructs the reconstruction matrix only along smooth reconstruction lines using an accuracy-preserving weighting strategy. This treatment effectively suppresses the persistent high-frequency oscillations observed in the original compact least-squares schemes and yields sharper profiles near discontinuities. Second, the method simplifies the shock-capturing procedure and improves efficiency by using a single set of polynomials, whereas the original compact least-squares schemes require both unlimited and limited polynomials. Combined with spectral optimization, the proposed method exhibits more favorable spectral properties than conventional weighted essentially non-oscillatory schemes. The smoothness indicators and penalty matrices are constructed through an efficient iterative procedure with modest additional cost. Numerical results for inviscid and viscous one-, two-, and three-dimensional flows demonstrate that the proposed method provides robust shock-capturing capability while maintaining high resolution in smooth regions and across contact discontinuities.

Drifting models are a recent class of one-step generative models that evolve the model distribution during training using a predefined sample-based drift field. Although they avoid iterative inference, their kernel-based drift fields induce frequency-dependent training dynamics: In the linearized regime, each Fourier mode of the density residual decays at a rate determined by the kernel spectrum, leading to slow recovery of fine-scale structure. We propose Second-Order Drifting Models, which lift drifting dynamics into phase space by augmenting generated samples with artificial velocity variables. We show that the resulting density perturbations obey accelerated second-order dynamics in Fourier space, connecting drifting models to the celebrated Nesterov acceleration from optimization theory. This provides a principled mechanism for mitigating the spectral stiffness of first-order drifting while preserving one-step inference. We derive a practical semi-implicit training algorithm and evaluate it on synthetic distribution matching, sequential data generation, and robotic control. Across these settings, the second-order drifting model improves convergence behavior and achieves competitive or superior performance over first-order drifting baselines.

We connect different degrees of smoothness of real-valued periodic functions to the number of qubits required for their high-precision Fourier-basis amplitude encodings. Our central observation is that the Gevrey hierarchy, which stratifies the space between smooth and analytic functions, provides a natural class for high-precision quantum algorithms. We use Gevrey regularity to analyze a quantum algorithm for solving general linear partial differential equations (PDEs) with periodic boundary conditions. We then extend the solver to PDEs with Dirichlet boundary conditions using thin Gevrey collars and a regularized inverse. Our resource bounds specify precision-dependent Fourier cutoffs and account for restricted discrete stability, elementary gate costs, and restriction to the Dirichlet domain. Finally, we apply our approach to isotropic and anisotropic Poisson equations with periodic and Dirichlet boundary conditions.

We establish sharp asymptotic rates for the two Barzilai–Borwein (BB) rules on uniformly positive quadratics and local nonlinear problems. In finite dimensions, for either fixed rule and an arbitrary positive first step, the gradient root factor is bounded by $(b_0-a_0)/(b_0+a_0)$, where $[a_0,b_0]$ is the initially active spectral interval. Hence the worst trajectory factor is $c_H=(\kappa(H)-1)/(\kappa(H)+1)$. When $H$ has at least two distinct eigenvalues, matched initialization and a balanced endpoint trajectory attain this value. Under matched initialization, the same constant is the optimal uniform-envelope threshold. For bounded, self-adjoint, uniformly positive operators on Hilbert space, scalar spectral measures yield the corresponding active-support bound and optimal matched uniform-envelope threshold, including continuous endpoint spectrum. Finally, if the gradient is strictly Fréchet differentiable at a stationary point and its derivative is self-adjoint and uniformly positive, every $\gamma\in(c_*,1)$, where $c_*=(\kappa(A_*)-1)/(\kappa(A_*)+1)$, is a uniform local envelope rate for either pure BB rule. Every well-defined trajectory converging to the stationary point has error and gradient root factors at most $c_*$ and objective-gap root factor at most $c_*^2$. Over the class of objectives with prescribed distinct derivative endpoints $m_*<M_*$, matched endpoint trajectories for quadratic and $C^\infty$ genuinely nonquadratic examples in $\mathbb R^2$ attain these factors.