The certification of machine learning methods has emerged as a critical challenge in scientific computing, particularly for physics-informed neural networks (PINNs) where reliability guarantees are essential for deployment in safety-critical applications. Significant progress has been made in deriving a posteriori upper bounds for PINN errors [1]. While upper bounds estimate the proximity of the exact and approximate solutions from above, lower bounds provide complementary information on the minimal separation between them. Complete certification therefore requires both upper and lower error enclosures to provide tight bounds on prediction quality. However, the complementary problem of establishing rigorous lower bounds remains largely unexplored.

Literature review. Neural network methods for solving ordinary and partial differential equations (ODEs and PDEs), especially PINNs, have become very popular in recent years, as discussed in [2]. PINNs incorporate physical constraints directly into the learning objective through differential operators, enabling them to leverage neural network expressivity while respecting physical consistency. This hybrid approach facilitates accurate predictions even with limited or noisy data by constructionally enforcing physical feasibility [3].

Certified error estimation for PINNs builds upon foundational developments in scientific machine learning [3, 4]. Existing approaches cover statistical uncertainty quantification and a deep confidence framework [5], dual-weighted residual methods [6], and goal-oriented error estimation [7], though these typically provide either statistical rather than rigorous bounds or require substantial training data.

Statistical error models quantify uncertainty only in a statistical sense: they rely on sufficiently rich and representative training data to construct confidence intervals for the prediction error and therefore do not, in general, provide deterministic rigorous certificates for a given trained PINN. In contrast, certified a posteriori error estimation aims at computable worst-case bounds for the prediction error of a specific network realization, without requiring knowledge of the exact solution and, in the sense of [1], even for inputs not used during training. In physics-informed problems, we often lack a large labeled dataset of true errors, so direct application of deep confidence is limited. Drawing from certification frameworks in reduced-order modeling [8] and classical stability theory for ODEs [9], this work establishes rigorous a posteriori error bounds that remain valid for unseen data without requiring a priori knowledge of the true solution. Theoretical analyses by [10] further advance PINN reliability by deriving rigorous, dimension-robust a priori approximation error bounds for both solution and operator-learning networks, while providing constructive and efficient guidelines that explicitly relate network architecture to achievable accuracy.

While the certification of PINNs for ODEs has received increasing attention, extending such guarantees to PDEs poses substantially greater theoretical and computational challenges. In the PDE setting, certification requires accounting for spatial derivatives, boundary and initial conditions, and stability properties of the underlying operator. Recent studies have begun to address these issues by deriving rigorous a posteriori error bounds and verification frameworks for PDE-defined PINNs; see [11, 12] and references therein. Hillebrecht and Unger [12, 13] established one of the first rigorous certification schemes for PINNs governed by classical PDEs such as the heat, transport, and Navier–Stokes equations. Their approach derives provable upper bounds on the prediction error using PDE stability estimates and residual norms, enabling verification without access to the exact solution. In parallel, [14] introduced an efficient certification method based on bounding the global residual, providing computable correctness conditions analogous to solver tolerances. Furthermore, [15] proposed a functional-type a posteriori framework that integrates classical majorant estimators into PINN training, yielding guaranteed error upper bounds with high computational efficiency. Guo and Haghighat [16] extended these ideas to elasticity problems, employing the constitutive relation error to obtain both global and goal-oriented error estimates. However, establishing corresponding a posteriori lower bounds, which would provide guaranteed measures of the minimal attainable error, remains largely unexplored in the current literature.

Practical applicability and scope. The theory developed below is designed for ODE models and operating regimes in which the vector field satisfies suitable local one-sided growth conditions on a certified state-space domain. Representative examples include early-growth epidemic models and their multigroup linearizations [17, 18, 19], low-incidence SIS-type reductions, and autocatalytic or other positive-feedback chemical kinetics in batch and flow reactors [20, 21, 22]. From a practical viewpoint, the relevant issue is therefore local rather than global: one verifies the required sign conditions for the symmetric part of the Jacobian on a certified domain containing both the exact and approximate trajectories. This motivates the localized framework adopted in the present paper.

Initial conditions and informativeness of lower bounds. The informativeness of a lower certificate depends crucially on how initial conditions are imposed in the PINN. In standard PINNs, the initial condition is typically enforced softly through a data-misfit term in the loss, rather than exactly by construction [3]. Hard-constraint formulations instead embed the condition into the network architecture or trial ansatz so that it is satisfied exactly by construction [23, 24]. In the latter case, the scalar lower bound can become trivial after nonnegativity post-processing because the anchoring error vanishes. We therefore do not claim that the scalar estimator alone can certify a positive lower error in this regime. Instead, when the model is linear, we use the signed residual itself through a finite set of coordinate probes which gives a full-interval lower diagnostic. This distinction is practically relevant because recent work on stiff time-dependent PINNs shows that exact hard enforcement of initial and boundary conditions can materially affect robustness and training efficiency [24].

Organization and contributions. Accordingly, the main contributions of this manuscript are:

• computable localized a posteriori lower bounds for PINN errors on certified state-space domains under a strong monotonicity condition,
• complementary localized a posteriori upper bounds under a one-sided Lipschitz condition, extending earlier certification results based on Lipschitz assumptions and yielding potentially sharper upper certificates on certified domains [1, 12],
• explicit formulas for linear time-invariant and time-varying systems in terms of the extremal eigenvalues of the symmetric part of the system matrix,
• a certified-containment criterion that makes precise when the local assumptions remain valid and when certification must terminate,
• a signed-residual finite-probe lower certificate for linear inhomogeneous systems, implemented with coordinate unit vectors and remaining informative when the scalar lower bound vanishes,
• a certificate-informed training formulation based on the propagated upper certificate and its auxiliary scalar ODE representation.

These two-sided certified error bands provide mathematically rigorous enclosures of the PINN approximation error. A narrow band indicates tight certification, whereas a wide band reflects either conservative local constants or limited physical fidelity of the trained network. The band width therefore serves as a quantitative indicator of certificate sharpness and practical reliability.

The paper is organized as follows. Section 2 introduces the problem formulation, the localized certification setting on state-space domains, and the flow-map notation used throughout the paper. Section 3 presents the main a posteriori error estimates, including localized lower and upper bounds for nonlinear systems, sharpened formulas for linear time-invariant and time-varying systems, and the RK4-based numerical post-processing used to evaluate the certification terms. Section 4 explains the connection with physics-informed neural-network training. Section 5 develops the certificate-informed training formulation. Section 6 reports numerical examples illustrating the theory. Section 7 summarizes the main conclusions and limitations.

Notation. We conclude this introductory section by introducing standard mathematical notation used throughout the paper. The time interval is denoted by \( \mathbb{T}=[0,T]\) with \( T>0\) . For vectors in \( \mathbb{R}^n\) , \( \|\cdot\|\) denotes the Euclidean norm and \( \cdot\) the standard inner product. For a matrix \( A\) , \( A^\top\) denotes its transpose. For a finite set \( A\) , we write \( |A|\) for its cardinality. Unless stated otherwise, all norms in the theoretical developments are understood in the Euclidean sense, although several arguments extend to other equivalent norms on finite-dimensional spaces.

We consider the initial value problem on the time interval \(\mathbb{T}=[0,T]\):

where \(f:\mathbb{T}\times\mathbb{R}^n\to\mathbb{R}^n\) is continuous and \(x_0\in\mathbb{R}^n\).

A standard sufficient condition for local existence and uniqueness of solutions to (1) is continuity in \(t\) together with local Lipschitz continuity of \(f\) with respect to \(x\); see, for example, the Picard–Lindel
specialChar{34}of theorem [25]. In the present paper, however, we work on a certified convex state-space domain \(D\subset\mathbb{R}^n\). All structural assumptions are therefore imposed only on \(\mathbb{T}\times D\), and the resulting certificates are valid on every time interval on which both the exact trajectory and the PINN approximation remain in \(D\).

For clarity, we first formulate the assumptions with constant certified coefficients on \(D\). Thus, for all \(u,v\in D\) and all \(t\in\mathbb{T}\), we assume

Here,

• \(\ell_D\in\mathbb{R}\) is a local lower one-sided Lipschitz coefficient,
• \(\mu_D\in\mathbb{R}\) is a local upper one-sided Lipschitz constant,
• \(L_D>0\) is a local Lipschitz constant.

The terminology is important: the inequality (2) is a lower one-sided Lipschitz estimate in general, and it becomes a strong monotonicity condition precisely in the case \(\ell_D>0\). This is the regime in which the lower a posteriori error bound is informative. By contrast, (3) is the relevant assumption for upper error bounds; if \(\mu_D<0\), then the dynamics are contractive on \(D\).

The role of the three coefficients is therefore different. The local Lipschitz bound (4), together with continuity in \(t\), provides the standard well-posedness assumption on \(D\). The lower one-sided bound (2) is used for lower error certificates, whereas the upper one-sided bound (3) is used for upper error certificates. Since (4) implies (3) with the same constant, the upper one-sided Lipschitz assumption is weaker than the Lipschitz assumption used in earlier upper-bound results such as [1, 12].

If \(f\) is continuously differentiable with respect to \(x\), then these certified coefficients can be chosen from the symmetric part of the Jacobian. Writing

one may define

In the scalar case \(n=1\), this reduces to

For each initial value \(\xi\in D\) such that the corresponding solution exists on \(\mathbb{T}\) and remains in \(D\), we denote that solution by \(t\mapsto\varphi(t,\xi)\). This defines the associated local flow map

on the subset of \(\mathbb{T}\times D\) where the trajectory stays in \(D\). In particular, for the prescribed initial value \(x_0\), the exact trajectory is

In general, the flow map is not available in closed form and must therefore be approximated numerically. In this work, we use a physics-informed neural network. More precisely, we introduce a parameterized candidate map

where \(\theta\in\mathbb{R}^k\) denotes the trainable parameter vector of the network. The architecture determines the approximation class, while \(\theta\) collects the weights and biases. Throughout the paper, we assume that \(\hat{\varphi}_\theta\) is sufficiently smooth with respect to \(t\) so that the differential-equation residual is well defined.

For the prescribed initial value \(x_0\), the corresponding PINN trajectory is the slice

The training procedure returns a parameter vector, denoted by \(\theta^\star\), and we write

Here, \(\theta^\star\) simply denotes the trained parameter vector produced by the chosen optimization procedure.

The approximation error is then

Our goal is to estimate this error a posteriori, using only computable quantities derived from the trained PINN. More precisely, the certificates below depend on

• the trained trajectory \(\hat{x}\),
• a computable majorant of the residual,
• certified coefficients \(\ell_D\), \(\mu_D\), and \(L_D\) on the domain \(D\),
• and the initial mismatch \(\|e(0)\|=\|x_0-\hat{x}(0)\|\).

In soft-constrained PINNs, this initial mismatch is controlled by the initial-condition part of the training loss and is directly computable from the trained network. In hard-constrained PINNs, one has \(\hat{x}(0)=x_0\) by construction and therefore \(e(0)=0\).

Main Problem. Given a trained PINN approximation \(\hat{\varphi}_{\theta^\star}\), or equivalently the associated trajectory \(\hat{x}(t)=\hat{\varphi}_{\theta^\star}(t,x_0)\), derive computable and rigorous a posteriori certificates for the error (13) on \(\mathbb{T}\), in particular lower and upper bounds that can be evaluated at low computational cost without using the exact solution.

In this section, we derive computable a posteriori error bounds for a trained PINN approximation. We treat the problem in two stages. First, we establish localized lower and upper error bounds for nonlinear systems on a certified state-space domain under one-sided growth assumptions. Second, we specialize these results to linear systems, where the relevant constants admit explicit spectral representations in terms of the symmetric part of the system matrix. Finally, we record the RK4-based numerical post-processing used to evaluate the certification integrals.

3.1 Nonlinear case: localized lower and upper bounds

3.2 Nontrivial lower information beyond scalar residual majorants

3.3 Linear case: sharpened bounds

3.4 Auxiliary numerical certification by the classical fourth-order Runge–Kutta method

PINNs incorporate the governing differential equation into the training objective through residual minimization [3, 4]. In the present setting, the aim is to approximate the local flow map

introduced in Section 2 by a neural network

where \(\theta\in\mathbb{R}^k\) denotes the trainable parameter vector.

Since exact solution data are typically unavailable, training is based on collocation in time and initial-value space. Let

be a finite set of collocation points \(y=(t,\xi)\). The physics loss is then defined by

This loss penalizes violations of the ODE and is the training analogue of the residual quantity appearing in the a posteriori certification results.

For initial-value problems, two enforcement strategies are common.

Soft-constrained PINNs. In the soft-constrained setting, the initial condition is imposed through an additional penalty term. Given training initial values

one defines

The total loss takes the form

where \(\gamma_{\mathrm{physics}},\gamma_{\mathrm{initial}}>0\) are weighting parameters.

Hard-constrained PINNs. In the hard-constrained setting, the network ansatz is chosen so that the initial condition is satisfied exactly by construction [23]. In that case,

and no separate initial-condition penalty is required. The training objective then reduces to the physics loss, possibly combined with additional regularization terms.

Connection with the certified bounds. After training, let \(\theta^*\) denote the learned parameter vector, and fix an initial value \(x_0\in D_{\mathrm{adm}}\). The corresponding PINN trajectory is

with residual

The certification results of Section 3 apply to this trajectory once a computable residual majorant \(\delta\) is available such that

In the soft-constrained case, the anchoring error

is directly computable from the trained network and enters both the lower and upper certificates. In the hard-constrained case, one has

As a consequence, the scalar lower certificate from Section 3 may become trivial after nonnegativity truncation. In the linear setting, this loss of lower information can be repaired by the signed-residual finite-probe certificate in Theorem 2.

Thus, the proposed certification framework is compatible with both common PINN training paradigms. Soft-constrained PINNs naturally provide a computable initial mismatch for the scalar two-sided formulas. Hard-constrained PINNs may require directional residual information, such as the finite-probe certificate in the linear case, to obtain nontrivial lower diagnostics.

The certificates above are primarily a posteriori evaluation tools. Nevertheless, the upper certificate can also be used constructively during training. This use must be separated from lower-bound evaluation. The upper estimator is monotone in both the initial mismatch and the residual majorant, whereas the scalar lower estimator can be driven to the trivial value zero or can reward a larger initial mismatch if it is maximized. Thus the lower certificate is kept as a diagnostic, while the upper certificate can be used as an auxiliary regularizer.

Let \(\hat x_\theta\) be a PINN approximation with residual

and let \(\delta_\theta(t)\ge \|R_\theta(t)\|\) be a smooth computable residual majorant. For a constant \(c\) equal either to an upper one-sided Lipschitz constant \(\mu_D\) or to a Lipschitz constant \(L_D\), define

The auxiliary ODE formulation gives the training loss

where, on a sorted auxiliary grid \(Y=\{t_i\}_{i=0}^N\), one may choose

The auxiliary state is recomputed from the current network at each training step. It is not an additional neural-network output.

6.1 Nonlinear radial-growth plus rotation system

We first consider the nonlinear two-dimensional system

\[ \begin{equation} \dot x(t)=f(x(t)),~~ x(0)=x_0, \end{equation} \]

(179)

where

\[ \begin{equation} f(x)=\bigl(1+\alpha\|x\|_2^2\bigr)x+\beta Jx,~~ J=\begin{pmatrix}0&-1\\ 1&0\end{pmatrix},~~ x_0=(0.2,0)^\top, \end{equation} \]

(180)

with

\[ \begin{equation} \alpha=0.4,~~ \beta=6,~~ T=1. \end{equation} \]

(181)

This example is useful because the rotational term changes the Euclidean Lipschitz constant but cancels from the symmetric part of the Jacobian. In polar coordinates,

\[ \begin{equation} r'(t)=\bigl(1+\alpha r(t)^2\bigr)r(t),~~ \theta'(t)=\beta, \end{equation} \]

(182)

and hence

\[ \begin{equation} r(t)=\bigl((r_0^{-2}+\alpha)e^{-2t}-\alpha\bigr)^{-1/2},~~ x(t)=r(t)\begin{pmatrix}\cos(\beta t)\\ \sin(\beta t)\end{pmatrix}. \end{equation} \]

(183)

The Jacobian is

\[ \begin{equation} Df(x)=\bigl(1+\alpha\|x\|_2^2\bigr)I+2\alpha xx^\top+\beta J, \end{equation} \]

(184)

and therefore

\[ \begin{equation} \operatorname{sym}(Df(x))=\bigl(1+\alpha\|x\|_2^2\bigr)I+2\alpha xx^\top. \end{equation} \]

(185)

On the ball \(B_R=\{x:\|x\|_2\le R\}\), the constant certified coefficients are

\[ \begin{equation} \ell_{B_R}=1,~~ \mu_{B_R}=1+3\alpha R^2,~~ L_{B_R}\le 1+3\alpha R^2+|\beta|. \end{equation} \]

(186)

The skew-symmetric part \(\beta J\) is responsible for the large Lipschitz constant but does not affect \(\ell_{B_R}\) and \(\mu_{B_R}\).

Pan/Zoom

Pan/Zoom

Figure 1 shows that both trajectories remain in the certified ball. At the final time, the true error is \(3.86\times10^{-2}\), the global-ball lower certificate is \(3.30\times10^{-2}\), the local-segment diagnostic lower curve is \(3.47\times10^{-2}\), the global one-sided upper certificate is \(6.01\times10^{-2}\), and the local-segment diagnostic upper curve is \(4.71\times10^{-2}\). The corresponding global Lipschitz upper certificate is valid but much larger, with final value about \(2.27\times10^{1}\), because it carries the rotational contribution \(|\beta|\). Figure 2 confirms that the segment and secant curves are useful for understanding sharpness, while the constant ball certificate remains the main computable certificate.

Remark 14 (Constant examples and time-dependent extensions)

The reported certified constants in the main examples are constant on the certified domain. The local-segment and secant curves in Figure 1, Figure 2 are diagnostic comparisons for this synthetic problem. The general extension to time-dependent coefficients follows the integrating-factor remark in Section 2.

6.2 High-dimensional stiff linear ODE

The second example is a high-dimensional linear system, and it is also the place where the different linear constants are made explicit. We consider the stiff diagonal system in dimension \(d=96\):

\[ \begin{equation} x_i'(t)=-\lambda_i x_i(t), ~~ \lambda_i=10^{4(i-1)/(d-1)}, ~~ i=1,\ldots,d, \end{equation} \]

(187)

with \(x_i(0)=(-1)^{i-1}/\sqrt d\) and \(T=1\). The stiffness ratio is \(10^4\). A full flow-map approximation on a box in \(\mathbb R^{96}\) would suffer from the curse of dimensionality, whereas the present trajectory-wise certificate only requires the trained path, residual evaluations, and matrix-vector operations.

For \(A=-\operatorname{diag}(\lambda_1,…,\lambda_d)\),

\[ \begin{equation} m(A)=-10^4,~~ M(A)=-1,~~ \|A\|_2=10^4. \end{equation} \]

(188)

These values illustrate why positivity should not be imposed on one-sided constants. The system is dissipative, so the sharp one-sided constants are negative; this is valid for upper certification, but it means that no positive strong-monotonicity lower constant is available. To display all admissibility classes in one example, we also report non-sharp comparison constants

\[ \begin{equation} \ell_{\rm adm}=-1.2\times10^4\le m(A),~ \mu_{\rm adm}=0\ge M(A),~ L_{\rm adm}=1.2\times10^4\ge\|A\|_2. \end{equation} \]

(189)

Thus the six displayed numbers are two lower one-sided coefficients, two upper one-sided coefficients, and two Lipschitz constants; only \(m(A)\), \(M(A)\), and \(\|A\|_2\) are sharp.

The scalar lower certificate is not useful here because the system is dissipative and the hard-constrained ansatz has \(e(0)=0\). The signed-residual finite-probe certificate from Theorem 2 remains nontrivial. We use only the coordinate unit vectors \(P=\{e_1,…,e_d\}\); hence the reported lower curve is \(\max_i |e_i(t)|\) reconstructed from the signed residual identity, not the full norm \(\|e(t)\|_2\).

Pan/Zoom

Table 1. Six constants for the stiff linear system. The one-sided constants may be negative; the Lipschitz constants are nonnegative.

Class	Constant	Value
One-sided (lower)	\( m(A)\)	\( -10^{4}\)
One-sided (lower, admissible)	\( \ell_{\rm {adm}}\)	\( -1.2\times 10^{4}\)
One-sided (upper)	\( M(A)\)	\( -1\)
One-sided (upper, admissible)	\( \mu_{\rm {adm}}\)	\( 0\)
Lipschitz	\( \|A\|_{2}\)	\( 10^{4}\)
Lipschitz (admissible)	\( L_{\rm {adm}}\)	\( 1.2\times 10^{4}\)

Pan/Zoom

Figure 3, Table 1, Figure 4 show the effect of stiffness on certification. The maximum true error over the grid is \(4.01\times10^{-7}\). At the final time, the true error is \(1.19\times10^{-8}\), while the finite-probe lower certificate is \(7.94\times10^{-9}\), about \(66.5\%\) of the true error. The one-sided upper certificate is \(6.51\times10^{-5}\), whereas the Lipschitz upper certificate overflows the useful plotting scale. This example is the main use of the signed-residual finite-probe lower repair in high dimension.

6.3 Certificate-informed training on an unstable oscillatory system

The third example uses the certificate-informed training formulation from Section 5. We consider

\[ \begin{equation} x'(t)=Ax(t)+g(t),~~ A=\begin{pmatrix}0.35&-20\\ 20&0.35\end{pmatrix},~~ x(0)=0, \end{equation} \]

(190)

with a smooth forcing \(g(t)\). The sharp upper one-sided constant and the Euclidean Lipschitz constant are

\[ \begin{equation} M(A)=0.35,~~ \|A\|_2=20.003. \end{equation} \]

(191)

This is an unstable oscillatory system: the symmetric part gives the true amplitude growth, whereas the large skew-symmetric part makes the Lipschitz upper constant much larger. This is therefore a more informative test for upper-certificate training than a scalar ODE.

The PINN is hard-constrained at the initial condition. Stage one uses the standard residual loss. Stage two continues from the same checkpoint with three branches: ordinary PINN continuation, one-sided upper-certificate training based on \(U_{M(A)}\), and Lipschitz upper-certificate training based on \(U_{\|A\|_2}\). The auxiliary certificate state is recomputed from the current residual at each training epoch, as described in Section 5. The lower certificate is not optimized; it is only reported after training as a finite-probe diagnostic.

Pan/Zoom

Pan/Zoom

Pan/Zoom

Figure 5, Figure 6, Figure 7 show that the one-sided upper-certificate branch is the most effective of the three continuations. Its final error is \(5.68\times10^{-4}\), compared with \(4.43\times10^{-3}\) for the standard continuation and \(9.79\times10^{-3}\) for the Lipschitz-refined branch. For the one-sided trained model, the final one-sided upper certificate is \(2.26\times10^{-2}\), whereas the final Lipschitz-based upper certificate is about \(5.70\times10^7\). This supports the practical recommendation to use the one-sided upper certificate for training whenever a reliable one-sided constant is available.

We developed a computable two-sided a posteriori certification framework for PINN approximations of ordinary differential equations. In contrast to certification results based only on global Lipschitz continuity, which provide upper error bounds, the present approach gives both lower and upper certificates under local one-sided growth assumptions. In the nonlinear setting, the bounds are governed by admissible lower and upper one-sided constants on a certified state-space domain. In the linear autonomous case, these constants reduce to the sharp symmetric-part quantities

The simultaneous availability of lower and upper bounds gives a geometric localization of the exact solution around the PINN prediction. Namely, if

then the exact state belongs to

Thus the certificate defines a union of intervals in one dimension, an annulus in two dimensions, a spherical shell in three dimensions, and a closed annular region in general dimension. The band width \( \overline E(t)-\underline E(t)\) therefore quantifies the sharpness of the localization.

The numerical examples illustrate these conclusions in three regimes. The nonlinear radial-growth plus rotation example shows that localization and one-sided constants are informative in a genuinely nonlinear multidimensional system; the rotational term enlarges the Lipschitz constant without affecting the symmetric-part growth. The high-dimensional stiff linear example displays the six admissible linear constants and shows that a finite-probe signed-residual lower certificate can remain informative when the scalar lower certificate is zero. The unstable oscillatory training example shows how the upper certificate can be used as a refinement term, with the lower certificate retained as a post-training diagnostic.

Limitations and admissible model classes