Papers

We consider boundary integral methods for Stokes mobility and resistance problems involving smooth axisymmetric particles. A primary numerical challenge is the accurate and efficient evaluation of layer potentials at off-surface points close to particle surfaces. We present a tolerance-driven workflow for evaluating the Stokes double layer potential at such evaluation points (targets) to prescribed accuracy while avoiding unnecessary computational cost. For each target–particle interaction, a fast classifier selects the least costly option estimated to meet the tolerance among standard, upsampled, and special quadrature. Geometry-dependent unit-density error indicators are precomputed and tabulated in reduced cylindrical coordinates, then combined on the fly with a local layer-density modifier, making its cost negligible relative to evaluating the potential. Targets requiring special quadrature are treated using a stabilized version of singularity swap surface quadrature: the periodic azimuthal integral is evaluated first using translated singularity swap quadrature to prevent severe cancellation near the surface, followed by adaptive Gauss–Legendre quadrature in the polar direction guided by error indicators. We integrate this workflow into a boundary integral solver with precomputed quadrature by expansion for on-surface self-interactions and demonstrate the workflow's performance for challenging configurations of spheroidal particles. Numerical results show that target classification is highly accurate. The prescribed tolerance is met for nearly all target–particle interactions and the error remains within a modest factor of the tolerance in the few remaining cases. Although the experiments focus on the Stokes double layer potential for spheroids, the off-surface framework applies to general smooth axisymmetric surfaces and can be easily adapted to other Stokes layer potentials.

We discuss some theoretical results and their applications to specific practical problems from wireless communications. We assume that we know the noisy version of the signal, which is sparse with respect to a given system of elements (dictionary), at a finite number of points and we want to approximately recover it. This problem of recovery of a noisy signal is closely related to the problem of establishing the Lebesgue-type inequalities for the corresponding algorithms and it motivates us to prove such inequalities. Under certain conditions on a dictionary (RIP-type condition, coherence condition) we obtain different kinds of the Lebesgue-type inequalities for the OMP and its version WOMP algorithms. The most important feature of our new theoretical result is the assumption that the dictionary has the RIP-type property instead of the assumption that it is the Riesz basis, which was used in the previous results. Our approach allows us to treat redundant (overcomplete) systems, which is important in applications. We consider the setting of sparse recovery for highly-coherent dictionaries that appear in OFDM setting for wireless communication. Our main example is the oversampled Fourier dictionary and the recovery of the frequency response of a sparse channel. We develop a general approach to such problems and we prove that under certain conditions the OMP algorithm recovers all dictionary elements with large coefficients, and estimate its accuracy (in NMSE metric) in terms of signal-to-noise ratio.

Linear programming plays a fundamental role in computer science, with applications in optimization, formal verification, SMT solving, and numerous other domains. When exactness guarantees are required, numerical inaccuracies arising from floating-point arithmetic can compromise soundness, whereas exact rational arithmetic often incurs significant computational overhead. In this paper, we investigate the numerical stability of the simplex algorithm and establish conditions under which a precision-boosting floating-point implementation provably produces the same pivot decisions and final basis as an exact rational implementation, from which the final result is reconstructed and certified exactly. Our analysis shows that correctness depends on both arithmetic precision and a careful handling of numerical tolerances. Based on these results, we develop a tolerance-aware, precision-boosting simplex algorithm with formal correctness guarantees. Finally, we introduce a delta-complete termination criterion that allows the algorithm to terminate once certified upper and lower bounds on the optimal objective differ by at most a user-specified threshold delta, providing a certified, user-controlled optimality gap.

Tree-encoded partitionings of graphs are fundamental tools for the decomposition and approximation of graph signals. For the efficient approximation of such graph signals, we develop strategies based on $hp$-refinement by combining domain decomposition with an improved local approximation using polynomials of higher degree. In this way, from a given graph partitioning tree, a more efficient subtree is extracted in which the cost of the signal approximation is considerably reduced by still maintaining the same total error. To this end, we interpret the refinement process as a binary knapsack problem to determine an enhanced partitioning tree. We further study an a-posteriori strategy which prunes the partitioning tree by optimizing the polynomial degrees over the subdomains. To make polynomial basis systems accessible for general graphs or high-dimensional data, we propose local embeddings of graphs into low dimensional Euclidean spaces. We underpin the efficiency of our algorithms with extensive numerical tests which carefully assess the impact of the applied refinements and optimization strategies.

Bressan and Shen (2021) jointly with Bressan and Murray (2020) proposed a simple example of a Cauchy problem on $\mathbb{R}^2$, characterized by two key parameters, intended to yield nonunique weak solutions of the incompressible Euler equations. They suggest that a limiting process applied to two specific approximations of their example initial data yields two distinct solutions. Here we provide evidence coming from a carefully designed numerical procedure to support if, or when, this intriguing claim is valid. Some of the numerical challenges encountered in the present study are: (i) the proposed initial data is unbounded in a neighborhood of the origin, (ii) the proposed problem is posed on all of $\mathbb{R}^2$, (iii) several essential symmetries present in the true solution must be captured, (iv) the proposed example depends on two parameters which may be critically linked. We demonstrate that an extremely large computational domain is mathematically necessary to capture the core symmetries of the exact free-space problem. Accordingly we implement a nested-grid strategy specifically designed for its particular suitability and efficacy for the problem at hand. Our strategy also allows us to achieve extreme grid refinement in the region where solution accuracy is most critical. Within our nested-grid framework, we design an underlying discretization scheme that is high-order, robust, and fast. We consider a broad collection of allowable parameter pairs, but among all of them we find only one that may satisfy the suggested claim.

Direct evaluation of the defining product in the real Cayley–Dickson algebra $A_n$, of dimension $N=2^n$, has quadratic arithmetic complexity. This paper gives a uniform algorithm for multiplication using $O(N\log N)$ real arithmetic operations and $O(N)$ auxiliary storage. The algorithm reduces multiplication to the alternating product on the imaginary subspace, then evaluates that product by a two-call recursion over one fixed quadratic coefficient extension. For $n\ge1$, the resulting bilinear algorithm uses at most $(9n-15)2^{n-1}+10$ input-dependent real multiplications, and for $n\ge3$, the specified arithmetic schedule uses $(34n-83)2^{n-1}+50$ real operations in total. Under this counting convention, the quasilinear schedule uses fewer operations than direct multiplication for $N\ge16$ and than the uniform Cariow–Cariowa method for $N\ge32$. The algorithm is implemented in the MIT-licensed C11 library fastCD, with a NumPy-backed Python interface, and its results are checked against an independent implementation of the defining recursion. In single-core benchmarks against direct multiplication and the uniform Cariow–Cariowa method, the quasilinear implementation had the lowest mean time of the three at every tested dimension $N\ge32$, for both single and batched products, and was roughly $16$ times faster than direct multiplication at $N=1024$.

This paper investigates finite element approximations of the one-dimensional viscous stochastic Burgers equation with additive trace-class noise. For the $P_2$ finite element spatial semi-discretization, we derive strong error estimates that are optimal with respect to regularity in \(L^p([0,T] \times \Omega;H_{D}^{\alpha,q})\) for \(p,q\in[2,\infty)\) and \(\alpha\in[-1,0]\), as well as an almost regularity-optimal estimate in \(L^p(\Omega;C([0,T];L^\infty(\mathcal{O})))\). Furthermore, we derive weak error estimates for moments of both terminal $L^q$-norms and space-time $L^p(0,T;L^q(\mathcal{O}))$-norms, with weak convergence rates (nearly) twice the corresponding strong ones. For the fully discrete scheme, which combines the \(P_2\) finite element method in space with a drift-implicit Euler–Maruyama scheme in time, we establish a strong temporal convergence rate of order \(\tau^{1/2-\varepsilon}\) in a discrete analogue of \(L^p(\Omega;C([0,T];L^\infty(\mathcal{O})))\), under the condition \(\tau \leqslant h^2\). Numerical experiments are presented to illustrate the theoretical convergence rates.

We analyze the mean square error of stochastic Runge-Kutta integrators for overdamped Langevin dynamics whose potential is strongly convex outside a bounded region. A decomposition splits the local error into a mean-zero term and a smaller remainder, and the discrete Poisson equation turns their moments into a bound on the error of a time average. We carry this out for two stochastic Runge-Kutta integrators proposed by Yang & Wang (2026), both of strong order $\frac32$ and weak order $2$, which evaluate the gradient of the potential and no higher derivative. We show that their laws approach the law of the exact solution at second order in Wasserstein-1 distance, up to a logarithm, uniformly in the number of steps. For a test function with bounded derivatives up to third order, we prove that the mean square error over $N$ steps with step size $h$ is $\mathcal O ( \frac{1}{Nh} + h^4 )$, which is the optimal order in the discretization. Experiments measure the strong and weak orders and the sampling bias on a nonconvex potential in $\mathbb R^2$, and compare the integrators on a diffusion model of CIFAR-10.

This paper introduces a novel method for approximating the inverse of discrete Fourier integral operators (FIOs). Given an $N \times N$ matrix representation $K$ of an FIO, the proposed algorithm consists of two stages. In the offline stage, we first construct a butterfly factorization (BF) $\tilde{K}$ of $K$, which enables fast forward matrix-vector multiplication. We then construct a hierarchically semiseparable (HSS) approximation $\tilde{G} \approx G$, where $G = K^{*} K$, using fast applications of $\tilde{K}$ and $\tilde{K}^{*}$ to random matrices. Finally, we apply the ULV factorization to the HSS matrix $\tilde{G}$ to obtain an approximation $\tilde{F} \approx G^{-1}$. Combining these approximations yields an approximate inverse $K^{-1} \approx \tilde{F} \tilde{K}^{*}$. The offline stage has complexity $O(N \log^{2} N)$ for 1D problems and $O(N^{1.5} \log N)$ for 2D problems. In the online stage, the proposed method approximates $K^{-1} u$ for a given input vector $u$ with complexity $O(N \log N)$ for both 1D and 2D problems. The proposed method can be used either as a direct solver or as a preconditioner for iterative methods. Numerical results for 1D and 2D FIOs demonstrate the effectiveness of the proposed method.

We construct Scott-Zhang-type space-time quasi-interpolation operators on tensor- product meshes for parabolic settings. Specifically, these are linear and uniformly bounded projections from the classical parabolic energy space onto a space-time tensor-product finite element space that preserve discrete Dirichlet boundary conditions on the lateral space-time boundary. We apply our operator in the context of space-time first-order system least-squares finite elements for parabolic equations to establish a quasi-optimal method for inhomogeneous Dirichlet boundary conditions.

This paper is the first in a three-part series on the construction of polynomial splines with the highest order of smoothness over hierarchical T-meshes, referred to as $\PHtwoT$-splines. For splines of bi-degree $(d,d)$, we study suitable refinement conditions for the subsequent basis construction, which requires dimensional stability of the underlying spline space. We present two groups of examples, considering unrestricted hierarchical refinement and refinement without vanishable T $l$-edges, respectively. The first setting permits new edges without additional degrees of freedom. In the second group, every refinement level excludes vanishable T $l$-edges and increases the dimension. Nevertheless, the dimension is unstable in both groups. We then introduce template translations to describe each refined region as a union of translates of a fixed template in the cell-index grid. Together with the known stability result under $(d-1)\times(d-1)$ template refinement, these examples justify this condition as a reasonable mesh assumption for the subsequent $\PHtwoT$-spline construction.

We stabilize stationary states of coupled multi-species McKean–Vlasov equations by feedback control, extending a spectral approach based on the Wasserstein Hessian and a Riccati equation from a single density to a product of Wasserstein spaces. We give necessary and sufficient conditions for local exponential stabilization at a prescribed rate when only one species is actuated, together with the minimum number of control inputs required. For two species, the condition reduces to a rank condition on the Fourier coefficients of the cross-interaction kernel over the eigenspaces of the unactuated species' linearized dynamics. A two-community noisy Kuramoto model and a cell-sorting model illustrate the results.

We introduce an amortized neural sampler that combines operator learning with flow methods for sampling. It maps SDE coefficient functions to pushforwards from a reference measure to the invariant measures, enabling efficient sampling across families of stochastic differential equations. Our framework shifts traditional sampling cost to an initial training phase, after which new SDE instances require only one encoder pass and a few ODE solver steps, independent of mixing time. To handle problems in high dimensions, we use Lagrangian trajectory sensors for the coefficient functions and cross attention in the architecture. We also theoretically establish the expressivity and resolution invariance of our framework. Experiments on 1D and 2D SDE families show competitive accuracy with substantial speedups over MCMC in regimes with slow mixing, transfer across sensor counts, and demonstration results on a 64D interacting particle SDE where traditional grid approaches are infeasible.

This paper describes Krylov subspace methods based on the tensor t-product for approximating quantities associated with third-order tensor functions. Specifically, it introduces the tensor-tubal Arnoldi method, which projects highdimensional tensor problems onto lower-dimensional tensor Krylov subspaces. The derived method exploits novel algebraic properties of the t-product to address large-scale tensor computations efficiently. Applications include the evaluation of parameter-dependent tensor functions and the solution of multidimensional ordinary differential equations. Numerical experiments illustrate the effectiveness and accuracy of the proposed methods compared with other approaches.

This article introduces four parallel numerical techniques for computing entries of the inverse of large sparse symmetric systems, all grounded in incomplete $LDL^T$ (ILDL) factorizations: (1) the selected inversion method (SelInv), which applies the $LDL^T$ factorization to recover entries of the matrix inverse within the sparsity pattern of the computed factors; (2) an approximate inversion method based on a truncated Neumann series expansion applied to the inverse of the L factor (NInv), providing an alternative at the cost of reduced accuracy; (3) a Mix approximation that merges the best of both SelInv and NInv; and (4) Mix-SPAI, which applies sparse approximate inverse (SPAI) refinement on the output of the Mix method to improve entry-level accuracy. To further improve accuracy while maintaining a stable sparsity pattern, we additionally employ a low-rank correction based on eigenvector updates, providing an alternative to tightening the drop tolerance. We report the performance of the proposed numerical techniques on a comprehensive collection of sparse matrices from scientific and industrial applications.

In this paper, we consider the classical time-optimal elevator problem in the presence of higher-order state constraints. The main contribution is a constructive indirect solution framework based on a Pontryagin Maximum Principle specifically formulated for higher-order constrained systems. First, we derive specialized optimality conditions and analyze the resulting structure of extremal trajectories. Second, we show that the infinite-dimensional optimal control problem can be transformed into a finite-dimensional system of nonlinear algebraic equations involving switching times, boundary-contact times, and multiplier parameters. This provides a computationally efficient procedure for calculating candidate optimal trajectories using standard nonlinear equation solvers. Third, we characterize the geometry of boundary contacts and explain the emergence of contact-chattering phenomena in higher-order constrained systems. In particular, we show that normal extremals cannot evolve along the constraint boundary for a positive amount of time, although sequences of boundary contacts may accumulate indefinitely. The proposed approach is demonstrated on a fourth-order elevator model, and the computed solutions are independently verified using a direct IPOPT-based optimization method.

For a wide class of Helmholtz scattering problems (namely, those fitting in the black-box framework of Sjöstrand–Zworski) we prove that the error incurred by approximating the radiation condition by a Cartesian perfectly-matched layer (PML) decreases exponentially in the width of the PML, the strength of the PML scaling, and the frequency.

We study early stopping for a continuous-time model of stochastic gradient descent (SGD) in ill-posed linear inverse problems. We consider an a posteriori stopping rule based on the discrepancy principle, which stops the dynamics once the residual reaches the noise level. Unlike deterministic gradient flow, the stochastic dynamics exhibits persistent multiplicative fluctuations whose amplitude scales with the residual itself. Under a source condition on the initial error and a polynomial bound on the spectrum of the empirical covariance operator, we prove that, with high probability, the error at this random stopping time achieves the minimax rate over the corresponding source class, up to a logarithmic factor. Our results therefore establish the discrepancy principle as an adaptive regularization strategy for continuous-time SGD that is minimax-optimal up to a logarithmic factor, despite the persistent fluctuations induced by stochastic sampling. To our knowledge, this is the first convergence-rate result for a continuous-time model of stochastic gradient descent in inverse problems under an a posteriori stopping rule: existing analyses, including those of variance-reduced variants, bound the stopping time but do not quantify the error attained there.

In a recent paper (H. Wang, Math. Comp. 93, 2861-2884, 2024), Wang has presented a number of results concerning the convergence of approximations based on generalised Laguerre polynomials, claiming to have provided "the first rigorous proof of root-exponential convergence of Laguerre approximations for analytic functions". In this note we argue that the proofs of the main results of Wang's paper are incomplete, because they rely on taking a limit under an integral sign, and this is not properly justified. We explain how the proofs could be completed using the dominated convergence theorem, provided that certain nontrivial inequalities involving generalised Laguerre polynomials and their weighted Cauchy transforms hold. We make a related conjecture, which, if true, would complete the proofs of Wang's results. This conjecture remains unproven, but appears plausible, based on numerical experiments.

A three-dimensional multiphysics finite element model for the simulation of high-frequency induction welding of tubes is presented. The model couples a time-harmonic $h-\phi$ formulation with a stabilized quasi-stationary advection-diffusion heat transport equation, enabling accurate prediction of electromagnetic and thermal fields under industrial operating conditions. The model incorporates parameterized geometry generation and semi-automated, physics-tailored mesh construction and is implemented using the open-source tools GetDP and Gmsh. Validation against measurements from a commercial induction welding line for AISI 304 stainless steel tubes demonstrates good agreement with operating data. The validated model is subsequently applied to investigate the influence of impeder material by comparing a conventional FeNiZnV ferrite with the soft magnetic composite Ferrotron 559H for the induction welding of AISI 304 stainless steel tubes.

We consider the problem of finding sparse solutions of an underdetermined linear system $Ax=b$. In contrast to conventional approaches based on greedy algorithms or convex relaxation, we reformulate sparse approximation as a structured system of polynomial equations and connect with the literature on tensor methods. We develop an eigenvalue decomposition based method that formally guarantees recovery of all sparse solutions if there is more than one. We also develop two optimization-based methods achieving favorable computational complexity. The new methods allow explicit control of the target sparsity. Numerical experiments illustrate the performance and compare to basis pursuit (denoising) and orthogonal matching pursuit.

We consider a class of sixth-order elliptic partial differential equations in two and three dimensions subject to simply supported and Cahn–Hilliard-type boundary conditions. Based on the Ciarlet–Raviart reformulation, we recast the sixth-order problem as an equivalent mixed system involving second- and fourth-order equations and establish its well-posedness under suitable assumptions. For the resulting mixed formulation, we propose and analyse two discretization frameworks: a fully non-conforming hybrid high-order (HHO) method on general polytopal meshes and a $C^0$-conforming HHO–finite element method on simplicial meshes. We prove stability and derive optimal-order a priori error estimates for the primary variables under appropriate regularity assumptions. For the fully non-conforming HHO method, we further derive a reliable residual-based a posteriori error estimator. Numerical experiments confirm the theoretical convergence rates and illustrate the robustness of the proposed methods for a range of mobility parameters.

Hyperbolic partial differential equations (PDEs) with time delay are essential mathematical tools used to model numerous physical systems where wave propagation or oscillations depend heavily on their historical states. As these applications scale in physical complexity, efficient parallel computing techniques become essential; however, developing highly scalable parallel solvers for delayed PDEs remains a significant computational challenge. To address this gap, this study extends the Dirichlet-Neumann Waveform Relaxation (DNWR) method to a multi-subdomain framework explicitly designed for solving time-delayed hyperbolic PDEs in parallel. We advance the underlying mathematical framework and provide a rigorous convergence analysis for both one-dimensional and two-dimensional spatial configurations. A central theoretical contribution of this work is the establishment of finite-step convergence for the multi-domain DNWR algorithm. These theoretical guarantees are firmly corroborated through comprehensive numerical experiments. Furthermore, a detailed comparative analysis against alternative domain decomposition techniques namely NNWR, OSWR, and Classical SWR, highlight the advantages of the proposed method. Finally, this approach provides a robust, highly parallelizable solver capable of efficiently handling complex hyperbolic PDEs with time delay.

Hyperbolic PDEs with time delay are essential for modeling a wide range of real-world applications. Since closed-form solutions are often not attainable, developing accurate and efficient numerical methods becomes crucial. In this work a novel substructuring waveform relaxation method namely Neumann-Neumann Waveform Relaxation (NNWR) has been implemented for numerically solving hyperbolic delay PDEs in the case of non-uniform subdomain partitioning. The convergence behavior is first analyzed using Fourier analysis to obtain an estimate, followed by a demonstration of finite-step convergence through Laplace transform techniques. Various test cases are presented as numerical illustrations.

We introduce a class of additive reduced basis preconditioners designed to accelerate the iterative solution of large-scale linear systems arising from discretized parametrized PDEs. The main idea is to regularize the inherently singular reduced-order approximation by adding simple correction terms: either a scaled identity correction or a projected correction based on a basic preconditioner. This yields nonsingular preconditioners under explicit and easily checked conditions. The construction and application of the preconditioners are integrated into an FGMRES framework through an offline strategy that dynamically builds the reduced-basis component by proper orthogonal decomposition at each FGMRES step. We establish sufficient conditions for the nonsingularity of the preconditioners and derive error bounds for the preconditioned Richardson iteration. Numerical results for convection-diffusion, anisotropic vortex, Stokes, and Helmholtz problems are provided to verify the efficiency and convergence of the proposed ARB preconditioners. The method consistently converges in a few iterations and substantially reduces online solve time, supporting its efficiency for multi-query engineering scenarios.

Diagonally implicit Runge-Kutta (DIRK) methods are a prominent class of numerical methods for solving stiff systems of ordinary differential equations (ODEs). Stiffness does not only impose stability challenges on Runge-Kutta methods; it can also degrade the order of convergence. This so-called order reduction phenomenon occurs when assumptions used for classical convergence analysis, e.g., an asymptotically small step size, fail to hold. In a prior paper by the authors, sharp order conditions and global error bounds for Runge-Kutta methods were developed, which hold uniformly with respect to stiffness when applied to a wide class of semilinear ODEs. In this work, those conditions are leveraged to construct the first DIRK methods of order four and five which satisfy these conditions and thus do not exhibit order reduction. Numerical results demonstrate that for a broad class of relevant nonlinear test problems, these new methods successfully mitigate order reduction, accurately estimate local error via an embedding for adaptive step size control, and can outperform classical DIRK methods.

This paper considers the lowest-order first type Nédélec edge element method (EEM) on Cartesian grids for the three-dimensional time-harmonic Maxwell equations with a large wave number. New polynomial preserving recovery (PPR) operators are proposed for the curl of the edge element solution and for the solution itself, respectively. Under the condition that $\kappa^3 h^2 C_{\mathrm{sol}}$ is sufficiently small, second-order superconvergence estimates are proved for both the recovered curl and the recovered solution, where $\kappa$ is the wave number, $h$ is the mesh size, and $C_{\mathrm{sol}}$ is a stability constant associated with the Maxwell solution operator. In particular, the analysis shows that the proposed PPR procedures cannot mitigate the well-known pollution effect inherent to the EEM. To reduce the pollution error, we further propose a new continuous interior penalty edge element method (CIP-EEM) that incorporates an additional normal-jump penalty term. It is shown that by appropriately choosing the penalty parameters, the new CIP-EEM can improve the phase error by two orders in $\kappa h$. Numerical experiments are presented to confirm the theoretical superconvergence results and to demonstrate that the CIP-EEM can effectively reduce the pollution error in the high-frequency regime.

We develop a method that uses Picard iterations to accelerate Newton continuation for the semi-implicit finite-difference discretization of forward-backward partial differential equation (PDE) systems arising in mean field games (MFGs). We first investigate the computational properties of the Picard and Newton methods, which are widely used separately in the MFG literature but whose comparative cost and robustness across different regimes remain insufficiently explored and documented. The Picard method uses an outer fixed-point iteration that alternates a forward Fokker-Planck solve and a backward Hamilton-Jacobi-Bellman solve. The Newton method instead applies Newton's method directly to the coupled nonlinear space-time system. Across one- and two-dimensional MFG benchmarks, Picard offers substantial computational savings in favorable regimes, but may fail at sufficiently low viscosity or require strong damping under temporal shocks, making Newton continuation preferable. With parameter continuation in the viscosity parameter, the Newton method is more robust in these regimes, at the cost of larger coupled linear systems. We relate these trade-offs to the residuals, Jacobian blocks, and sparsity structures produced by separable, local nonseparable, and nonlocal Hamiltonians. We then demonstrate how to combine inexpensive Picard iterations with Newton continuation in a hybrid method to reduce the total computational cost for a two- dimensional double-well problem.

Exploring high dimensional energy landscape is a challenging problem in many applications where both ground states and transition states offer important information about the underlying physical systems. To facilitate the computational exploration, we propose a cubic-regularized second-order method for locating critical points of prescribed index, including both energy minima (with an index $0$) and transition states (saddle points with positive index). At each step of the iteration, the proposed method involves several ingredients, such as updating the eigen-directions of the Hessian according to the prescribed index, reflecting the gradient and the Hessian along selected unstable eigen-directions, if any, and constructing the next step by solving a cubic-regularized subproblem. Under suitable assumptions, we show that, when the eigenspace and the cubic subproblem are solved exactly, the iteration converges locally and quadratically to a critical point of the prescribed index. We further establish convergence results in the inexact setting, where both the eigenspace computation and the solution of the cubic subproblem are performed approximately. We show that the resulting iteration retains local linear convergence, with the convergence rate depending on the accuracy of the inexact computations. Moreover, we prove that the fixed points of the algorithm are precisely critical points of the prescribed index. In addition, we introduce an adaptive strategy for updating the cubic regularization parameter to improve robustness and computational efficiency of the proposed algorithm. Numerical experiments confirm the predicted convergence behavior, demonstrate the effectiveness of the adaptive strategy, and illustrate the capability of the proposed method for locating critical points and constructing solution landscapes.

Stress localization in heterogeneous solids is governed by the bimaterial interface, where the displacement field remains $C^0$-continuous, while in-plane stresses jump due to the stiffness mismatch. Coordinate-based physics-informed neural networks (PINNs) represent this jump via a prescribed regularization width or a weighted interface penalty, making their accuracy sensitive to how phase-contrast changes are handled. This work presents a variational, label-free physics-informed graph neural network (PI-GNN) in which the heterogeneity is carried by the discretization rather than by the trial field. The solver operates on a conforming adaptive mesh graph, assigns constitutive behavior per element, and minimizes the discrete total potential energy as a single unweighted objective in which only first derivatives appear. The discrete energy on piecewise-linear elements coincides with the finite element (FE) Ritz functional. Dirichlet conditions are enforced by construction, with no penalty term, no interface weight, and no prescribed transition width. Using one fixed architecture, optimizer, and loss across small-strain elasticity and finite-strain Neo-Hookean hyperelasticity in two and three dimensions, the von Mises error remains below $3.58\%$ across a stiffness-contrast sweep spanning $(E_{\mathrm{inc}}/E_{\mathrm{mat}}\in[10^{-2},10^{2}])$, where a strong-form PINN degrades to $5.58\%$, and its displacement error reaches $7.66\%$ against $0.49\%$ for the PI-GNN. A trained network halves the ($\sigma_{xx}$) error of an energy-based PINN ($5.01\%$ versus $10.94\%$). Training cost exceeds a single FE solve by more than an order of magnitude, so the construction is a variationally consistent, penalty-free interface representation for parametric surrogates and inverse identification rather than a replacement for a one-off FE analysis.