Papers

We prove a three-term asymptotic formula for the eigenvalues of the fractional Laplacian on the bounded interval $(-1,1)$. This improves the eigenvalue asymptotics of Kulczycki–Kwaśnicki–Małecki–Stós and Kwaśnicki, and confirms the conjectural $O_\alpha(n^{-2})$ remainder suggested by the numerical simulations of Kaleta–Kwaśnicki–Małecki. Moreover, we prove that the normalized eigenfunctions are bounded uniformly in the eigenvalue index $n$ and the fractional order $\alpha$. This settles the conjecture proposed by Kwaśnicki through numerical experiments. Furthermore, we prove that the $n$-th eigenfunction has exactly $n-1$ zeros in the interval $(-1,1)$ and every zero is simple, and hence there are exactly $n$ nodal domains. A key ingredient in the proof is an explicit representation of the eigenfunction.

Transform-based tensor products, including the T-product and its more general form, namely the higher-order tensor-tensor product, have become fundamental tools for multilinear data analysis in applications such as image processing, signal reconstruction, and robotics. While invertible transforms enable tensor computations to be carried out via matrix operations in the transform domain, the resulting storage and computational costs remain prohibitive for high-dimensional, higher-order tensors. To address this challenge, we integrate low-rank tensor decomposition techniques, specifically tensor train decomposition (TTD) and hierarchical Tucker decomposition (HTD), into transform-based multilinear algebra to improve computational and memory efficiency. In particular, we develop TTD- and HTD-based formulations for the T-product and its associated key algebra, such as block diagonalization and tensor singular value decomposition, by operating directly on the factor matrices or tensors of the decompositions. The framework is further generalized to the higher-order tensor-tensor product and applied to multilinear model order reduction problems. We demonstrate the effectiveness and efficiency of our framework with numerical examples.

We employ the machinery of linear relations to the study of optimization problems in linear algebra. We first show that the relational version of the pseudo-inverse can be realized through a generalization of the least-squares problem. This allows one to prove that the pseudo-inverse realizes the solution of certain relational optimization problems. Our main result is showing that a certain truncation of this pseudo-inverse defines a solution to a relational version of the classical low-rank approximation problem which recovers both the Eckart-Young Theorem and several optimization problems involving pairs of matrices and vector spaces.

An explicit sequence of tetrahedra generated by the longest-edge bisection algorithm is shown to degenerate. The example violates shape regularity and both the minimum- and maximum-angle conditions, demonstrating that arbitrary tie-breaking among longest edges does not guarantee nondegeneration.

Bond-based peridynamics provides a non-local framework for modelling fracture without requiring spatial derivatives of the displacement field. However, when spatially varying horizons are used together with non-uniform discretisations, the classical single-horizon bond-based peridynamics formulation leads to asymmetric interactions between material points. These asymmetric interactions violate balance laws and can introduce non-physical artefacts such as ghost forces and spurious wave reflections. In this work, we develop a variational formulation for bond-based peridynamics with spatially varying horizons. Starting from the Lagrange-d'Alembert principle, we derive the governing equations of motion and show that the dual-horizon peridynamics formulation emerges naturally from the variation of the internal energy. Building on this variational structure, we construct asynchronous variational integrators that allow different time step sizes in different regions of the domain. This is particularly useful for dynamic fracture simulations with local refinement, where small time steps are required only near regions of high resolution or expected crack growth. Numerical examples involving wave propagation, a pre-cracked plate under tension, and the Kalthoff-Winkler impact experiment demonstrate that the proposed framework removes spurious reflections caused by non-uniform horizons, preserves physically consistent fracture patterns, and achieves results comparable to uniformly refined simulations. At the same time, the asynchronous variational integrator reduces the number of internal force evaluations compared to the standard velocity-Verlet method. The proposed approach therefore provides a consistent variational foundation and an efficient time-integration strategy for bond-based peridynamic simulations with spatially varying horizons.

Kellogg's checkerboard interface problem is a classical benchmark for robust adaptive finite element methods. Its successful adaptive meshes are radial: they refine strongly toward the interface crossing but show no angular structure, despite the large contrast and the asymmetric solution. We explain this by proving that the singular solution $u(r,\theta)=r^\gamma\mu(\theta)$ satisfies the exact identities $\kappa|\nabla u|^2=\Lambda r^{2\gamma-2}$ and $\kappa|\nabla^2 u|_F^2=2(1-\gamma)^2\Lambda r^{2\gamma-4}$, with $\Lambda=\gamma^2\cos^2(\pi\gamma/4)$ and the Hessian taken separately in each quadrant. The point is what has disappeared: the right-hand sides depend on $r$ alone, although $\kappa$ and $u$ each depend on the angle as well. Combined with equal discretization-error distribution, this shows that the target element density is radial, so a correct mesh should display nothing but refinement toward the center, and a Kellogg mesh that is not radial is visible evidence that the computation is not following the coefficient-weighted local difficulty. The reading is specific to this benchmark: on a second interface problem the same estimator correctly produces a strongly material-biased mesh, with a computed element-count ratio of $3.934$ against the predicted $4$. A byproduct gives the benchmark constants in closed form, so the problem data can be generated from $\gamma$ alone at any precision.

Accurate simulation of three-dimensional time-harmonic wave propagation over many wavelengths requires control of phase error and efficient solution of large indefinite systems. We develop a three-grid solver based on the compact 27-point interpolated optimized finite-difference (IOFD) discretization. Its wavenumber-dependent stencil supports a fine-grid resolution of six points per shortest wavelength and an unshifted physical correction on the \(2h\) grid at only three points per shortest wavelength. The method retains unshifted IOFD operators on the \(h\) and \(2h\) grids, while a complex-shifted \(2h\)–\(4h\) auxiliary cycle preconditions a factorization-free iterative approximation of the coarse inverse. Restricting the shift to this auxiliary cycle preserves the propagative character of the coarse correction. Comparison with the outgoing Green function confirms phase and relative-amplitude accuracy on a sequence of meshes up to \revision{\(10240^3\)}, \revision{exceeding one trillion unknowns,} with the largest problem spanning approximately \revision{1704} wavelengths per coordinate. The same fixed solver configuration retains robust convergence across smooth, discontinuous, high-contrast, and geophysical velocity models and exhibits scalable parallel performance. In particular, a \revision{\(2560^3\)} problem spanning approximately \revision{425} wavelengths in each coordinate direction is solved in \revision{43.8 seconds} on just 64 NVIDIA A100 GPUs.

Three-dimensional time-harmonic Maxwell simulations generate massive complex indefinite systems whose mesh coarsening is strictly limited by phase accuracy. Although matrix-free finite element kernels utilize GPU throughput efficiently, standard multilevel solvers are ultimately bottlenecked by the memory and communication costs of exact coarse-grid factorizations. We present a fully matrix-free, factorization-free three-grid preconditioner for curl-conforming Nédelec discretizations with perfectly matched layers (PML) and optimally blended quadrature. The method employs an outer FGMRES to solve the unshifted fine-grid equation, while an intermediate-grid correction is computed by a fixed-work FGMRES preconditioned with a complex-shifted $2h$–$4h$ cycle. This strategically confines the complex shift to an auxiliary preconditioner, preserving the physical Maxwell operator. A local Fourier analysis derives the blended Maxwell branches and compatible edge transfers, identifying robust shift and Jacobi damping parameters. Validated against the analytical Maxwell Green tensor, our approach demonstrates extreme scalability: using a single solver configuration, both homogeneous and highly heterogeneous systems with approximately 10.89 billion complex edge unknowns are solved in 42.0–72.0 seconds on just 64 NVIDIA A100 GPUs.

The collision operator for the Boltzmann equation is a nonlinear nonlocal operator. When lifted in the extended 2-particle space, it is viewed as the projection of a collisional linear operator. In this paper, we propose a particle method that samples the post-collision relative velocity directly from the Green's function of this operator (the transition probability of the generated time-continuous Markov chain). The normalizing flow amortized sampling is then proposed to reduce the sampling complexity. The resulted method takes $O(N)$ each time where $N$ is the particle number, and conserves momentum and energy exactly. This method does not require the boundedness of the kernel and, more importantly, it allows learning the Green's function directly from the scattering data without selecting the kernel in a specified family.

This work derives convergence guarantees for mirror descent and proximal mirror descent algorithms when a logarithmic barrier is used as a distance-generating function. Standard approaches cannot be applied when the solution lies on the boundary, where the Bregman divergence blows up. We show that, in a specific setting, both methods enjoy an $O(\log k / k)$ rate, which is also tight. In addition, our contributions include: (i) a new technique for handling the blow-up; (ii) a resolution of a gap in the theory of relative smoothness; and (iii) a comparison of the proposed approach with interior-point methods.

Stone (Ph.D.\ thesis, Department of Operations Research, Stanford University, 1981) proved that every matrix in $U \cap Q_0$ is a $P_0$-matrix and conjectured that the same conclusion holds for the larger class $E_0^f \cap Q_0$ of fully semimonotone $Q_0$-matrices. Murthy and Parthasarathy [SIAM J.\ Matrix Anal.\ Appl.\ 16 (1995), 1268–1286] verified the conjecture for matrices of order up to $4 \times 4$, for $5 \times 5$ and $6 \times 6$ matrices under additional hypotheses, and for several special subclasses of arbitrary order, but the conjecture remains open in general. In this paper we prove that every $E_0^f$-matrix with positive determinant is a $P_0$-matrix, for matrices of arbitrary order $n$; our proof proceeds by induction on $n$, via an algebraic analysis of principal minors under principal pivotal transforms. We further exhibit a matrix $A \in E_0^f$ with $\det A > 0$ that fails to belong to $Q_0$, showing that the hypothesis $\det A > 0$ used in our theorem cannot, by itself, be deduced from membership in $Q_0$, and hence does not on its own yield a proof of Stone's conjecture. Stone's conjecture itself remains open.

We develop mixed finite element methods for the linear stress-gradient elasticity problem based on symmetric Taylor–Hood elements. We establish the stability of the symmetric Taylor–Hood pair on simplicial meshes in both two and three dimensions, thereby resolving the stability left open by Brezzi, Fortin, and Marini in 1993. Combing the Nitsche's method, we further construct finite element schemes for linear stress problem in general boundary condition. Numerical experiments confirm the theoretical results.

A companion paper \cite{ItkinDF2026} introduced the Diagonal Frog (DF) positivity-preserving schemes for anisotropic Fokker–Planck equations, advancing each directional substep by a Krylov-computed matrix exponential, which dominates the cost. Replacing that exponential by a rational map $r(\gamma L)$ reduces the substep to a banded solve, but the positivity argument no longer applies. We prove that for eventually exponentially positive generators the entrywise sign of $r(\gamma L)$ at large steps is decided by a single number, the value $r(\infty)$ taken on infinitely stiff modes. Nonnegativity holds above a computable threshold when $0\le r(\infty)<1$, and at most on a bounded interval, empty or vanishingly narrow in all our tests, when $r(\infty)<0$. The criterion rejects the Crank–Nicolson (trapezoidal) method, where $r(\infty)=-1$, and selects the subdiagonal Padé$(0,2)$ method, which is second order, L-stable and provably positive above an explicit threshold. The resulting DF-ADI scheme costs $O(N)$ per step, keeps the implicit factorized mixed derivative unchanged, is second order in space and time, and conserves discrete mass exactly. In the strong cross-diffusion regime, however, the directional factors demand a step larger than the mixed derivative permits, so the composite second-order scheme is only empirically positive there, and the criterion serves to discriminate the well-behaved multiplicative factors from the stabilizing-correction schemes rather than to guarantee positivity. Against the Krylov exponential it runs ten to thirty-two times faster at matched accuracy in our tests with the gain growing with the mesh. We extend the construction to the backward Kolmogorov equation and to jump-diffusion models.

Stochastic reaction networks are continuous-time Markov chain models for interacting populations, with applications in biochemistry, epidemiology, ecology, and related areas. We study finite-difference sensitivity estimation when a single estimator requires several nearby parameterized paths. Existing variance-reducing couplings are typically pairwise, so that repeated use is either inefficient or requires application-specific choices in multi-path settings. We introduce the multi-path stacked coupling (MSC), a space-time Poisson construction that jointly generates any finite collection of parameterized paths. Each pairwise marginal of MSC has the same law as the corresponding split coupling pair, allowing existing variance bounds to transfer directly; in finite-state settings, we also obtain first-order expansions for the mean and second moment of finite-difference numerators. We apply MSC in three settings of practical importance: estimating many first derivatives simultaneously, estimating a single first derivative using a wider finite-difference stencil, and estimating higher-order derivatives. Numerical experiments on a processive phosphorylation network demonstrate strong performance in each of the three application areas considered, consistent with the theoretical advantages of MSC: across all three applications, MSC achieves the smallest root mean square error (RMSE) among the methods considered over the tested computational budgets.

Gradient random walk (GRW) methods represent a solution's spatial derivative by weighted particles and recover the solution by cumulative summation. We examine their accuracy for the one-dimensional heat, scalar cubic reaction-diffusion, and Burgers' equations. Thirty-seed ensembles, paired reconstructions, and deterministic controls identify the sources of two apparent accuracy limits. For the heat equation, comparing cumulative bin sums at bin centers introduces a first-order alignment error, which the exact reflected-particle distribution predicts and bin-edge comparison removes. An identity for the expected squared error predicts the remaining floor and the particle count for a target root-mean-square accuracy. With unequal weights, a count chosen in advance gives an error within 1% of its prediction. Reaction-front position, speed, and shape errors decrease under refinement. For Burgers' equation, exact transformed boundary data reduce the deterministic Cole-Hopf error about 490-fold, while fixed-endpoint particle runs approach their deterministic reference.

Graph neural networks (GNNs) have emerged as a powerful framework for learning from graph-structured data. However, their efficient training remains challenging, particularly in distributed computing environments. This challenge arises from the use of message passing, which couples all graph nodes, leading to expensive optimization steps, high memory requirements, and substantial communication overhead. To alleviate these limitations, we propose a novel domain-decomposition (DD) variant of AG2m, an AdaGrad method enhanced with second-order curvature information and momentum, denoted by DD-AG2m. The proposed DD-AG2m alternates between AG2m optimization on the original (global) graph and AG2m optimization on the partitioned graphs. To incorporate global information at reduced cost, we further introduce a two-level variant (2DD-AG2m) that performs global optimization steps on a coarse graph obtained by randomly subsampling nodes within each subdomain. Numerical experiments spanning graph classification, node-level regression, and spatiotemporal forecasting tasks demonstrate that the proposed DD methods reduce the computational cost required to achieve the same predictive performance by a factor of 4-8. Moreover, for the fixed computational cost, they improve the predictive performance of GNNs by up to 22% compared with the baseline AG2m.

We consider the compression of kernel matrices on the unit interval $[0,1]$ by interpolets that have sufficiently many vanishing moments. We define a compression rule which discards most matrix coefficients without compromising the accuracy offered by the underlying discretization. Since interpolets can be scaled such that the compressed kernel matrices are well conditioned, we derive a fully discrete scheme that solves a kernel interpolation problem under consideration in linear overall complexity. We finally generalize this approach to the unit $n$-cube $[0,1]^n$ by means of the sparse grid combination technique. Numerical experiments are carried out to validate the theoretical findings.

Theoretical guarantees for kernel regression are typically formulated in terms of smoothness, but practical accuracy depends critically on how the design resolution compares with the target and kernel lengthscales. We develop a finite-sample theory for Matérn regression on periodic domains with quasi-uniform designs, covering noiseless interpolation and noisy kernel ridge regression. Our minimax result shows that accurate recovery requires a design dense enough to resolve the target lengthscale and sufficient information at that scale to overcome noise. We prove that Matérn interpolation additionally requires the design to resolve the kernel lengthscale: if the kernel is too short relative to the point spacing, an additional error remains even when the target is well resolved. For noisy Matérn regression, we derive a three-term fixed-ridge risk characterization consisting of target-scale bias, kernel-scale bias, and variance, and show that optimizing over the ridge parameter yields four distinct contributions. When the target is no more than twice as smooth as the kernel, choosing a kernel lengthscale longer than the target does not worsen the oracle risk, whereas choosing one too short can. Thus these resolution conditions determine when accurate recovery becomes possible, while smoothness determines how rapidly the error decreases thereafter.

This article establishes the quasi-optimality of adaptive nonconforming finite element methods for a class of optimal control problems involving Dirac measures governed by the biharmonic equation. The nonconforming Morley finite elements are employed for discretising both the state and adjoint variables. A modified right-hand side through a companion operator that maps Morley finite elements toaconformingspacehelpstoovercomethechallengeinhandlingpointsourcesontheright-hand side. A priori and a posteriori error estimates for the optimal control problems are derived. Further,optimal convergence rates for adaptive finite element methods are established using an axiomatic framework: by proving key properties such as stability, reduction, discrete reliability, and quasi-orthogonality. Numerical experiments for three types of optimal control problems are discussed extensively and they validate the theoretical results.

The Jin–Xin relaxation of a nonlinear hyperbolic conservation law introduces a relaxation parameter that controls the width of the internal layer resolving a shock; the discontinuity of the limiting conservation law emerges only in the singular limit as this width vanishes. Physics-informed neural networks (PINNs) use smooth network approximations and are therefore not well suited to this limit, while relaxation PINNs with a fixed parameter resolve only a single scale and cannot follow the multiscale transition toward the limiting solution. We propose the Jin–Xin relaxation gradual convergence method (JXRGCM), which treats the relaxation parameter as a continuation variable, annealing it to zero along a schedule and warm-starting each stage from the previous one, so that the approximation follows the relaxation profile through progressively sharper scales. Under the sub-characteristic condition we establish a stability estimate whose constant is independent of the relaxation parameter; combined with the relaxation limit, it yields for scalar conservation laws an $L^2$ convergence rate of $\mathcal{O}(\varepsilon^{1/4})$ toward the entropy solution. Numerical experiments on the Burgers equation, the shallow-water dam-break problem, and the Sod shock tube show that JXRGCM improves shock and rarefaction resolution compared with fixed-parameter relaxation PINNs and other physics-informed approaches.

This work presents a concise, unified stability theory of high-order decoupled \lan{implicit-explicit linear multistep (IELM)} methods based on the well-known consistent splitting technique for the incompressible Navier-Stokes equation. With the help of the recent semi-generating function approach and the global discrete energy analysis, one can establish the unconditional stability of a consistent splitting IELM method with respect to the $\ell^{\infty}(H^1)\cap \ell^{2}(H^2)$ norm if the associated implicit-explicit controllability intensity is larger than $\sqrt{2}/2$, a constant determined by the Stokes pressure estimate. It is shown that the $\beta$-parameterized GBDF-$\rmk$ ($2\le \rmk\le5$) schemes and $\gamma$-parameterized SIELM-$\rmk$ ($2\le \rmk\le9$) schemes can fulfill this requirement of implicit-explicit controllability intensity by choosing proper parameters so that they can theoretically maintain the unconditional stability of the associated consistent splitting IELM methods. Numerical experiments are also included to support our theory.

We propose and analyze an IMEX-SDIRK2 mean-reverting concurrent-correction scalar auxiliary variable (mr-ccSAV) method for the forced two-dimensional periodic Navier–Stokes equations in vorticity form. The viscous term is treated by Alexander's SDIRK2 method and advection explicitly. Each stage requires two elliptic solves with the same shifted Laplacian and the solution of either a cubic or a linear scalar algebraic equation. For initial vorticity in $\dot L^s(\Omega)$, $s>2$, a stage solution exists for every positive time step; uniqueness is established separately under an explicit small-step condition. The principal result is a direct, unconditional uniform-in-time enstrophy bound for arbitrary positive time steps. For persistently bounded forcing, this estimate is absorbing: the influence of the initial data decays, and the forcing contribution does not accumulate in time. Under additional regularity, uniformly bounded step sizes, and bounded neighboring step ratios, we also establish uniform-in-time $H^1$ and $H^2$ vorticity bounds without a small-step condition. For smooth solutions, the method converges optimally at second order. Numerical experiments confirm its accuracy, long-time robustness, and effectiveness of a companion embedded time-step selector.

In this paper, we propose an extension of the well-known Robinson implicit function theorem for generalized equations from the upper Lipschitzian case to the upper Hölderian case. Explicit exponents dependence between the generalized equation and its linearization is determined. Applications to ball constrained least squares problems, including linear least squares and separable nonlinear least squares, are studied. In particular, we establish that the solution mapping of ball constrained linear least squares under linear perturbation is locally upper Hölder continuous with exponent $1/3$, which is of independent interest. Ultilizing the upper Hölderian version of the implicit function theorem, we show the local upper Hölderian of the solution mapping of parametric ball constrained linear least squares.

The smallest and largest eigenvalues of the Gram matrix induced by the differential neural tangent kernel (DNTK) play a pivotal role in the analysis of over-parameterized PINNs trained by gradient type algorithms. However, a theoretical analysis of the extreme eigenvalues remains completely absent due to the challenge posed by the presence of multiple differential operators. In this work, we provide explicit lower and upper bounds for the extreme eigenvalues of the infinite DNTK matrix for the Poisson equation with the Dirichlet boundary condition for two-layer RePU neural networks without the bias term. The setting is fairly general with respect to the sampling points and input dimension \(d\): \(\delta\)-separated and additionally \(d\geq 3\) when deriving the lower bound of the smallest eigenvalue. These results extend that for the neural tangent kernel, and to the best of our knowledge, represent the first results on the spectrum of the DNTK.

In this paper we focus on the stability of positive semidefinite matrix sensing maps $\Phi_{\mathcal{A}}(X)=(\langle A_i,X\rangle)_{i=1}^m$ where $A_i\succeq 0$, $X\succeq0$ and $\operatorname{rank}(X)\le r$. We introduce the bi-Lipschitz constants of $\Phi_{\mathcal{A}}(X)$ and define the global condition numbers as the ratio of upper and lower Lipschitz constants. We give deterministic universal lower bounds for these condition numbers, that depend only on the rank $r$ and on the underlying field. We then investigate the random rank-one Gaussian measurements and show that our lower bounds on condition numbers are asymptotically sharp, and therefore the random rank-one Gaussian measurements are asymptotically optimal. As an application, we derive the stability guarantees for an $\ell_1$-residual PhaseLift-type estimator at the optimal sampling scale.

The low-rank alternating direction implicit (ADI) method is an efficient solver for large-scale Stein equations with low-rank solutions. This paper shows that, as in the continuous-time Lyapunov equation case, the low-rank Cholesky factor ADI (LRCF-ADI) method for Stein equations implicitly performs $\mathcal{H}_2$-pseudo-optimal model order reduction for discrete-time systems. This observation leads to an automatic shift-generation strategy, allowing LRCF-ADI to select subsequent shifts without user intervention. The standard LRCF-ADI method requires shifts outside the unit circle. We generalize the method to allow shifts anywhere in the complex plane, including on the unit circle. This extension enables numerical integration for frequency-limited Stein equations by interpolating the integrand at points on the unit circle. It also enables non-intrusive, data-driven balanced truncation and frequency-limited balanced truncation using experimentally measurable transfer function samples on the unit circle, without requiring access to a state-space realization. Numerical results for large-scale models demonstrate the effectiveness of the proposed methods as low-rank Stein equation solvers and data-driven model order reduction methods.

This work presents a mathematical and computational framework for computing the dispersion relations of shear horizontal (SH) waves in continuously varying heterogeneous layered structures. The approach isolates the contribution of the heterogeneous layer from the complete dispersion relation, learns this contribution using a physics-informed neural network (PINN) and subsequently incorporates the trained model to determine the complete dispersion relation. The mathematical properties of the discretized problem are investigated, including the singularity of the finite-difference system and the oscillatory behavior of the layer solution, while a generalization-error estimate is established for the PINN approximation. The framework is first tested on a seismological configuration consisting of a heterogeneous sandstone layer over a granite half-space, where exponential heterogeneity is considered with independent variation rates in shear modulus and density. The proposed approach is validated against analytical solutions in special cases, while for general configurations the Haskell matrix method demonstrates convergence toward the continuously varying dispersion relation predicted by the PINN as the number of homogeneous sublayers increases. Parametric studies further confirm consistency with the underlying physics. An important feature of the method is that the heterogeneous-layer equation can be trained independently and reused in multiple settings. To demonstrate this, the same trained network is coupled with piezoelectric and piezomagnetic substrates governed by fundamentally different physical laws, highlighting the potential of the PINN framework as a reusable computational module for dispersion analysis in heterogeneous layered media.

Parameter estimation is a central task in data-driven learning of dynamical systems. It aims to recover the underlying physical parameters from observed time-series data, thereby providing interpretable insights into the physical mechanisms governing the system. Gradient/derivative matching methods based on Gaussian process provide an efficient way to perform parameter estimation. Those methods avoid repeated numerical integration and enforce local derivative consistency. However, such local matching may result in global inconsistency with the governing flow map, particularly under scarce and noisy observations. To address this limitation, we propose a framework based on Gaussian process learning with flow map refinement (GPL-FMR), a two-stage parameter estimation framework. The first stage is based on Gaussian process learning algorithm and the posterior obtained from which is transferred as an informative prior to the second stage based on flow-map refinement. The second stage further improves the parameter estimation via optimisation based on global dynamical constraints. We demonstrate and analyse its performance on multiple numerical examples, including the Van der Pol oscillator, the Lotka-Volterra model, and the Lorenz-63 system. The results show that the proposed framework consistently improves parameter estimation accuracy, particularly under scarce and noisy observations.

We analyze coefficient design in a one-parameter family of explicit palindromic electric–magnetic splittings for the time-domain Maxwell equations. After fourth-order staggered spatial discretization, the Fourier amplification matrix depends on the single scalar $g_2=a(1-2a)/2$. We prove that $a=1/4$ is the unique real coefficient maximizing the spectral CFL interval, with threshold $s_*=12/(7\sqrt d)$. We then identify a real-coefficient obstruction to higher phase accuracy: cancellation of the leading temporal phase defect requires $g_2=1/12$, whereas every real member satisfies $g_2\le 1/16$. The resulting complex-conjugate coefficients give fourth-order temporal phase accuracy for each fixed semidiscrete Fourier mode and have threshold $6\sqrt3/(7\sqrt d)$, while the complete field update remains globally second order in time. For real Maxwell data, the physical output is the real projection of the complex trajectory; this projection is branch independent and preserves the second-order error bound. We further give an exactly equivalent doubled real-arithmetic realization, which clarifies the role of the auxiliary imaginary component without changing the numerical method. A semidiscrete convergence result and numerical experiments confirm the distinction between stability optimization and phase optimization.

We study the natural superconvergence points and asymptotic expansions of one-dimensional spline finite element approximations. For a spline space of degree $k$ and any smoothness $0\le\mu\le k-1$, we prove that the $s$-th derivative of the error exhibits enhanced convergence of order $O(h^{k+2-s})$ at points where $k-s$ is even, provided the mesh is symmetric within a region of size $Ch|\ln h|$ around the point. This condition is known to be optimal for the cases of low derivative order $s=0,1$; the present analysis shows that the same local condition is sufficient for all admissible $s$. Moreover, by expanding the error in Legendre polynomials, a closure theorem determines the leading-order Legendre coefficients (the asymptotic expansion of the error) by combining the Galerkin orthogonality with the superconvergence conditions. For $\mu=k-1$ (B-splines) and $\mu=k-2$, the Galerkin orthogonality conditions vanish and the coefficients are determined solely by the superconvergence conditions. The asymptotic expansion can be expressed through a simple antiderivative recurrence on Legendre polynomials. The resulting polynomial's zeros encode the complete set of superconvergence points for all derivative orders. Numerical experiments for selected $(k,\mu)$ pairs confirm the theoretical predictions.