1 Introduction↩︎

Physics-informed neural networks (PINNs) are commonly trained using sampled PDE residuals, boundary terms, and data constraints [1], [2]. However, the numerical values of the corresponding loss terms do not constitute a certificate for the predicted spatial field. In this work, we provide deterministic interval certifications for values of the exact PDE solution at specified query points. Such deterministic certificate must guarantee that the exact solution value is contained in the reported interval. A small empirical residual loss at finitely many samples does not provide such a guarantee; rigorous error control for PINN-type methods depends on stability, sampling, approximation, optimization, and loss consistency assumptions rather than on the sampled residual value alone [3], [4].

The inability to construct such certificate from loss values computed at fixed collocation points is structural. For an unprojected collocation PINN trained using finitely many automatic-differentiation residual samples, the field may be altered away from sampled locations without changing the sampled loss, while an unsampled point value may be changed arbitrarily. Finite collocation sampling therefore does not control unsampled point values deterministically, unless additional structure is imposed. 1 formalizes this limitation, and the compatibility condition introduced next provides the finite-dimensional structure needed for certification.

The certified model class is a mesh-based, piecewise-linear numerical PINN. The trainable model outputs an admissible coefficient vector in a finite-dimensional space. The resulting continuous field is the finite-element interpolation of that vector. The certifying residual that is used to derive the error bounds is the residual of the finite-dimensional numerical system applied to the same vector. This alignment is called compatibility. Compatibility is the central structural condition, since it makes the certified object well defined: the field, the residual, and the query all act on one admissible discrete state.

The error estimate is conducted as follows. First, for a compatible linear residual system, the error admits an exact discrete adjoint Green identity relative to a discrete numerical target. The signed correction recovers that target exactly, while norm-based, localized, randomized, and inexact variants yield computable approximate intervals when full correction is impractical. The transfer from the discrete target to the exact PDE solution is handled separately by rigorous comparison estimates. Finally, the same procedure extends to nonlinear residual systems. The stationary-target variant is given in 12. The linear square-system formulation assumes an invertible reduced operator after admissible reduction, covering the standard conforming Galerkin setting. A full-column-rank overdetermined extension is given in 11.2. Rank-deficient reduced operators are outside the certified scope.

The main contributions are as follows.

1. Compatibility is identified as the condition that makes deterministic pointwise certification well defined for a predicted numerical PINN field.

2. The square-system linear certificate is derived from the discrete target, through the exact identity and transfer step, to the final certified interval; a full-column-rank overdetermined extension is given in 11.2.

3. Discrete certification is separated from discrete-to-continuous transfer, so that the final error bound decomposes into distinct error contributions, namely: discrete correction or residual, query sensitivity, transfer, and representation mismatch.

4. A fully computable a posteriori transfer theorem is proved for a concrete coercive class, namely one-dimensional reaction–diffusion with conforming \(P_1\) Galerkin discretization, explicit residual estimators, and querywise constants.

5. The same procedure is extended to localized, randomized, and nonlinear certificates, with certified inexact-solve tolerances and a stationary-target extension provided in the appendices.

6. The theoretical claims are evaluated on manufactured pointwise benchmarks in two and three dimensions, on a public three-dimensional elasticity benchmark, and on projected neural load families, using scalar compatible queries, query-dependent adjoint Green sensitivities, corrected centers, raw-output intervals after projection, and comparison identities based on the released benchmark field.

For theoretical reasons stated in 1, the proposed method does not provide a certificate for unprojected collocation PINNs trained only on finitely many sampled residuals. The certificate applies after a compatible finite-dimensional numerical system is specified and the predicted field is represented within that system. Second, in a single linear instance, exact correction recovers the target of the chosen numerical discretization. The numerical PINN is therefore not presented as a single-instance substitute for a direct classical resolution, but as a source of admissible states that can be certified, approximately corrected, or reused across parametric, many-instance, and nonlinear settings.

The paper is organized as follows. 2 positions the work relative to the closest PINN, finite-element PINN, and goal-oriented finite-element literatures. 3 defines the certified object, the notation, and the scope assumptions. 4 gives the square-system linear certificate. 5 gives downstream computational variants and finite-element specializations, including a fully computable a posteriori transfer theorem for one-dimensional coercive reaction–diffusion. 11.2 gives the full-column-rank overdetermined linear extension. 6 gives the exact-root nonlinear certificate, while the stationary-target extension is given in 12. 7 provides an implementation procedure. 8 organizes the experiments by theorem-level claim. Additional experiments and reporting details are given in 10. Auxiliary analytical material is given in 11, and the proofs are given in 13.

2 Related Work↩︎

The problem that is addressed here is the computation of certified error bounds for PINN predictions. Recent works develop residual-based analysis of physics-informed models, and rigorous control of neural PDE solvers. Generalization, stability, and residual-to-solution estimates are developed in [3], [5]–[7], with Navier–Stokes extensions in [8]. For elliptic problems, [4] show that provable error control depends critically on loss consistency. Rigorous a posteriori certification has also begun to appear, including [9], [10]. These studies are closest in scope because they connect residual information to rigorous solution control. Other works, such as [11], leverage statistical learning theory to derive probabilistic error bounds, rather than deterministic guarantees. Finally, the authors of [12] derive an error decomposition for neural differential equation solvers, separating the contribution of the time discretization from that of the approximation error. However, their analysis focuses on neural ordinary differential equation solvers embedded within a numerical time-stepping framework, rather than on direct predictions produced by PINNs.

A second related line introduces variational, finite-element, or mesh-based structure into PINNs. hp-VPINNs, cPINNs, and related weak formulations build neural trial spaces into weighted-residual or domain-decomposition schemes [13], [14]. Quadrature and test-space effects are analyzed in [15], and reliable residual-type a posteriori estimators in the energy norm are given in [16]. Finite element informed neural networks (FEINNs) and their variants associate PINNs with increasingly explicit finite-element structure [17]–[21]. The Deep Ritz literature provides a related variational approach with its own approximation and certification results [22]–[24]. These works endow neural approximations with numerical structure. The present work is different in objective: it does not propose a new mesh-based neural solver, but a deterministic pointwise certification layer for a released admissible numerical field whose residual and query map are compatible with the same finite-dimensional system.

A separate issue is whether a mesh-based numerical PINN has a role beyond a direct finite-element solve. For a single well-posed coercive linear instance requiring high accuracy, no advantage over a direct finite-element solve is claimed; under exact correction, the certificate recovers the target of the chosen numerical discretization. The relevant regime is one in which a trainable model produces admissible states that are reused across loads, parameters, geometries, or queries, corrected only approximately, or used as local base points for nonlinear verification. In that regime, the finite-element structure is essential. It supplies the admissible coordinates, query map, and compatible residual needed to define which numerical object is certified. This object is not defined by an unprojected collocation PINN alone.

The transfer layer and the pointwise viewpoint draw on classical finite-element analysis rather than on PINN training theory. Goal-oriented and dual-weighted-residual methods for pointwise or scalar quantities of interest were developed in [25]–[27], while classical pointwise finite-element estimates date to [28]. These techniques serve a different purpose in the present setting. The approximate state is supplied by a numerical PINN, the discrete certificate is attached to the compatible residual of the same discrete system that defines the certified field. The adaptation to the continuous problem is completed only after a separate discrete-to-continuous transfer step.

The main contribution of our work, relative to the aforementioned studies, can be articulated along three complementary aspects. First, the certified quantity is defined through the compatibility between the field, the residual, and the query evaluated on a common admissible vector. Second, the passage from the discrete target to the exact PDE solution is formulated as an independent certified layer. Finally, our primary objective is the deterministic certification of scalar compatible queries, including pointwise evaluations, rather than the sole control of residuals or energy norms.

3 Framework and Compatibility↩︎

This section defines the scope used in this work to derive error bounds. First, the discrete, target solution is projected onto the exact PDE solution space. The error bounds are derived from this projection. An alternative construction based on direct continuous residual-to-point stability is given in 11.1.

3.1 Definition and objective↩︎

Let \(\Omega\subset\mathbb{R}^d\) be a bounded domain with closure \(\overline{\Omega}\). The goal is to certify the values of the exact PDE solution \(u^\star\) at prescribed query points \[\mathcal{Q}=\{x_1,\dots,x_K\}\subset\overline{\Omega}.\] For each query point, define the point-evaluation functional \[\lambda_i(u):=u(x_i), \qquad i=1,\dots,K.\] The desired deterministic certificate is an interval \(I_i\subset\mathbb{R}\) satisfying \[u^\star(x_i)\in I_i, \qquad i=1,\dots,K. \label{eq:deterministic95target}\tag{1}\] When randomized Green-sensitivity estimation is employed, the corresponding simultaneous certificate is \[\mathbb{P}_{\Xi} \left\{ u^\star(x_i)\in I_i\;\text{for all } i=1,\dots,K \right\} \ge 1-\delta, \label{eq:random95target}\tag{2}\] where the probability is taken only over the auxiliary randomization \(\Xi\). This probability statement is not a probabilistic model for \(u^\star\).

The certification problem therefore has two stages. A compatible discrete target is certified first; the resulting statement is then transferred to the continuous PDE solution through independently certified comparison bounds.

3.2 Certified numerical object and compatibility↩︎

The certificate is attached to one admissible coefficient vector and to one discrete query/residual system built on the same coordinates. The algebra is stated first for a general scalar compatible query. Nodal and off-grid \(P_1\) point values are then immediate special cases.

Let \(V_h\) be a finite-dimensional trial space, let \(\mathcal{U}_h\subseteq\mathbb{R}^{N_h}\) be the admissible coefficient set, let \[\mathcal{I}_h:\mathbb{R}^{N_h}\to V_h\] be a linear reconstruction map, and let \[F_h:\mathcal{U}_h\to\mathbb{R}^{M_h}\] be the discrete residual map.

For each prescribed query, assume there exists a vector \(\ell_{i,h}\in\mathbb{R}^{N_h}\) such that \[\Lambda_{i,h}(U):=\ell_{i,h}^{\top}U \qquad\text{for every }U\in\mathcal{U}_h. \label{eq:general95discrete95query}\tag{3}\] In the main pointwise setting, \[\Lambda_{i,h}(U)=(\mathcal{I}_hU)(x_i), \qquad i=1,\dots,K, \label{eq:point95query95special95case}\tag{4}\] Other scalar compatible queries can be assembled in the same basis.

This definition specifies the certified scope: admissible trainable models must output, or be explicitly mapped to, a coefficient vector in a specified finite-dimensional space. The certificate is attached to that vector and to the numerical system built on it, not to the training loss alone.

The linear results in 4 specialize to \(F_h(U)=A_hU-b_h\). The nonlinear results in 6 use the same compatibility requirement after reduction to admissible coordinates.

A mesh-based numerical PINN is any trainable map that outputs such an admissible vector \(U_\theta\in\mathcal{U}_h\). The certified field is \[u_{h,\theta}:=\mathcal{I}_hU_\theta. \label{eq:certified95field95general}\tag{5}\] If an unprojected neural field \(u_\theta\) is trained first, then an explicit extraction map \[\Pi_h:u_\theta\mapsto U_\theta\in\mathcal{U}_h\] must be specified. Common choices include nodal interpolation, projection, and constrained projection. The certificate is attached to \(u_{h,\theta}=\mathcal{I}_hU_\theta\), not to the unprojected neural field.

For the mesh-based \(P_1\) case, we use standard conforming simplicial finite-element notation [29]. Let \(\mathcal{T}_h\) be a conforming, shape-regular simplicial mesh of \(\Omega\) with nodal set \[\mathcal{N}_h=\{a_1,\dots,a_{N_h}\}.\] If a fitted mesh is not used and the computation is performed on \(\Omega_h\neq\Omega\), the resulting geometric consistency terms are absorbed later into the transfer bounds.

Let \[V_h^{P_1} := \left\{ v_h\in C^0(\overline{\Omega}): v_h|_T\in P_1(T)\;\text{for every }T\in\mathcal{T}_h \right\},\] with nodal basis \(\{\phi_j\}_{j=1}^{N_h}\). The reconstruction is \[\mathcal{I}_hU=\sum_{j=1}^{N_h}U_j\phi_j. \label{eq:p195reconstruction95map}\tag{6}\] Hence the point-query vector at \(x\in\overline{\Omega}\) is \[\ell_h(x):=\bigl(\phi_1(x),\dots,\phi_{N_h}(x)\bigr)^\top, \qquad (\mathcal{I}_hU)(x)=\ell_h(x)^\top U. \label{eq:p195query95vector}\tag{7}\] Consequently:

1. if \(x=a_k\) is a mesh node, then \(\ell_h(x)=e_k\) and the query is a nodal value;

2. if \(x\) lies in an element \(T=\operatorname{conv}\{a_{j_0},\dots,a_{j_d}\}\), then the nonzero entries of \(\ell_h(x)\) are exactly the barycentric coordinates of \(x\) in \(T\), so the query is an off-grid point value;

3. if a different scalar linear functional is needed later, its compatible query vector is assembled once in the same coefficient basis and used by the same algebra.

For the prescribed point set \(\mathcal{Q}=\{x_1,\dots,x_K\}\), define \[\ell_{i,h}:=\ell_h(x_i), \qquad i=1,\dots,K. \label{eq:query95vectors}\tag{8}\]

The compatible residual of the certified vector is \[r_h:=F_h(U_\theta). \label{eq:compatible95discrete95residual}\tag{9}\] In the linear setting this becomes \(r_h=A_hU_\theta-b_h\); in the nonlinear setting the reduced residual is introduced in 6. A sampled automatic-differentiation residual may replace 9 only if it is proved to equal the residual of this same compatible finite-dimensional system applied to an admissible vector.

If the final reported center is an unprojected neural output \(u_\theta(x_i)\) rather than the certified reconstruction \(u_{h,\theta}(x_i)=\ell_{i,h}^{\top}U_\theta\), we define the representation mismatch \[\zeta_{i,h}^{\theta} := \abs{u_\theta(x_i)-\ell_{i,h}^{\top}U_\theta}. \label{eq:representation95mismatch}\tag{10}\]

The compatibility chain is summarized in 1.

For linear systems, the primary formulation assumes square operators after admissible reduction; the full-column-rank overdetermined extension is given in 11.2.

For vector-valued problems, such as elasticity, the coefficient-level algebra is unchanged: replace \(V_h\) by the corresponding vector-valued finite-element space, stack the coefficients in \(U\), and choose \(\ell_{i,h}\) to extract the desired scalar component or scalar linear functional.

3.3 Finite collocation residuals evaluation do not control unsampled point values↩︎

The following proposition justifies the impossibility of constructing deterministic error bounds from finitely many collocation points, and therefore justifies the framework proposed here.

This proposition identifies a potential failure mode at the level of the underlying smooth problem, even if the specific perturbation \(u+\beta\varphi\) may not be representable by a given PINN architecture.

3.4 Continuous PDE target, transfer interface, and scope↩︎

The exact certification target is the continuous PDE solution \(u^\star\). The main theorems first certify a compatible discrete target and only then transfer that statement to \(u^\star\). Only the transfer interface required for the subsequent interval results is stated here; the alternative direct route based on continuous residual-to-point stability is given in 11.1.

Let \(V\) denote the admissible continuous solution class. For each reported query, the corresponding linear functional must be well defined on the class used in the transfer argument. For point queries in \(d\ge 2\), this is not provided by \(H_0^1(\Omega)\) alone and therefore requires additional regularity or a transfer argument stated in a stronger class. Smoothed scalar linear functionals, such as fixed-radius averages, do not require that pointwise boundedness assumption. Let \[\mathcal{L}:V\to Y_\Omega, \qquad \mathcal{B}:V\to Y_{\partial\Omega},\] and consider the boundary-value problem \[\mathcal{L}(u^\star)=f\quad\text{in }\Omega, \qquad \mathcal{B}(u^\star)=g\quad\text{on }\partial\Omega. \label{eq:continuous95pde}\tag{11}\]

The linear and nonlinear interval theorems below certify \(u^\star(x_i)\) by combining a discrete certificate for the compatible numerical target with a separate transfer term. That transfer term is supplied by an admissible comparison state together with rigorous residual and evaluation bounds; see 4.3 for the linear construction and 9 for the nonlinear exact-root analogue.

3.5 Main scope assumptions↩︎

The following assumptions determine the certified scope of the results.

These assumptions are part of the certified scope. They define the class of released numerical PINN fields for which the certificates apply.

4 Linear Deterministic Certification↩︎

This section presents the square reduced linear certificate: the compatible discrete target is defined, the exact residual-to-query identity is derived, the target is transferred to the continuous PDE solution, and the deterministic interval theorem is stated. Localized, randomized, and inexact linear variants are given in 5; the nonlinear generalization is given in 6. The proofs of the propositions and theorems of this section are given in 13.2.

4.1 Admissible coefficient space and discrete target↩︎

Let the discrete admissible set be the affine space \[\mathcal{U}_h:=U_{0,h}+Z_h\mathbb{R}^{n_h}\subseteq\mathbb{R}^{N_h}, \label{eq:admissible95affine95space}\tag{12}\] where \(Z_h\in\mathbb{R}^{N_h\times n_h}\) has full column rank. This representation covers strongly imposed boundary conditions, gauge constraints, quotient-space parametrizations, and nullspace removal. In the unconstrained case, \(U_{0,h}=0\), \(Z_h=I_{N_h}\), and \(n_h=N_h\).

Assume that the trained numerical PINN vector is admissible: \[U_\theta=U_{0,h}+Z_hy_\theta.\] If the unconstrained PINN vector is not admissible, it must first be replaced by a specified admissible projection, and the resulting computable mismatch must be included in the final interval.

For the linear discrete residual \[F_h(U)=A_hU-b_h, \qquad A_h\in\mathbb{R}^{M_h\times N_h}, \qquad b_h\in\mathbb{R}^{M_h},\] define \[\widehat A_h:=A_hZ_h, \qquad \widehat b_h:=b_h-A_hU_{0,h}.\] Then \[F_h(U_{0,h}+Z_hy)=\widehat A_hy-\widehat b_h.\] Let \(W_h\in\mathbb{R}^{M_h\times M_h}\) be symmetric positive definite and define the residual inner product \[\langle r,s\rangle_{W_h}:=r^\top W_h s, \qquad \norm{r}_{W_h}:=(r^\top W_h r)^{1/2}.\] Set \[\widehat H_h:=\widehat A_h^\top W_h\widehat A_h. \label{eq:reduced95normal95matrix95spd}\tag{13}\]

Different choices of \(W_h\) specify different residual geometries (Euclidean, quadrature-weighted, or energy-compatible) but in the square reduced setting they do not alter the exact identities below.

Define the discrete numerical target by \[y_h^\star := \widehat A_h^{-1}\widehat b_h, \qquad U_h^\star := U_{0,h}+Z_hy_h^\star. \label{eq:admissible95discrete95target}\tag{14}\] The compatible residual of the certified vector is \[r_h:=A_hU_\theta-b_h = \widehat A_hy_\theta-\widehat b_h. \label{eq:admissible95pinn95residual}\tag{15}\]

4.2 Exact query identities and sensitivities↩︎

For a discrete quantity-of-interest vector \(\ell_h\in\mathbb{R}^{N_h}\), we consider its restriction to the fixed admissible space. Define \[\widehat\ell_h:=Z_h^\top\ell_h\in\mathbb{R}^{n_h}. \label{eq:reduced95qoi95vector}\tag{16}\] Let \(q_{\ell,h}\in\mathbb{R}^{n_h}\) solve \[\widehat H_hq_{\ell,h}=\widehat\ell_h, \label{eq:reduced95adjoint95equation}\tag{17}\] and define the residual-side Green vector \[g_{\ell,h}:=W_h\widehat A_hq_{\ell,h}\in\mathbb{R}^{M_h}. \label{eq:reduced95residual95green95vector}\tag{18}\] Also set \[\sigma_{\ell,h} := \norm{g_{\ell,h}}_{W_h^{-1}}. \label{eq:general95sigma95definition}\tag{19}\] In the square reduced setting, \(g_{\ell,h}=\widehat A_h^{-\top}\widehat\ell_h\), independently of \(W_h\); the appendix gives the overdetermined full-column-rank interpretation.

For the query set \(\mathcal{Q}\), define the reduced query quantities by \[\begin{align} \widehat\ell_{i,h}&:=Z_h^\top\ell_{i,h}, &q_{i,h}&:=q_{\ell_{i,h},h},\\ g_{i,h}&:=g_{\ell_{i,h},h}, &\sigma_{i,h}&:=\sigma_{\ell_{i,h},h}. \end{align} \label{eq:query95green95quantities}\tag{20}\] Equivalently, \[\sigma_{i,h} = \left( \widehat\ell_{i,h}^{\top}\widehat H_h^{-1}\widehat\ell_{i,h} \right)^{1/2}. \label{eq:query95sigma95reduced}\tag{21}\]

4.3 Transfer to the continuous target↩︎

To transfer the discrete certificate to the continuous PDE solution, choose any admissible comparison vector \[U_h^\pi = U_{0,h}+Z_h y_h^\pi \in \mathcal{U}_h. \label{eq:admissible95comparison95vector}\tag{22}\] This vector need not arise from a canonical coefficient map of \(u^\star\). It may be an interpolant, a constrained lift, a projection, or any other admissible numerical representation for which certified comparison estimates are available.

Define its consistency residual by \[\tau_h^\pi := A_h U_h^\pi - b_h = \widehat A_h y_h^\pi - \widehat b_h. \label{eq:comparison95consistency95residual}\tag{23}\] Assume that certified bounds are available such that \[T_h\ge\norm{\tau_h^\pi}_{W_h}, \label{eq:certified95consistency95bound}\tag{24}\] and \[\varepsilon_{i,h}^{\mathrm{eval}} \ge \abs{\ell_{i,h}^{\top}U_h^\pi-u^\star(x_i)}. \label{eq:certified95eval95error}\tag{25}\] Define \[\mu_{i,h}:=\sigma_{i,h}T_h+\varepsilon_{i,h}^{\mathrm{eval}}. \label{eq:admissible95mu}\tag{26}\]

4.4 Deterministic pointwise intervals↩︎

5 Linear Variants and Specializations↩︎

This section gives computational variants for the square reduced setting: norm-only certificates, a standard coercive Galerkin specialization, localization, randomized sensitivity estimation, and inexact linear algebra. The effective-residual projector and the full-column-rank overdetermined extension are given in 11.2.

5.1 Norm-only certificates and width decomposition↩︎

In the square reduced setting, \[\abs{g_{i,h}^{\top}r_h} \le \sigma_{i,h}\norm{r_h}_{W_h}. \label{eq:projected95norm95bound}\tag{32}\]

Define the three radii \[\begin{align} \rho_{i,h}^{\mathrm{corr}} &:= \sigma_{i,h}T_h+\varepsilon_{i,h}^{\mathrm{eval}}, \tag{33}\\ \rho_{i,h}^{\mathrm{pred}} &:= \abs{g_{i,h}^{\top}r_h} + \sigma_{i,h}T_h + \varepsilon_{i,h}^{\mathrm{eval}}, \tag{34}\\ \rho_{i,h}^{\mathrm{norm}} &:= \sigma_{i,h}\norm{r_h}_{W_h} + \sigma_{i,h}T_h + \varepsilon_{i,h}^{\mathrm{eval}}. \tag{35} \end{align}\] Then \[\begin{align} \abs{c_{i,h}^{\mathrm{corr}}-u^\star(x_i)} &\le \rho_{i,h}^{\mathrm{corr}}, \tag{36}\\ \abs{\ell_{i,h}^{\top}U_\theta-u^\star(x_i)} &\le \rho_{i,h}^{\mathrm{pred}} \le \rho_{i,h}^{\mathrm{norm}}, \tag{37}\\ \abs{u_\theta(x_i)-u^\star(x_i)} &\le \zeta_{i,h}^{\theta}+\rho_{i,h}^{\mathrm{norm}}. \tag{38} \end{align}\] If only an a priori transfer rate \(T_h\le C_\tau h^\alpha\) is available, it can be directly substituted into 33 –35 .

The width has three mathematically distinct sources: \[\begin{multlined} \rho_{i,h}^{\mathrm{norm}} = \underbrace{\sigma_{i,h}\norm{r_h}_{W_h}}_{\text{discrete compatible error bound}}\\ + \underbrace{\sigma_{i,h}T_h}_{\text{discrete-to-continuous transfer}} + \underbrace{\varepsilon_{i,h}^{\mathrm{eval}}}_{\text{evaluation comparison}}. \end{multlined} \label{eq:norm95radius95source95split}\tag{39}\] The width also includes the additional term \(\zeta_{i,h}^{\theta}\) only when the reported center is an unprojected neural output rather than the certified reconstruction. For the signed prediction-centered interval, replace the first term in 39 by \(\abs{g_{i,h}^{\top}r_h}\).

A large interval radius may result from large sensitivity \(\sigma_{i,h}\), large discrete compatible error, large transfer residual \(T_h\), or large representation/evaluation mismatch. These contributions should be reported separately. The sharper overdetermined refinement based on the effective-residual projector \(P_h^{\mathrm{eff}}\) is given in 11.2.

5.2 Finite-element application for coercive Dirichlet elliptic problems↩︎

The abstract transfer term becomes explicit after \(T_h\) and \(\varepsilon_{i,h}^{\mathrm{eval}}\) are instantiated for a PDE class. The finite-element ingredients used below are standard conforming Galerkin theory, Céa-type quasi-optimality, and \(P_1\) interpolation estimates on shape-regular meshes [29]. Classical pointwise finite-element estimates and goal-oriented a posteriori estimators may supply pointwise transfer bounds [25]–[28]. A standard linear specialization is the conforming \(P_1\) Galerkin discretization of a coercive second-order Dirichlet problem.

Assume \(V=H_0^1(\Omega)\) and that the continuous problem is \[a(u^\star,v)=L(v) \qquad\text{for all }v\in V, \label{eq:coercive95dirichlet95weak95form}\tag{40}\] with \[a(u,v) := \int_\Omega \mathbf{A}(x)\nabla u\cdot\nabla v + c(x)uv\,\,\mathrm{d}x, \label{eq:coercive95bilinear95form}\tag{41}\] where \(\mathbf{A}\in L^\infty(\Omega;\mathbb{R}_{\mathrm{sym}}^{d\times d})\) is uniformly elliptic and \(c\in L^\infty(\Omega)\) satisfies \(c\ge 0\). Assume \[\alpha_a \norm{v}_{H_0^1(\Omega)}^2 \le a(v,v), \qquad \abs{a(u,v)} \le M_a \norm{u}_{H_0^1(\Omega)}\norm{v}_{H_0^1(\Omega)}.\] Let \(u_h^{\mathrm{FE}}\in V_h\) denote the conforming \(P_1\) Galerkin solution satisfying \[a(u_h^{\mathrm{FE}},v_h)=L(v_h) \qquad\text{for all }v_h\in V_h. \label{eq:galerkin95fe95problem}\tag{42}\] Write \[u_h^{\mathrm{FE}}=\mathcal{I}_hU_h^{\mathrm{FE}} = \mathcal{I}_h(U_{0,h}+Z_h y_h^{\mathrm{FE}}).\] After elimination of essential constraints, the reduced Galerkin matrix \(\widehat A_h\) is square and symmetric positive definite. Set \[\norm{v}_{a}:=a(v,v)^{1/2}.\]

Thus, in the exact-correction regime, the remaining continuous-transfer task is the pointwise finite-element error problem for the chosen PDE class. The certificate uses such bounds as external transfer inputs and does not alter the finite-dimensional discrete correction.

Let \(\mathcal{J}_h:C^0(\overline{\Omega})\to V_h\) denote the nodal interpolation operator and set \[u_h^\pi:=\mathcal{J}_h u^\star, \qquad U_h^\pi \text{ the coefficient vector of }u_h^\pi.\]

The next bound is an a priori instantiation of the transfer interface obtained by combining the comparison identity above with the standard finite-element estimates cited at the start of this subsection.

Asymptotic contraction criteria are presented in 11.3. They aid interpretation of the radii but are not needed for the main interval construction.

5.3 Localized residual certificates↩︎

Localization is query-dependent. In the square reduced setting, the exact discrete error is represented directly by \(r_h\), so the local/tail split is written in terms of the full compatible residual.

Fix a query point \(x_i\) and a radius parameter \(R\). Let \[\Pi_{i,R}:\mathbb{R}^{M_h}\to\mathbb{R}^{M_h}, \qquad \Pi_{i,R}^2=\Pi_{i,R},\] be any linear projector whose range is the subspace of residual rows designated as local to \(x_i\) inside radius \(R\). Define \[g_{i,h}^{\mathrm{loc}}(R):=\Pi_{i,R}^\top g_{i,h}, \qquad g_{i,h}^{\mathrm{tail}}(R):=(I-\Pi_{i,R})^\top g_{i,h},\] and \[r_h^{\mathrm{loc}}(R):=\Pi_{i,R}r_h, \qquad r_h^{\mathrm{tail}}(R):=(I-\Pi_{i,R})r_h.\]

Define the local and tail sensitivity masses by \[\chi_{i,h}(R):=\norm{g_{i,h}^{\mathrm{loc}}(R)}_{W_h^{-1}}, \qquad \vartheta_{i,h}(R):=\norm{g_{i,h}^{\mathrm{tail}}(R)}_{W_h^{-1}}. \label{eq:local95tail95sensitivity95masses}\tag{50}\] The quantity \(\chi_{i,h}(R)\) is a rigorous sensitivity index for the local residual block. The product \(\vartheta_{i,h}(R)\norm{r_h^{\mathrm{tail}}(R)}_{W_h}\) is the required deterministic tail contribution. Omitting the tail term removes the deterministic certificate and leaves only a heuristic indicator.

If \(\Pi_{i,R}\) is \(W_h\)-orthogonal, for example a coordinate selector when \(W_h\) is diagonal or block diagonal by residual cell, then \[\sigma_{i,h}^2 = \chi_{i,h}(R)^2+\vartheta_{i,h}(R)^2. \label{eq:orthogonal95local95tail95sigma95split}\tag{51}\] Indeed, the cross term vanishes because \[\bigl(g_{i,h}^{\mathrm{loc}}(R)\bigr)^\top W_h^{-1} g_{i,h}^{\mathrm{tail}}(R)=0\] whenever \(\Pi_{i,R}^\top W_h=W_h\Pi_{i,R}\).

The full-column-rank overdetermined version, where the exact discrete error depends on \(P_h^{\mathrm{eff}}r_h\), is given in 11.2.

5.4 Randomized sensitivity estimation↩︎

When the number of queries is large, exact computation of all \(\sigma_{i,h}\) may dominate cost. Random probing for diagonal and quadratic-form estimation is standard in randomized numerical linear algebra [30]. Here Gaussian probes are used because they give an exact \(\chi^2\) calibration for the scalar quadratic forms below. Gaussian probing estimates \[\sigma_{i,h}^2 = \widehat\ell_{i,h}^{\top}\widehat H_h^{-1}\widehat\ell_{i,h}.\] Let \(W_h^{1/2}\) denote the symmetric positive-definite square root of \(W_h\). Let \(z^{(t)}\sim\mathcal{N}(0,I_{M_h})\) be independent, and solve \[\widehat H_h \psi^{(t)} = \widehat A_h^\top W_h^{1/2}z^{(t)}, \qquad t=1,\dots,m. \label{eq:reduced95random95probe95solve}\tag{52}\] Then \[\operatorname{Cov}(\psi^{(t)})=\widehat H_h^{-1},\] and therefore \[\widehat\ell_{i,h}^{\top}\psi^{(t)} \sim \mathcal{N}(0,\sigma_{i,h}^2).\] Define \[\widehat\sigma_{i,h}^{\,2} := \frac{1}{m} \sum_{t=1}^m \left(\widehat\ell_{i,h}^{\top}\psi^{(t)}\right)^2. \label{eq:reduced95random95sigma95estimator}\tag{53}\] When \(\sigma_{i,h}>0\), \[\frac{m\widehat\sigma_{i,h}^{\,2}}{\sigma_{i,h}^2} \sim \chi_m^2.\] The zero-sensitivity case is immediate.

Let \(q_{\alpha,m}\) be the lower \(\alpha\)-quantile of \(\chi_m^2\): \[\mathbb{P}\{\chi_m^2\le q_{\alpha,m}\}=\alpha.\] For \(K\) query points, define \[\widehat\sigma_{i,h}^{+} := \left( \frac{m\widehat\sigma_{i,h}^{\,2}}{q_{\delta/K,m}} \right)^{1/2}. \label{eq:simultaneous95sigma95upper95bound}\tag{54}\] By the union bound, \[\mathbb{P} \left\{ \sigma_{i,h}\le\widehat\sigma_{i,h}^{+} \;\text{for all }i=1,\dots,K \right\} \ge 1-\delta. \label{eq:simultaneous95sigma95event}\tag{55}\]

For sharp simultaneous intervals, we compute the exact discrete correction once by solving \[\widehat A_h e_\theta = r_h. \label{eq:one95primal95correction95solve}\tag{56}\] Then \[U_\theta-U_h^\star=Z_he_\theta,\] and \[c_{i,h}^{\mathrm{corr}} = \ell_{i,h}^{\top}(U_\theta-Z_he_\theta). \label{eq:randomized95corrected95center}\tag{57}\] On the event 55 , the intervals \[\begin{multlined} I_i^{\mathrm{corr,rand}} :=\\ \left[ c_{i,h}^{\mathrm{corr}} - \widehat\sigma_{i,h}^{+}T_h - \varepsilon_{i,h}^{\mathrm{eval}},\right.\\ \left. c_{i,h}^{\mathrm{corr}} + \widehat\sigma_{i,h}^{+}T_h + \varepsilon_{i,h}^{\mathrm{eval}} \right] \end{multlined} \label{eq:randomized95corrected95interval}\tag{58}\] satisfy \[\mathbb{P} \left\{ u^\star(x_i)\in I_i^{\mathrm{corr,rand}} \;\text{for all }i=1,\dots,K \right\} \ge 1-\delta.\]

Without the correction solve 56 , the looser representation-centered interval is \[\begin{multlined} I_i^{\mathrm{norm,rand}} :=\\ \left[ \ell_{i,h}^{\top}U_\theta - \widehat\sigma_{i,h}^{+} \bigl(\norm{r_h}_{W_h}+T_h\bigr) - \varepsilon_{i,h}^{\mathrm{eval}},\right.\\ \left. \ell_{i,h}^{\top}U_\theta + \widehat\sigma_{i,h}^{+} \bigl(\norm{r_h}_{W_h}+T_h\bigr) + \varepsilon_{i,h}^{\mathrm{eval}} \right]. \end{multlined} \label{eq:randomized95norm95interval}\tag{59}\] Rademacher probes remain unbiased for \(\sigma_{i,h}^2\), but the exact \(\chi_m^2\) calibration above no longer applies; separate concentration bounds are required in that setting.

5.5 Inexact linear algebra↩︎

Certified allowance formulas for inexact correction are given in 11.3. In the square reduced setting, primal correction and adjoint solves should use \(\widehat A_h\) and \(\widehat A_h^\top\) directly, not the normal matrix \(\widehat H_h\). The full-column-rank overdetermined extension and its matrix-free least-squares solvers are given in 11.2.

6 Nonlinear Certification by Verified Local Roots↩︎

This section presents the exact-root certificate. A local discrete exact root is certified first; the resulting discrete target is then transferred to the PDE solution. The stationary-target variant is given in 12, and one explicit verified-ball construction for square systems is given in 11.3.

All nonlinear statements are written after admissible constraints have been eliminated. Fix a norm \(\norm{\cdot}_{X_h}\) on \(\mathbb{R}^{n_h}\); all local radii and operator norms below are measured with respect to \(\norm{\cdot}_{X_h}\) on the state space, \(\norm{\cdot}_{W_h}\) on the residual space \(\mathbb{R}^{M_h}\), and the dual norm \(\norm{\cdot}_{X_h^\ast}\) on \((\mathbb{R}^{n_h})^\ast\). The reduced unknown is \(y\in\mathbb{R}^{n_h}\), the corresponding admissible full vector is \[U(y):=U_{0,h}+Z_h y,\] and the query map at \(x_i\) is \[L_{i,h}(y):=\ell_{i,h}^{\top}U(y) = \ell_{i,h}^{\top}U_{0,h}+\widehat\ell_{i,h}^{\top}y, \qquad \widehat\ell_{i,h}:=Z_h^\top \ell_{i,h}. \label{eq:nonlinear95reduced95query95map}\tag{60}\] This is the reduced-coordinate form of the compatible query \(\Lambda_{i,h}\). Let \[F_h:\mathbb{R}^{n_h}\to\mathbb{R}^{M_h}\] be continuously differentiable on a neighborhood of the relevant points, and let the trained PINN induce \(y_\theta\in\mathbb{R}^{n_h}\) with residual \[r_h:=F_h(y_\theta).\]

6.1 Certified exact-root targets↩︎

A discrete exact-root target is a vector \(y_h^\star\) satisfying \[F_h(y_h^\star)=0.\] When \(M_h>n_h\), full column rank of \(DF_h(y)\) does not certify existence of such a root; exact-root certificates require a separately verified local zero.

For any base point \(y\) at which \[J_h(y):=DF_h(y)\] has full column rank, define \[H_h^{\mathrm{root}}(y):=J_h(y)^\top W_h J_h(y).\] For a reduced functional \(\widehat\ell_h\in\mathbb{R}^{n_h}\), let \(q_{\ell,h}^{\mathrm{root}}(y)\) solve \[H_h^{\mathrm{root}}(y)q_{\ell,h}^{\mathrm{root}}(y)=\widehat\ell_h,\] and set \[g_{\ell,h}^{\mathrm{root}}(y):=W_hJ_h(y)q_{\ell,h}^{\mathrm{root}}(y), \qquad \sigma_{\ell,h}^{\mathrm{root}}(y):=\norm{g_{\ell,h}^{\mathrm{root}}(y)}_{W_h^{-1}}.\] Equivalently, \[\sigma_{\ell,h}^{\mathrm{root}}(y) = \left( \widehat\ell_h^\top H_h^{\mathrm{root}}(y)^{-1}\widehat\ell_h \right)^{1/2}.\] For a deviation \(d\in\mathbb{R}^{n_h}\), define the base-point Taylor remainder \[R_h^{\mathrm{root}}(y;d) := F_h(y-d)-F_h(y)+J_h(y)d. \label{eq:base95point95root95remainder}\tag{61}\] If \(DF_h\) is Lipschitz on the ball \(B_{X_h}(y,t)\) with constant \(L_h(y,t)\) in the operator norm from \((\mathbb{R}^{n_h},\norm{\cdot}_{X_h})\) to \((\mathbb{R}^{M_h},\norm{\cdot}_{W_h})\), then \[\norm{R_h^{\mathrm{root}}(y;d)}_{W_h} \le \frac{L_h(y,t)}{2}\norm{d}_{X_h}^2, \qquad \norm{d}_{X_h}\le t, \label{eq:root95remainder95lipschitz95bound}\tag{62}\] because \[R_h^{\mathrm{root}}(y;d) = \int_0^1 \bigl(J_h(y)-J_h(y-sd)\bigr)d\,\,\mathrm{d}s.\]

Such neighborhoods can be certified by standard Newton–Kantorovich verification arguments [31], [32]. For square systems, 11.3 records the corresponding verified ball in the norms used here.

6.2 Discrete-to-continuous comparison for exact-root targets↩︎

6.3 Nonlinear exact-root certificate↩︎

7 Computational Procedure and Reporting for \(P_1\) Numerical PINNs↩︎

For the linear compatible \(P_1\) case, the certificate can be assembled as follows.

The following procedure applies to the square reduced setting.

1. Assemble the compatible system. Build the mesh \(\mathcal{T}_h\), the admissible set \(\mathcal{U}_h=U_{0,h}+Z_h\mathbb{R}^{n_h}\), the reconstruction \(\mathcal{I}_h\), and the compatible residual \(F_h(U)=A_hU-b_h\). Form \[\widehat A_h=A_hZ_h, \qquad \widehat b_h=b_h-A_hU_{0,h},\] and form \(\widehat H_h=\widehat A_h^\top W_h\widehat A_h\) only when it is specifically needed, for example in the randomized sensitivity procedure.

2. Represent the trained model by an admissible coefficient vector. Obtain \(U_\theta=U_{0,h}+Z_hy_\theta\in\mathcal{U}_h\), either directly or through an explicit extraction map \(\Pi_h\). Compute the compatible residual \[r_h=A_hU_\theta-b_h.\]

3. Assemble the query vectors once. For point queries, use \(\ell_h(x)=\bigl(\phi_1(x),\dots,\phi_{N_h}(x)\bigr)^\top\). Nodal values are the special case \(\ell_h(a_k)=e_k\). Off-grid point values use barycentric weights on the containing element. Any other scalar linear functional is treated by assembling its coefficient vector \(\ell_{i,h}\) in the same basis.

4. Choose the discrete correction regime. Use the direct square-system formulas from 4. Primal correction and adjoint solves should use \(\widehat A_h\) and \(\widehat A_h^\top\) directly; explicit formation of \(\widehat H_h\) is unnecessary except when the randomized sensitivity procedure is used. The full-column-rank overdetermined extension is given in 11.2.

5. Supply the transfer layer. Choose an admissible comparison vector \(U_h^\pi\in\mathcal{U}_h\) and rigorous bounds \(T_h\) and \(\varepsilon_{i,h}^{\mathrm{eval}}\) as in 4.3. For coercive conforming Galerkin problems, use the specialization in 5.2. In the fully computable one-dimensional reaction–diffusion specialization of 3, take \(U_h^\pi=U_h^\star=U_h^{\mathrm{FE}}\), so \(T_h=0\) and \(\varepsilon_{i,h}^{\mathrm{eval}}\) is given explicitly by 48 .

6. Report the interval. Use the corrected interval from 1 when exact correction is available, the randomized intervals from 5.4 when many queries are present, the localized form ?? when locality is exploited, or the norm-only radius 35 when only residual-norm information is desired. Certified inexact-solve allowances are given in 11.3. Add \(\zeta_{i,h}^{\theta}\) only when the reported center is an unprojected neural output rather than the certified reconstruction.

8 Experiments↩︎

The experiments are organized by claim rather than by benchmark. The objective is to evaluate certification mechanisms rather than prediction error in isolation. Manufactured tests provide direct validation for nodal and off-grid point queries and isolate the refinement behavior of point queries in three dimensions. The SimJEB experiments address whether the compatible-query algebra remains interpretable and computable in a large three-dimensional elasticity benchmark, both for compatible linear states and for raw neural predictions after projection into the admissible space. For that reason, this section reports the SimJEB query geometry, adjoint Green structure, corrected centers, released-field comparison identity, projected-neural raw-output intervals, and randomized sensitivity results; secondary SimJEB diagnostics such as raw-versus-corrected fields, localization sweeps, inexact correction, and additional load-family studies are given in 10. Details on the main claims tested by each experiment are gathered in ¿tbl:tab:experiment95map?.

All linear experiments in the square-system formulation use square reduced systems after admissible reduction. The full-column-rank overdetermined extension is illustrated in 10.2.

Unless stated otherwise, all reported linear intervals use the compatible residual \(r_h=A_hU_\theta-b_h\), the corresponding query vectors \(\ell_{i,h}\), and the transfer inputs from 4.3. In the SimJEB released-field experiment, the released nodal displacement field is not treated as a certified continuum reference solution. It is used only as an external admissible comparison vector \(U_h^\pi\) on the same mesh, so that the discrete comparison identity from 4 can be checked directly against a public benchmark field.

max width=

8.1 Necessity of compatibility↩︎

3 illustrates a failure mode excluded by the theory, and discussed in 1. A smooth bump is inserted between collocation points while the field remains unchanged at every sampled location. The sampled loss is therefore unchanged, but the hidden query value changes substantially. A compatible certificate reflects this change because the certified object is a mesh-based numerical field whose reported value and residual are defined on the same degrees of freedom. The behavior is consistent with 1: the sampled loss remains unchanged, while the hidden-point error increases and the compatible mesh and unprojected-output radii reflect the hidden error.

None

8.2 Linear calibration on a manufactured benchmark↩︎

The first controlled linear test evaluates the identities from 4 and the radius decomposition from 5.1. The benchmark includes both nodal and off-grid queries and is evaluated over a refinement sequence. 4 shows three effects. First, both the reconstruction-centered and corrected-center certificates contract under refinement and remain above the corresponding errors. Second, the corrected center agrees with the compatible discrete target to the resolution shown, as predicted by 5. Third, the median width decomposition separates the signed discrete correction, the transfer/comparison term, and the off-grid evaluation term instead of aggregating them into a non-decomposed interval width. The full-column-rank overdetermined extension, including the effective-residual refinement, is given in 10.2. Together, these results are consistent with calibration to the chosen numerical system rather than only empirical association with the observed error.

Certified early stopping for primal and adjoint solves for this benchmark is reported separately in 13.

None

8.3 Point queries in three dimensions↩︎

A manufactured three-dimensional refinement study was used to separate the point-query behavior of the signed transfer term from that of the norm route in 5.2. The reported quantities are interior nodal point queries and fixed-radius averages of the same scalar displacement component. For the point queries, the comparison state is chosen so that the evaluation term vanishes at the queried nodes, leaving the signed term \(\abs{g_{i,h}^{\top}\tau_h^\pi}\) and the norm route \(\sigma_{i,h}T_h\) as the relevant transfer quantities.

5 shows three effects. Panel A reports the maximum sensitivity over the query set. The point-query sensitivity grows under refinement, with observed slope \(-0.55\), while the fixed-radius averages remain nearly flat over the same meshes. Panel B compares the maximum pointwise bias \(E_{i,h}^{\mathrm{pt}}\) with the maximum signed term \(\abs{g_{i,h}^{\top}\tau_h^\pi}\); the largest discrepancy over the refinement sequence is \(7.16\times 10^{-14}\). Panel C compares the same pointwise bias with the maximum norm route \(\sigma_{i,h}T_h\). The norm route remains valid but wider, with ratio \(17.1\) on the finest mesh. This is the point-query behavior that motivates the use of fixed-radius averages in the large-scale three-dimensional study below.

8.4 Large-scale SimJEB validation: query geometry, Green sensitivities, correction, and transfer↩︎

The pointwise claims are tested on manufactured nodal and off-grid queries, including the refinement study in 8.3. The SimJEB study then addresses whether the compatible-query algebra remains query-specific and computable in a large three-dimensional benchmark when the reported quantities are patch-average displacement functionals assembled in the same admissible coefficient basis. Those quantities fall directly under the scalar compatible-query setting of 3 .

The benchmark is SimJEB, a public finite-element simulation dataset of crowdsourced jet-engine bracket designs; each entry includes a CAD file, tetrahedral mesh, triangular surface mesh, and structural simulation results for four load cases [33]. We use public bracket 148. A compatible \(P_1\) tetrahedral small-strain linear-elasticity system was reconstructed from the released mesh, material parameters and load specification, support conditions, and nodal displacement file. After support elimination, the reduced system is square. The vertical load case is used as the primary test problem. The four reported scalar queries are patch-average surface \(u_z\) values on a loaded interface patch, a near-support patch, a high-deflection bridge-arm patch, and a remote low-sensitivity patch. The benchmark geometry and the query-dependent adjoint Green maps are given directly in 6; the raw-versus-corrected displacement fields are provided in 17.

None

The Green maps quantify the querywise sensitivities used in the decomposition, localization, and randomized-sensitivity results. They indicate the physical structure of the querywise sensitivities \(\sigma_{i,h}\) used later in the analysis. The load-interface and high-deflection queries exhibit broader and stronger concentration, whereas the near-support and remote queries are more localized.

Starting from zero displacement, \(12\) CG iterations on the reduced compatible system produced the inexact state \(U_\theta\). The compatible residual \(r_h=K_hU_\theta-f_h\) was then formed using the same assembled stiffness matrix, eliminated support constraints, and reconstructed load vector. For each patch query, an adjoint solve gave the signed correction \(g_i^\top r_h\) and the corresponding norm-only certificate.

7 reports the patch-averaged discrete error, the signed correction, and the wider norm-only width for the four patch queries. The signed correction matches the actual discrete error at all four patches, as predicted by 3. The norm-only width is consistently conservative, and its conservatism is query-dependent in direct agreement with the heterogeneous Green concentration in 6. The near-support and remote patches have smaller norm widths, while the load-interface and high-deflection patches have larger sensitivities.

None

The transfer layer was then tested using the released SimJEB vertical-case nodal displacement field as an external admissible comparison state \(U_h^{\mathrm{rel}}\). With \[\tau_h^{\mathrm{rel}}:=K_hU_h^{\mathrm{rel}}-f_h,\] the comparison identity from 4 gives \[\ell_i^\top(U_h^\star-U_h^{\mathrm{rel}}) = -g_i^\top\tau_h^{\mathrm{rel}}.\] This does not certify the released field as a continuum reference solution. It tests the discrete comparison identity against an independently released admissible field on the same mesh. Because the figure reports absolute values, the sign convention does not affect the displayed comparison.

8 shows that the signed comparison identity holds to numerical precision for all four patch queries. The maximum identity gap is \(2.67\times10^{-12}\,\mathrm{mm}\). The worst released-field patch bias is \(8.11\times10^{-8}\,\mathrm{mm}\), attained at the high-deflection patch. The norm comparison bound remains conservative: the four norm bounds are \(2.62\times10^{-4}\), \(2.02\times10^{-5}\), \(3.69\times10^{-4}\), and \(1.73\times10^{-5}\,\mathrm{mm}\) for the load-interface, near-support, high-deflection, and remote patches, respectively. These results separate two effects. The signed comparison identity is sharp for the released comparison field, whereas the Cauchy–Schwarz norm bound is valid but substantially more conservative.

None

8.5 Projected neural certification on the SimJEB load family↩︎

The previous SimJEB study starts from an inexact compatible solve and evaluates the coefficient-space identities on one load case. The present test instead starts from a raw neural predictor and evaluates the raw-output interval 31 after explicit projection into the admissible space. The benchmark is again SimJEB bracket 148, now with the four-load family \[F_\mu = \mu_1 F_{\mathrm{vertical}} + \mu_2 F_{\mathrm{horizontal}} + \mu_3 F_{\mathrm{diagonal}} + \mu_4 F_{\mathrm{torsion}}.\] For each load, the neural model produces a raw coefficient vector \(\widetilde{U}_\theta(\mu)\), which is then projected to an admissible compatible state \(U_\theta(\mu)\in\mathcal{U}_h\). Certification uses the residual of this projected state, while the reported interval remains centered at the unprojected neural output. The test set contains \(256\) loads and \(28\) scalar queries: four patch-average surface queries, sixteen off-grid surface-point queries, and eight off-grid interior-point queries, for a total of \(7168\) query evaluations. For evaluation, a fresh compatible CG solve on the same reduced system provides the discrete target \(U_h^\star(\mu)\). This experiment therefore evaluates the discrete stage of the raw-output certificate, rather than the continuous transfer layer.

9 summarizes the resulting behavior. Panel A compares the raw neural error with the certified raw-output radius over all test query evaluations. The points lie on or above the line \(y=x\), and most are concentrated near that line. This is consistent with the corrected center matching the compatible discrete target to solver tolerance, so the numerical width in 30 is largely governed by \(\abs{u_\theta(x_i)-c_{i,h}^{\mathrm{corr}}}\). Panel B reports coverage \(1.0000\) for patch-average, surface-point, and interior-point queries. Panel C shows that the representation mismatch is small for most evaluations but not uniformly negligible; its upper tail is larger for interior queries than for the other query classes. Panel D compares the per-load runtime with a fresh compatible CG solve.

The aggregate values are reported in ¿tbl:tab:simjeb-projected-neural?. Over all \(7168\) test query evaluations, the raw-output interval contains the compatible discrete target in every case. The median raw error and median raw-output radius are \(3.6451\times 10^{-4}\) and \(3.6479\times 10^{-4}\), while the corresponding q95 values are \(3.9033\times 10^{-3}\) and \(3.9033\times 10^{-3}\). The q95 representation mismatch is \(1.1042\times 10^{-3}\). The corrected-center error remains at solver tolerance, with maximum \(3.5227\times 10^{-12}\). The median projected-neural certification runtime is \(8.61\times 10^{-2}\,\mathrm{s}\), compared with \(1.59\times 10^{1}\,\mathrm{s}\) for a fresh compatible CG solve, corresponding to a median speedup of \(184.3\times\).

8.6 Localization and the necessity of an explicit tail term↩︎

10 uses a dumbbell-domain benchmark to test the projector-based local/tail split from 7. The three query locations are chosen to induce qualitatively different adjoint geometries: near the source, in the bridge, and in the remote lobe. The local-only estimate is not uniformly reliable across these cases. Near the source, the local contribution dominates and the tail contribution decays rapidly. In the bridge and remote-lobe cases, a nontrivial tail remains necessary over a substantial range of radii. The certified local+tail estimate remains consistent with the actual point error across the sweep, whereas the local-only and tail-only terms are each insufficient when used alone. This is the computational implication of 7: locality is query-dependent and a deterministic certificate requires an explicit tail term.

None

The same mechanism appears on the large-scale SimJEB queries. The Green maps in 6 already show the corresponding query dependence, while the bound-term density maps and full SimJEB radius sweeps are retained in [fig:simjeb-bound-term-maps,fig:simjeb-localization].

8.7 Simultaneous randomized sensitivity estimation↩︎

The exact computation of all query sensitivities becomes increasingly expensive as the number of queries grows. 11 therefore evaluates the Gaussian probing construction from 5.4 on the SimJEB bracket. The simultaneous upper bound attains the nominal \(95\%\) all-query coverage across the tested probe counts. The inflation factor decreases rapidly when moving from four to sixteen probes and then plateaus, indicating that additional probes give limited further reduction. The representative trial shows the expected one-sided bound: the randomized quantities remain above the exact sensitivities. The runtime panel shows the corresponding cost of the probe solves. The randomized construction therefore reduces sensitivity estimation to a fixed number of probe solves while retaining a simultaneous coverage statement. The querywise spread of the randomized inflation factors is consistent with the heterogeneous Green concentration already visible in 6.

Additional synthetic probing diagnostics and direct runtime scaling against exact adjoints are reported in [fig:toy-randomized,fig:probe-scaling].

None

8.8 Nonlinear certification via verified exact roots↩︎

The nonlinear exact-root experiment tests 4. 12 shows that the corrected certificate contracts under refinement and remains substantially tighter than the norm-only version. The width decomposition further shows that, on this problem, the comparison-to-truth term dominates the nonlinear remainder term by several orders of magnitude. The remainder is nevertheless retained because it is required for verification in the certified root neighborhood. The verified-ball diagnostics show that the certified root neighborhood stays small relative to the base residual scale, and the per-query plot confirms that the reported radii remain comparable to the corrected errors query by query. The stationary-target variant exhibits the same qualitative behavior and is provided in 16.

None

Across these tests, the empirical results are consistent with the certificate chain. Compatibility is necessary. The linear corrected center recovers the compatible discrete target. The width decomposition separates discrete error, sensitivity, and transfer effects instead of aggregating them into a single scalar. On the SimJEB load family, the projected-neural study shows that the raw-output interval remains valid after projection into the admissible space, with a nonzero upper-tail mismatch term and per-load certification cost well below that of a fresh compatible solve. Localization requires an explicit tail. On SimJEB bracket 148, the released-field comparison experiment verifies the signed transfer identity against an external public displacement field and shows the expected conservatism of the norm comparison bound. Gaussian probing preserves simultaneous coverage at lower cost. The nonlinear exact-root construction remains non-vacuous after the remainder term is included. 10 reports corresponding results for inexact linear algebra, nonlinear stationary targets, supplementary SimJEB diagnostics, and additional load-family studies.

9 Conclusion↩︎

In this work, we introduce a constructive framework for certified error bounds of numerical PINNs. The analysis shows that pointwise residual minimization alone is insufficient to guarantee convergence over the full domain, motivating the need for compatibility between the learned solution, the discrete numerical system, and the query mechanism used for evaluation. This leads naturally to discretized domains with sampled query points and certified transfer operators to the continuous setting. Mesh-based \(P_1\) finite element spaces provide a particularly effective realization of this principle, although the proposed framework is not restricted to this choice.

The resulting certification decomposes the total error into explicit and interpretable contributions associated with the discrete residual, the discrete-to-continuous transfer, and the evaluation comparison. For linear problems, these quantities can be computed exactly and deterministically, while localized, randomized, and nonlinear extensions broaden the framework to more challenging settings. The theoretical results are supported by experiments on synthetic benchmarks and on SimJEB, a public dataset containing geometrically complex three-dimensional domains.

More broadly, this work establishes a bridge between neural PDE solvers and certified numerical analysis by combining learned approximation with rigorous verification. Rather than replacing classical numerical methods, the framework identifies how neural operators can be integrated into reliable scientific computing pipelines while retaining explicit control of approximation error. Promising directions for future work include combining the computed sensitivities with conformal prediction for sharper calibration, extending the framework to ill-posed and low-regularity systems, and coupling the certification layer with geometry-aware graph, hypergraph, and sheaf neural architectures. In particular, tangent-bundle and sheaf neural constructions suggest ways to generate vector-valued admissible fields on curved or heterogeneous geometries [34], while transferability results for hypergraph neural networks suggest a possible route toward certified amortized prediction across changing mesh, load, and interaction structures [35]. A further long-term direction is to express compatibility itself in a sheaf-theoretic local-to-global language, especially for multiphysics or heterogeneous-domain problems where local admissible data must agree across interfaces [36].

10 Additional Experiments and Reporting Details↩︎

This appendix gives experimental material that supports the certificate results and clarifies the numerical interpretation of the intervals.

10.1 Certification under inexact linear algebra↩︎

13 isolates the algebraic allowance terms from [prop:inexact_primal_correction,prop:inexact_adjoint_certificate]. On the manufactured benchmark, the certified primal and adjoint allowances decay in step with the actual corrected-center shift as the CG tolerance is tightened. The sensitivity inflation is close to one once the solves exit the low-accuracy regime, and the full-column-rank overdetermined example shows the same conservative behavior.

15 repeats the primal-correction study on the SimJEB bracket. The observed center shift decreases with the number of correction iterations, and the certified allowance remains above it at each patch query and for the max-over-queries statistic. These plots illustrate certified early stopping: the remaining linear-algebra error can be reported explicitly rather than absorbed without decomposition into the interval.

None

10.2 Full-column-rank overdetermined extension↩︎

This subsection presents the manufactured linear example that isolates the full-column-rank overdetermined extension from 11.2. This panel isolates the effective-residual refinement in the overdetermined setting.

None

None

10.3 Nonlinear least-squares stationary targets↩︎

16 evaluates the stationary-target certificate from 5. The corrected certificate contracts with refinement and is substantially tighter than the norm-only version. The stationarity remainder is small on this problem relative to the comparison-to-truth term, but it is retained because the theorem is stated around a verified stationary point rather than an exact discrete root. The verified stationary-ball radii remain small, and the per-query radii stay aligned with the corrected errors.

None

10.4 Supplementary SimJEB field and localization diagnostics↩︎

8.4 gives the bracket/query geometry and query-dependent Green maps for the SimJEB benchmark. This subsection gives the raw-versus-corrected field comparison, bound-term density maps, and full localization sweeps, which support interpretation but are not required for the primary theorem-level narrative.

17 compares the uncorrected early-stopped field with the corrected field that coincides with the compatible discrete target. The difference field is spatially structured rather than uniform, which is why the decomposition is reported patch by patch.

19 performs a radius sweep for the four patch queries. The certified local+tail estimate remains above the actual error across the sweep, whereas the local-only and tail-only terms are each inadequate on parts of the radius range. The load-interface and high-deflection queries retain a sizeable tail over long radii, while the near-support and remote queries localize more rapidly.

None

None

None

10.5 Supplementary diagnostics for randomized sensitivity estimation↩︎

The synthetic probing tests complement the SimJEB results. 20 shows the same two qualitative trends on a controlled benchmark: the simultaneous upper bound attains the nominal coverage and the inflation factor decays with the number of probes. 21 isolates the computational tradeoff. Exact adjoints scale with the number of queries, whereas probing cost is essentially fixed by the number of probe solves and becomes advantageous for large query sets. The second panel shows the corresponding cost–tightness tradeoff: more probes cost more but reduce the inflation of the \(95\)th-percentile bound.

None

None

10.6 Amortized load families↩︎

The load-family experiments address the many-instance setting discussed in 3. 22 studies a linear reduced load-family surrogate. Increasing the output rank reduces both the certified test quantities and the residual energy, with the signed correction remaining close to the actual 95th-percentile maximum query error and the norm-only width remaining wider.

23 reports the corresponding physics-trained neural load-family model. The training and validation losses stabilize, the certified 95th-percentile maximum signed correction tracks the actual 95th-percentile maximum query error, and rank three gives the smallest reported certified test quantity among the tested neural models while substantially reducing residual energy relative to rank two. These results characterize certificate behavior in the amortized regime, where the neural model produces many compatible coefficient vectors and the certificate verifies each one.

None

None

10.7 Experimental configuration summary↩︎

¿tbl:tab:app-experimental-configuration? gathers the configuration used in the numerical experiments throughout the paper.

For the main linear manufactured benchmark, the data are generated on \(\Omega=(0,1)^2\) with homogeneous Dirichlet boundary conditions from \[-\Delta u^\star=f,\] using the manufactured solution \[\begin{align} u^\star(x,y) &= x(1-x)y(1-y) \exp\!\left( -\frac{(x-0.7)^2+(y-0.35)^2}{0.12^2} \right) \\ &\quad +0.1\sin(3\pi x)\sin(2\pi y). \end{align}\] The forcing is generated by applying the operator to this solution, \(f=-\Delta u^\star\), so \(u^\star\) is the exact continuous solution by construction. The compatible \(P_1\) systems are assembled on structured meshes with \(8\), \(16\), \(24\), and \(32\) cells per coordinate direction. For each mesh, we use both the finite-element solution \(u_h^{\mathrm{FE}}\) and the nodal interpolant \(u_h^\pi=\mathcal{J}_hu^\star\). The trained boundary-vanishing coordinate MLP is evaluated through its nodal values in the admissible \(P_1\) space, and the fixed query set contains \(12\) interior nodal queries and \(12\) off-grid point queries.

max width=

11 Auxiliary Analytical Material↩︎

11.1 Continuous residual-to-point stability route↩︎

This subsection gives an alternative route that is not used in the certificate chain. Instead of certifying a compatible discrete target and then transferring to the continuous solution, one can work directly with a continuous residual norm strong enough to control point evaluation.

Assume in this subsection that \(\mathcal{A}:V\to Y\) is linear and bounded, where \[Y:=Y_\Omega\times Y_{\partial\Omega}\] is a Hilbert residual space with norm \(\norm{\cdot}_Y\) and inner product \(\langle\cdot,\cdot\rangle_Y\). Define \[\mathcal{A}u:=(\mathcal{L}u,\mathcal{B}u)\in Y, \qquad b:=(f,g)\in Y.\] Also define the continuous residual map \[\mathcal{F}(u):=\mathcal{A}u-b.\] For a fixed query point \(x\in\overline{\Omega}\), let \[\ell_x(u):=u(x).\] If point evaluation is controlled directly by the continuous residual, define \[S_x := \sup_{\substack{w\in V\\ \mathcal{A}w\neq 0}} \frac{\abs{\ell_x(w)}}{\norm{\mathcal{A}w}_Y}.\]

If, in addition, there exists a residual-side Green representative \(G_x\in Y\) such that \[\ell_x(w)=\langle G_x,\mathcal{A}w\rangle_Y \qquad \text{for all }w\in V,\] then \(S_x\le \norm{G_x}_Y\), and the pointwise error admits the exact representation \[u(x)-u^\star(x) = \langle G_x,\mathcal{F}(u)\rangle_Y.\] For any orthogonal projector \(\Pi_{x,R}:Y\to Y\), this yields the local/tail split \[\begin{multlined} \abs{u(x)-u^\star(x)} \le \abs{\langle \Pi_{x,R}G_x,\mathcal{F}(u)\rangle_Y}\\ + \norm{(I-\Pi_{x,R})G_x}_Y\, \norm{(I-\Pi_{x,R})\mathcal{F}(u)}_Y. \end{multlined} \label{eq:dual95green95local95tail95bound}\tag{66}\] This is the direct continuous analogue of the discrete local/tail split used in the primary formulation.

11.2 Full-column-rank overdetermined linear extension↩︎

This subsection gives the full-column-rank overdetermined linear extension. The compatible discrete target is now defined by least squares after admissible reduction, while rank-deficient reduced operators remain outside the certified scope.

Define the compatible discrete target by \[y_h^\star := \operatorname*{arg\,min}_{y\in\mathbb{R}^{n_h}} \frac{1}{2}\norm{\widehat A_hy-\widehat b_h}_{W_h}^2, \qquad U_h^\star:=U_{0,h}+Z_hy_h^\star. \label{eq:rectangular95discrete95target}\tag{67}\] The associated normal equations are \[\widehat A_h^\top W_h(\widehat A_hy_h^\star-\widehat b_h)=0. \label{eq:rectangular95normal95equations}\tag{68}\]

For inexact linear algebra in this full-column-rank overdetermined setting, the least-squares definition of the compatible discrete target is intrinsic, but correction and adjoint quantities should still be computed matrix-free by weighted QR, LSQR/LSMR, or preconditioned Krylov methods rather than by explicitly forming \(\widehat H_h^{-1}\). The certified allowance formulas in 11.3 can then be applied directly.

11.3 Auxiliary results for linear certification↩︎

The overdetermined projector identity from 11.2 also implies the following coefficient-norm consequence: \[\norm{U_\theta-U_h^\star}_2 \le \norm{Z_h}_2\, \lambda_{\min}(\widehat H_h)^{-1/2} \norm{r_h}_{W_h}. \label{eq:euclidean95coefficient95error95bound}\tag{69}\]

11.3.1 Inexact linear algebra for certification↩︎

Inexact linear algebra should distinguish the square reduced square-system setting from the full-column-rank overdetermined extension. If \(\widehat A_h\) is square and invertible, then the square-system formulas show that the exact primal correction and adjoint sensitivities are obtained from \[\widehat A_h e_\theta = r_h, \qquad \widehat A_h^\top z_{i,h} = \widehat\ell_{i,h},\] with \(g_{i,h}=z_{i,h}\). In that regime the computation should not be performed through \(\widehat H_h=\widehat A_h^\top W_h\widehat A_h\), because that squares the condition number without changing the exact corrected center.

If \(\widehat A_h\) is full-column-rank overdetermined, the normal matrix is intrinsic to the compatible discrete target, but it should still be accessed matrix-free through weighted QR, LSQR/LSMR, preconditioned Krylov methods, or the residual-side constrained problem of 13, rather than by explicitly forming \(\widehat H_h^{-1}\).

Any upper bound \(\bar\sigma_{i,h}\ge\sigma_{i,h}\) may replace \(\sigma_{i,h}\) in ?? ; this is useful when \(\sigma_{i,h}\) itself is obtained by an inexact or randomized procedure.

A sufficient computable choice for \(\delta_h^{\mathrm{prim}}\) or \(\delta_{i,h}^{\mathrm{adj}}\) follows from any coercivity bound \[\widehat H_h \succeq \alpha_h P_h,\] where \(P_h\) is symmetric positive definite and \(\alpha_h>0\). Since \[\widehat H_h^{-1}\preceq \alpha_h^{-1}P_h^{-1},\] one has \[\left( (\eta_h^{\mathrm{prim}})^\top\widehat H_h^{-1}\eta_h^{\mathrm{prim}} \right)^{1/2} \le \alpha_h^{-1/2}\norm{\eta_h^{\mathrm{prim}}}_{P_h^{-1}},\] and \[\left( s_{i,h}^\top\widehat H_h^{-1}s_{i,h} \right)^{1/2} \le \alpha_h^{-1/2}\norm{s_{i,h}}_{P_h^{-1}}.\] Thus the right-hand sides are valid certified choices for \(\delta_h^{\mathrm{prim}}\) and \(\delta_{i,h}^{\mathrm{adj}}\). In particular, with \(P_h=I\), any verified lower eigenvalue bound \(\lambda_{\min}(\widehat H_h)\ge \alpha_h>0\) is sufficient.

11.3.2 Verified exact-root ball for square systems↩︎

12 Nonlinear Stationary-Target Variant↩︎

This appendix isolates the nonlinear least-squares stationary-target construction so that the exact-root chain remains focused.

12.1 Certification of nonlinear least-squares stationary targets↩︎

For this section, assume in addition that \(F_h\) is twice continuously differentiable on the relevant neighborhood. Define \[\Phi_h(y):=\frac{1}{2}\norm{F_h(y)}_{W_h}^2, \qquad S_h(y):=DF_h(y)^\top W_hF_h(y)\in\mathbb{R}^{n_h}. \label{eq:nonlinear95stationarity95map}\tag{71}\] A nonlinear least-squares target \(y_{h,\mathrm{ls}}^\star\) satisfies \[S_h(y_{h,\mathrm{ls}}^\star)=0.\]

For any base point \(y\), define \[K_h(y):=DS_h(y).\] Explicitly, for \(v\in\mathbb{R}^{n_h}\), \[K_h(y)v = J_h(y)^\top W_hJ_h(y)v + \bigl[D^2F_h(y)[v]\bigr]^\top W_hF_h(y). \label{eq:nonlinear95ls95hessian}\tag{72}\] If \(K_h(y)\) is invertible, define \(q_{\ell,h}^{\mathrm{ls}}(y)\) by \[K_h(y)^\top q_{\ell,h}^{\mathrm{ls}}(y)=\widehat\ell_h.\] For a deviation \(d\in\mathbb{R}^{n_h}\), define \[R_h^{\mathrm{ls}}(y;d) := S_h(y-d)-S_h(y)+K_h(y)d. \label{eq:ls95base95point95remainder}\tag{73}\] If \(DS_h\) is Lipschitz on \(B_{X_h}(y,t)\) with constant \(L_h^{\mathrm{ls}}(y,t)\) in the operator norm from \((\mathbb{R}^{n_h},\norm{\cdot}_{X_h})\) to \((\mathbb{R}^{n_h},\norm{\cdot}_{X_h^\ast})\), then \[\norm{R_h^{\mathrm{ls}}(y;d)}_{X_h^\ast} \le \frac{L_h^{\mathrm{ls}}(y,t)}{2}\norm{d}_{X_h}^2, \qquad \norm{d}_{X_h}\le t,\] because \[R_h^{\mathrm{ls}}(y;d) = \int_0^1 \bigl(K_h(y)-K_h(y-sd)\bigr)d\,\,\mathrm{d}s.\]

A radius \(t_h^{\mathrm{ls}}(y_\theta)\) can be certified by applying a square Newton argument to the stationarity map \(S_h\). If \(K_h(y)\) remains positive definite on the verified ball, the stationary point is a strict local minimizer of \(\Phi_h\) in that ball.

13 Proofs↩︎

13.1 Proofs for 3 and 11↩︎

13.1.1 Proof of 1↩︎

13.1.2 Proof of 10↩︎

13.1.3 Proof of 16↩︎

13.1.4 Proof of 13↩︎

13.2 Proofs for 4 and 11.2↩︎

13.2.1 Proof of 11↩︎

13.2.2 Proof of 12↩︎

13.2.3 Proof of 2↩︎

13.2.4 Proof of 3↩︎

13.2.5 Proof of 4↩︎

13.2.6 Proof of 1↩︎

13.2.7 Proof of 5↩︎

13.3 Proofs for 5 and 11↩︎

13.3.1 Proof of 14↩︎

13.3.2 Proof of 6↩︎

13.3.3 Proof of 2↩︎

13.3.4 Proof of 1↩︎

13.3.5 Proof of 3↩︎

13.3.6 Proof of 2↩︎

13.3.7 Proof of 7↩︎

13.3.8 Proof of 15↩︎

13.3.9 Proof of 17↩︎

13.3.10 Proof of 18↩︎

13.3.11 Proof of 19↩︎

13.4 Proofs for 6↩︎

13.4.1 Proof of 8↩︎

13.4.2 Proof of 9↩︎

13.4.3 Proof of 4↩︎

13.5 Proofs for 12↩︎

13.5.1 Proof of 20↩︎

13.5.2 Proof of 21↩︎

13.5.3 Proof of 5↩︎

References↩︎