Residual-based error bound for physics-informed neural networks
Neural networks are universal approximators and are studied for their use in
solving differential equations. However, a major criticism is the lack of error
bounds for obtained solutions. This paper proposes a technique to rigorously
evaluate the error bound of Physics-Informed Neural Networks (PINNs) on most
linear ordinary differential equations (ODEs), certain nonlinear ODEs, and
first-order linear partial differential equations (PDEs). The error bound is
based purely on equation structure and residual information and does not depend
on assumptions of how well the networks are trained. We propose algorithms that
bound the error efficiently. Some proposed algorithms provide tighter bounds
than others at the cost of longer run time.
Shuheng Liu, Xiyue Huang, Pavlos Protopapas
# Residual-Based Error Bound for Physics-Informed Neural Networks

###### Abstract

Neural networks are universal approximators and are studied for their use in solving differential equations. However, a major criticism is the lack of error bounds for obtained solutions. This paper proposes a technique to rigorously evaluate the error bound of Physics-Informed Neural Networks (PINNs) on most linear ordinary differential equations (ODEs), certain nonlinear ODEs, and first-order linear partial differential equations (PDEs). The error bound is based purely on equation structure and residual information and does not depend on assumptions of how well the networks are trained. We propose algorithms that bound the error efficiently. Some proposed algorithms provide tighter bounds than others at the cost of longer run time.

## 1 Introduction

Differential equations (DEs) are a useful mathematical tool for describing various phenomena in natural sciences, engineering, and humanity studies. As universal approximators, neural networks are powerful in approximating unknown functions. With back-propagation and modern computing devices, neural networks are convenient to differentiate, making them an ideal choice for solving differential equations.

However, a major criticism of neural network solutions to DEs is the lack of error bound. Traditional numerical methods, such as the finite difference method (FDM) and the finite element method (FEM), compute numerical solutions with known error bounds. Unlike traditional methods, the error bounds of neural network solutions are not well-studied. Therefore, solving DEs with neural networks requires ad hoc customization and empirical hyperparameter finetuning. If the error of _any_ given network can be bounded, we can train neural networks until the error falls below a specified tolerance threshold.

Our contribution is that we propose rigorous error-bounding algorithms for any neural network solution to certain classes of equations, including linear ODEs, certain nonlinear ODEs, and first-order linear PDEs. These algorithms can also be extended to bound the error of other classes of equations as well. The proposed algorithms only use residual information and equation structure as inputs and do not rely on assumptions of finetuning.

Section 2 introduces the symbols and notations adopted in this paper. Section 3 reviews the literature on emerging areas of research that are relevant to solving DEs with neural networks. Section 4 explains the existing effort to bounding the error of neural network DE solutions. Sections 5 and 6 propose various algorithms for the error bound of ODEs and PDEs, respectively. Section 7 uses the method of manufactured solution to verify the validity of each error-bounding algorithm and provides visualization of the tightness of the bounds.

## 2 Symbols and Notations

DEs in this paper are posed w.r.t. unknown function \(v\),

\[\mathcal{D}v=f,\]

where \(\mathcal{D}\) is a possibly nonlinear differential operator and \(f\) is some forcing function. Unlike the exact solution \(v(\cdot)\), a neural network solution \(u(\cdot)\) does not strictly satisfy the equation. Instead, it incurs an additional residual term, \(r\), which the network aims to minimize, to the equation,

\[\mathcal{D}u=f+r.\]

The input to \(v\), \(u\), \(f\), and \(r\) is time \(t\) for ODEs and spatial coordinates \((x,y)\) for PDEs. We limit our reasoning to 2-dimensional PDEs in this work. In cases with multiple unknown functions, we use vector notations \(\mathbf{v}\), \(\mathbf{u}\), and \(\mathbf{r}\) instead of the scalar notations \(v\), \(u\), and \(r\).
The loss function of the network solution is defined as the \(L^{2}\) norm of residual \(r\) over the domain of interest,

\[\mathrm{Loss}(u):=\frac{1}{|I|}\int_{I}\|r\|^{2}\mathrm{d}I=\frac{1}{|I|}\int_{I }\|\mathcal{D}u-f\|^{2}\mathrm{d}I, \tag{1}\]

where a spatial domain \(\Omega\) is substituted for the temporal domain \(I\) in the case of a PDE.

### Initial and Boundary Conditions

For a neural network to satisfy initial or boundary conditions, we apply a technique called _parametrization_. As an intuitive example, the parametrization \(u(t)=(1-e^{-t})\mathrm{Net}(t)+v(0)\) guarantees that \(u(t)\) satisfies the initial condition \(u(0)=v(0)\) regardless of the network \(\mathrm{Net}(\cdot)\). This does not affect the capability of \(\mathrm{Net}(\cdot)\) to learn any solution.

The parametrization is more complicated for higher-order ODEs and most PDEs and has been extensively studied by Lagaris et al. (1998); Lagaris et al. (2000); McFall and Mahan (2009); Lagari et al. (2020), and Sukumar and Srivastava (2021). In this work, we assume all initial and boundary conditions are exactly satisfied.

### Error and Error Bound

The error of a network solution \(u\) is defined as

\[\eta:=u-v. \tag{2}\]

We are interested in _bounding_ the error with a scalar function \(\mathcal{B}\) such that

\[\|\eta(t)\|\leq\mathcal{B}(t)\quad\text{or}\quad\|\eta(x,y)\|\leq\mathcal{B}( x,y) \tag{3}\]

where \(\|\eta\|=\|u-v\|\) is the _absolute error_. If \(\mathcal{B}\) takes on the same value \(B\in\mathbb{R}^{+}\) over the domain, it can be replaced with a constant \(B\).

Notice that multiple bounds \(\mathcal{B}\) exist for the same network solution \(u\). For example, \(|\eta(t)|\leq\mathcal{B}^{(1)}(t)\leq\mathcal{B}^{(2)}(t)\leq\cdots\leq B\) are bounds in decreasing order of tightness. Tighter bounds incur a higher computational cost, and looser bounds (such as constant \(B\)) are faster to compute.

A summary of the applicability, restraints, run-time complexity, and relative tightness of all proposed algorithms is listed in Table 1.

## 3 Literature Review

Hornik et al. (1989) showed that neural networks are universal function approximators. Lagaris et al. (1998) first studied the application of neural networks in solving DEs. The term _physics-informed neural networks_, or PINNs, was first introduced by Raissi et al. (2019) to name neural networks that satisfy DEs while fitting observed data points. Although we train PINNs only to solve DEs without any observed data in this work, the error-bounding algorithms we propose work for any given neural network, regardless of the training process.

Flamant et al. (2020) and Desai et al. (2021) showed that one main advantage of neural networks over traditional numerical methods, such as FDM and FEM, is that neural networks can potentially learn the structure of the solution space and give a bundle of solutions \(u(\mathbf{x};\Theta)\) for different equation setup and initial/boundary conditions parameterized by \(\Theta\). For traditional methods, a new solution must be recomputed for any slight changes in equation setup or initial/boundary conditions.

Some effort has been made to redefine the objective loss function. Yu et al. (2017) applied the Ritz method to a particular class of variational problems. Mattheakis et al. (2020) incorporated an additional constraint to force the network to learn solutions with energy conservation. Parwani and Protopapas (2021) used an adversarial network for sampling in particular areas of the domain where the residual is large.

There are also works that study the failure modes of PINNs and quantify the error of PINN solutions in recent years. Graf et al. (2021) worked on quantifying the uncertainty of PINNs using the Bayesian framework. Krishnapriyan et al. (2021) characterized possible failure modes of PINNs by studying the performance of PINNs on simple problems and analyzing their loss landscape. Krishnapriyan et al. (2021) also concluded that optimization difficulty is the essential cause of failure.

Our work uncovers the mathematical relationship between residual information and the error of PINNs on several classes of ODEs and PDEs. We propose different algorithms for various classes of equations and experimentally validate these algorithms.

## 4 Existing Work

Sirignano and Spiliopoulos (2018) showed that for a class of quasi-linear parabolic PDEs, a neural network with a single hidden layer and sufficiently many hidden units could arbitrarily approximate the exact solutions. Guo and Haghighat (2022) proposed an energy-based _constitutive relation error_ bound for elasticity problems.

De Ryck and Mishra (2022) derived an error bound for ReLU networks on parametric hyperbolic conservation laws. De Ryck and Mishra (2022) showed that there exists some PINN with arbitrarily small residual for Kolmogorov PDEs. De Ryck and Mishra (2022) derived an error bound for operator learning with PINNs. The works of De Ryck and Mishra mentioned above did not bound the error of every given network. Instead, they mathematically proved the existence of a network with errors below a specified bound,
under certain assumptions of network architecture, including width, depth, and activation functions. The question remaining to be answered is how to overcome optimization difficulties and find such a neural network.

Our work differs from the above in that we bound the error of _any_ neural network regardless of finetuning, even networks with randomly initialized weights. Our algorithms only depend on inputs of residual information \(r\), often used as training loss, and equations structure \(\mathcal{D}v=f\). The output is a (possibly constant) function that guarantees to bound the error at any point in domain.

## 5 Error Bound for ODE

This section considers both linear and nonlinear ODEs over the temporal domain \(I=[0,T]\). Initial conditions are imposed on \(\frac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}v(t=0)\) for \(k=0,\ldots,(n-1)\), where \(n\) is the highest order of derivative terms in ODE.

### Error Bound for Linear ODE

Consider the linear ODE \(\mathcal{L}v(t)=f(t)\), where \(\mathcal{L}\) is a linear differential operator. Its neural network solution \(u\) satisfies \(\mathcal{L}u(t)=f(t)+r(t)\). Since error \(\eta:=u-v\), there is

\[\mathcal{L}\eta(t)=r(t). \tag{4}\]

With the assumption in Section 2.1 that \(u\) satisfies the initial conditions at \(t=0\), there is

\[\eta(0)=0,\quad\frac{\mathrm{d}}{\mathrm{d}t}\eta(0)=0,\quad\frac{\mathrm{d} ^{2}}{\mathrm{d}t^{2}}\eta(0)=0,\quad\ldots \tag{5}\]

With initial conditions 5 known, a unique inverse transform \(\mathcal{L}^{-1}\)to \(\mathcal{L}\) exists. Applying \(\mathcal{L}^{-1}\) to Eq. 4, there is

\[\eta(t)=\mathcal{L}^{-1}r(t). \tag{6}\]

Hence, bounding the absolute error \(|\eta|\) is equivalent to bounding \(\big{|}\mathcal{L}^{-1}r\big{|}\). Notice that only a) the equation structure \(\mathcal{L}\) and b) the residual information \(r\) are relevant to estimating the error bound. All other factors, including parameters of the neural network \(u\), forcing function \(f\), and initial conditions, do not affect the error bound at all.

#### 5.1.1 Single Linear ODE with Constant Coefficients

Consider the case where \(\mathcal{L}=\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}+\sum_{j=0}^{n-1}a_{j}\frac {\mathrm{d}^{j}}{\mathrm{d}t^{j}}\) consists of only constant coefficients \(a_{0},a_{1},\ldots,\in\mathbb{R}\). Its characteristic polynomial (defined below) can be factorized into

\[\lambda^{n}+a_{n-1}\lambda^{n-1}+\cdots+a_{0}=\prod_{j=1}^{n}(\lambda-\lambda _{j}), \tag{7}\]

where \(\lambda_{1},\ldots,\lambda_{n}\in\mathbb{C}\) are the characteristic roots.

It can be shown that, for a semi-stable system (\(\mathcal{R}e\left(\lambda_{j}\right)\leq 0\) for all \(\lambda_{j}\)), an error bound can be formulated as

\[|\eta(t)|\leq\mathcal{B}_{loose}(t):=C_{\lambda_{1:n}}\;R_{\max}\,t^{Z}, \tag{8}\]

where \(0\leq Z\leq n\) is the number of \(\lambda_{j}\) whose real part is \(0\), \(C_{\lambda_{1:n}}:=\frac{1}{Z!}\prod_{j=1;\lambda_{j}\neq 0}^{n}\frac{1}{ \mathcal{R}e(-\lambda_{j})}\) is a constant coefficient, and \(R_{\max}:=\max_{t\in I}|r(t)|\) is the maximum absolute residual. Knowing bound 8 is sufficient to qualitatively estimate the error for applications where only the order of error is concerned. See Alg. 1 for reference.

An issue with Eq. 8 and Alg. 1 is that they assume \(\mathcal{R}e\left(\lambda_{j}\right)\leq 0\) for all characteristic roots \(\lambda_{j}\). To address this issue, we propose an alternative error-bounding Alg. 2, which requires more computation but does not require the system to be semi-stable and provides a tighter bound.

Notice that the bounds of \(\eta\) in Eq. 6 can be estimated if the inverse operator \(\mathcal{L}^{-1}\) is known. Let Eq. 7 be the factorization of characteristic polynomial of \(\mathcal{L}\). Define operator \(\mathcal{I}_{\lambda}\) as 1

Footnote 1: This paper assumes the network solution exactly satisfies the initial conditions as discussed in Section 2.1. However, all our algorithms can be extended to cases where the network solution differs from the exact solution by some value. This is achieved by replacing \(\mathcal{I}_{\lambda}\psi(t)\) with \(\mathcal{I}_{\lambda,\delta}\psi(t)=\mathcal{I}_{\lambda}\psi(t)+\delta e^{ \lambda t}\) in Eq. 9.

\[\mathcal{I}_{\lambda}\psi(t):=e^{\lambda t}\int_{0}^{t}e^{-\lambda\tau}\psi(\tau )\mathrm{d}\tau,\quad\forall\psi:I\to\mathbb{C}. \tag{9}\]

We show in supplementary material that \(\mathcal{L}^{-1}=\mathcal{I}_{\lambda_{n}}\circ\mathcal{I}_{\lambda_{n-1}}\circ \cdots\circ\mathcal{I}_{\lambda_{1}}\) and that \(|\mathcal{I}_{\lambda}\psi|\;\leq\mathcal{I}_{\mathcal{R}e(\lambda)}|\psi|\) for any \(\lambda\in\mathbb{C}\) and function \(\psi\). Hence, another error bound can be formulated as

\[\mathcal{B}_{tight}(t):=\left(\mathcal{I}_{\mathcal{R}e(\lambda_{n})}\circ \cdots\circ\mathcal{I}_{\mathcal{R}e(\lambda_{1})}\right)|r(t)|. \tag{10}\]

\begin{table}
\begin{tabular}{l c c c c} \hline \hline
**Algorithm** & **Applicable to** & **Restraint** & **Run-Time** & **Comment** \\ \hline Algorithm 1 & Linear ODE & Semi-stable & \(O(L)\) & Looser than Alg 2 \\ Algorithm 2 & Linear ODE & & \(O(nL)\) & Tighter than Alg 1 \\ Algorithm 3 & Linear ODE System & & \(O(n^{3}L)\) & Norm and elementwise bounds \\ Algorithm 4 & Nonlinear ODE & Nonlinear term is \(\varepsilon v^{k}\) & \(O(JnL)\) & Bounded solution for family of DEs \\ Algorithm 5 & Linear 1st-Order PDE & Coeff. \(c\neq 0\) over domain & \(O(\text{mesh})\) & Constant bound; Looser than Alg 6 \\ Algorithm 6 & Linear 1st-Order PDE & Solvable characteristics & \(O(KL)\) & Tighter than Alg 5 if computable \\ \hline \hline \end{tabular}
\end{table}
Table 1: **Overview of Proposed Algorithms.** The symbols in run-time analysis are defined and explained in detail in Sections 5 and 6, with the exception of \(K\), which is the number of steps used in each numerical integration.
It is also proven in supplementary material that \(\mathcal{B}_{tight}\) is tighter than \(\mathcal{B}_{loose}\) when \(\mathcal{B}_{loose}\) is applicable,

\[|\eta(t)|\leq\mathcal{B}_{tight}(t)\leq\mathcal{B}_{loose}(t)\quad\forall t\in I. \tag{11}\]

Based on Eq. 10, we propose Alg. 2, which computes \(\mathcal{B}_{tight}\) by repeatedly evaluating integrals in 9 using the cumulative trapezoidal rule.

``` Input: Coefficients \(\{a_{j}\}_{j=0}^{n-1}\) for operator \(\mathcal{L}\), residual information \(r(\cdot)\), domain of interest \(I=[0,T]\), and a sequence of time points \(\{t_{\ell}\}_{\ell=1}^{L}\) where error bound is to be evaluated Output: Error bound at given time points \(\{\mathcal{B}(t_{\ell})\}_{\ell=1}^{L}\) \(\mathcal{L
meantime, a _norm bound_ (scalar) \(\mathcal{B}(t)\) also exists

\[\|\boldsymbol{\eta}(t)\|\leq\mathcal{B}(t):=\mathrm{cond}(P)\left\|\boldsymbol{ \mathcal{I}}\big{[}\|\mathbf{r}\|\mathbf{1}\big{]}(t)\right\| \tag{16}\]

where \(\mathrm{cond}(P)\) is the conditional number of \(P\) w.r.t. induced matrix norm and \(\mathbf{1}\) is an \(n\times 1\) column vector of \(1\)s. Proof of Eq. 15 and Eq. 16 can be found in in supplementary material. See Alg. 3 for implementation.

**Input:** Coefficient matrix \(A\in\mathbb{R}^{n\times n}\), residual vector \(\mathbf{r}(t)\), and a sequence of points \(\{t_{\ell}\}_{\ell=1}^{L}\) where error is to be bounded

**Output:** Norm bound (scalar) \(\{\mathcal{B}(t_{\ell})\}_{\ell=1}^{L}\) and component-wise bound (vector) \(\{\boldsymbol{\mathcal{B}}(t_{\ell})\}_{\ell=1}^{L}\) at given time points

**Ensure:**: \(\|\eta(t_{\ell})\|\leq\mathcal{B}(t_{\ell})\) and \(\eta(t_{\ell})\preceq\boldsymbol{\mathcal{B}}(t_{\ell})\) for all \(\ell\)

\(J,P\gets\) Jordan canonicalization of \(A=PJP^{-1}\)

**for** each Jordan block \(J_{k}\) of shape \(n_{k}\times n_{k}\)**do**

\(\mathbf{I}_{k}\leftarrow\) construct operator block using Eq. 14

**end for**

\(\boldsymbol{\mathcal{I}}\leftarrow\) diag\((\mathbf{I}_{1},\mathbf{I}_{2},\ldots)\)

\(\{\boldsymbol{\mathcal{B}}(t_{\ell})\}_{\ell=1}^{L}\leftarrow\{P^{|\cdot|} \boldsymbol{\mathcal{I}}\big{[}(P^{-1})^{|\cdot|}\boldsymbol{\cdot}^{|\cdot|} \big{]}(t_{\ell})\}_{\ell=1}^{L}\)

\(\{\mathcal{B}(t_{\ell})\}_{\ell=1}^{L}\leftarrow\{\mathrm{cond}(P)\left\| \boldsymbol{\mathcal{I}}\big{[}\|\mathbf{r}\|\mathbf{1}\big{]}(t_{\ell})\right\| \}_{\ell=1}^{L}\)

**return**\(\{\mathcal{B}(t_{\ell})\}_{\ell=1}^{L}\), \(\{\boldsymbol{\mathcal{B}}(t_{\ell})\}_{\ell=1}^{L}\)

**Algorithm 3** ODE System Bound (norm and elementwise)

### Nonlinear ODE

Nonlinear ODEs are hard to solve in general. In this work, we only deal with nonlinear ODEs with a single nonlinear term of the form \(\varepsilon v^{k}(t)\), where \(\varepsilon\in\mathbb{R}\) is a small number. Ideally, \(|\varepsilon|\ll 1\). The exact requirement for \(\epsilon\) is given in Section 5.2.2. The value of \(\varepsilon\) can vary within a certain range or be fixed. With the perturbation technique, we obtain a family of solutions \(v(t;\varepsilon)\) parameterized by \(\varepsilon\) at the cost of solving a (countable) collection of equations. As explained below in Section 5.2.1, we train finitely many networks, each approximately solving an equation in the collection.

#### 5.2.1 Perturbation Theory

Consider the nonlinear ODE with nonlinear term \(\varepsilon v^{k}(t)\),

\[\mathcal{L}v(t)+\varepsilon v^{k}(t)=f(t), \tag{17}\]

where \(\mathcal{L}\) is a linear differential operator discussed in 5.1 and initial conditions are specified for the system at time \(t=0\). Notice that each \(\varepsilon\in\mathbb{R}\) corresponds to a solution \(v(t;\varepsilon)\). We expand the solution \(v(t;\varepsilon)\) in terms of \(\varepsilon\)

\[v(t;\varepsilon)=\sum_{j=0}^{\infty}\varepsilon^{j}v_{j}(t)=v_{0}(t)+ \varepsilon v_{1}(t)+\ldots \tag{18}\]

Only \(v_{0}(t)\) is subject to the original initial conditions at \(t=0\), while other components, \(v_{1}\), \(v_{2}\),..., have initial conditions of \(0\) at \(t=0\). Substituting Eq. 18 into Eq. 17,

\[\mathcal{L}\sum_{j=0}^{\infty}\varepsilon^{j}v_{j}+\varepsilon \left(\sum_{j=0}^{\infty}\varepsilon^{j}v_{j}\right)^{k}=f \tag{19}\] \[\sum_{j=0}^{\infty}\varepsilon^{j}\mathcal{L}v_{j}+\sum_{j=0}^{ \infty}\varepsilon^{j+1}\sum_{\begin{subarray}{c}j_{1}+\cdots j_{k}=j\\ j_{1},\ldots,j_{k}\geq 0\end{subarray}}v_{j_{1}}\ldots v_{j_{k}}=f\] (20) \[\mathcal{L}v_{0}+\sum_{j=1}^{\infty}\varepsilon^{j}\bigg{(} \mathcal{L}v_{j}+\sum_{\begin{subarray}{c}j_{1}+\cdots+j_{k}=j-1\\ j_{1},\ldots,j_{k}\geq 0\end{subarray}}v_{j_{1}}\ldots v_{j_{k}}\bigg{)}=f. \tag{21}\]

In order for Eq. 21 to hold true for all \(\varepsilon\), the coefficients for each \(\varepsilon^{j}\) must match on both sides of Eq. 21. Hence,

\[\mathcal{L}v_{0} =f \tag{22}\] \[\mathcal{L}v_{1}+v_{0}^{k} =0\] (23) \[\mathcal{L}v_{2}+kv_{0}^{k-1}v_{1} =0\] (24) \[\mathcal{L}v_{3}+\frac{k(k-1)}{2}v_{0}^{k-2}v_{1}^{2}+kv_{0}^{k-1 }v_{2} =0\] (25) \[\vdots\]

For \(\varepsilon=0\), Eq. 18 is reduced to \(v_{0}(t)\), which solves the linear problem \(\mathcal{L}v=f\).

The above system can be solved in a _sequential_ manner, either analytically or using neural networks,

1. Eq. 22 is linear in \(v_{0}\) and can be solved first.
2. With \(v_{0}\) known, Eq. 23 is linear in \(v_{1}\) and can be solved for \(v_{1}\).
3. Similarly, with \(v_{0}\) and \(v_{1}\) known, Eq. 24 is linear in \(v_{2}\) and can be solved for \(v_{2}\).
4. The process can be repeated for Eq. 25 and beyond. Only a linear ODE is solved each time.

To solve the system with PINNs, we approximate exact solutions \(\{v_{j}(t)\}_{j=1}^{\infty}\) with neural network solutions \(\{u_{j}(t)\}_{j=0}^{J}\) trained sequentially on Eq. 22, Eq. 23, and beyond. In practice, we only consider components up to order \(J\) to avoid the infinity in expansion 18. Ideally, \(J\) should be large enough so that higher order residuals in expansion 18 can be neglected.

After obtaining \(\{u_{j}(t)\}_{j=0}^{J}\), we can reconstruct the solution \(u(t;\varepsilon)=\sum_{j=0}^{J}\varepsilon^{j}u_{j}(t)\) to the original nonlinear equation 17 for varying \(\varepsilon\). See Alg. 4 for details.

#### 5.2.2 Expansion of Bounds

The absolute error \(|\eta(t;\varepsilon)|=|u(t;\varepsilon)-v(t;\varepsilon)|\)is given by

\[|\eta(t;\varepsilon)|=\left|\sum_{j=0}^{J}\varepsilon^{j}\Big{(}u_{j}(t)-v_{j}(t) \Big{)}-\sum_{j=J+1}^{\infty}\varepsilon^{j}v_{j}(t)\right|\]
\[\leq\sum_{j=0}^{J}\Big{|}\eta_{j}(t)\Big{|}|\varepsilon|^{j}+\left|\sum_{j=J+1}^{ \infty}\varepsilon^{j}v_{j}(t)\right| \tag{26}\]

where \(\eta_{j}(t):=u_{j}(t)-v_{j}(t)\) is the _component error_ between \(u_{j}(t)\) and \(v_{j}(t)\). Let \(\mathcal{B}_{j}\) denote the _bound component_ such that \(|\eta_{j}(t)|\leq\mathcal{B}_{j}(t)\). Assuming \(J\) is large and \(\varepsilon\) is small such that higher order terms \(\left|\sum_{j=J+1}^{\infty}\varepsilon^{j}v_{j}(t)\right|\) are negligible, there is

\[\Big{|}\eta(t;\varepsilon)\Big{|}\leq\mathcal{B}(t;\varepsilon):=\sum_{j=0}^{J }\mathcal{B}_{j}(t)\,|\varepsilon|^{j} \tag{27}\]

where each bound component \(\mathcal{B}_{j}\) can be evaluated using the techinque in Section 5.1. See Alg. 4 for details.

**Input:** Linear operator \(\mathcal{L}\), nonlinear degree \(k\), domain \(I=[0,T]\), highest order \(J\) for expansion, and a sequence \(\{(t_{\ell},\varepsilon_{\ell})\}_{\ell=1}^{L}\) where solution \(u(t;\varepsilon)\) and error bound \(\mathcal{B}(t;\varepsilon)\) are to be evaluated

**Output:** Solution \(\{u(t_{\ell};\varepsilon_{\ell})\}_{\ell=1}^{L}\) and error bound \(\{\mathcal{B}(t_{\ell};\varepsilon_{\ell})\}_{\ell=1}^{L}\)

**Require:**\(t_{\ell}\in I\), and \(|\varepsilon_{\ell}|\) to be small (ideally \(|\varepsilon_{\ell}|\ll 1\))

**Ensure:**\(\eta(t_{\ell};\varepsilon_{\ell})\leq\mathcal{B}(t_{\ell};\varepsilon_{\ell})\)

\(u_{0},r_{0},\leftarrow\) net solution, residual of \(\mathcal{L}u_{0}=\underset{\{\mathcal{B}(0_{\ell})\}_{\ell=1}^{L}}{\mathcal{B }(t_{\ell})}\leftarrow\) bound of \(\left|\mathcal{L}^{-1}r_{0}\right|\) at \(\{t_{\ell}\}_{\ell=1}^{L}\)

**for**\(j\gets 1\ldots J\)**do**

Macro \(\text{NL}_{j}[\phi]\leftarrow\sum_{\begin{subarray}{c}j_{1}+\cdots+j_{k}=j-1 \end{subarray}}\phi_{j_{1}}\ldots\phi_{j_{k}}\)

\(u_{j},r_{j}\leftarrow\) net solution, residual of \(\mathcal{L}u_{j}+\text{NL}_{j}[u]=0\)

\(\mathcal{B}_{\text{NL}}\leftarrow\) upper bound of \(\left|\text{NL}_{j}[u]-\text{NL}_{j}[v]\right|\)

\(\left\{\mathcal{B}_{j}(t_{\ell})\right\}_{\ell=1}^{L}\leftarrow\) bound of \(\left|\mathcal{L}^{-1}r_{j}\right|+\left|\mathcal{L}^{-1}\mathcal{B}_{\text{NL}}\right|\)

**end for**

\(\{u(t_{\ell};\varepsilon_{\ell})\}_{\ell=1}^{L}\leftarrow\)\(\big{\{}\sum_{j=0}^{J}\varepsilon_{\ell}^{j}u_{j}(t_{\ell})\big{\}}_{\ell=1}^{L}\)

\(\{\mathcal{B}(t_{\ell};\varepsilon_{\ell})\}_{\ell=1}^{L}=\)\(\big{\{}\sum_{j=0}^{J}\varepsilon_{\ell}^{j}\mathcal{B}_{j}(t_{\ell})\big{\}}_{ \ell=1}^{L}\)

**return \(\{u(t_{\ell};\varepsilon_{\ell})\}_{\ell=1}^{L},\{\mathcal{B}(t_{\ell}; \varepsilon_{\ell})\}_{\ell=1}^{L}\)**

**Note 1**: \(\mathcal{B}_{0}\) and \(\mathcal{B}_{1:J}\) can be evaluated using Alg. 1 or 2.

**Note 2**: \(\mathcal{B}_{\text{NL}}\) can be estimated even though exact solutions \(v_{0:j-1}\) are unknown. This is because \(v_{i}\in[u_{i}-\mathcal{B}_{i},u_{i}+\mathcal{B}_{i}]\) for all \(i\), and \(u_{0:j-1}\), \(\mathcal{B}_{0:j-1}\) are known from previous iterations.

## 6 Error Bound for PDE

This section considers first-order linear PDEs defined on a 2-dimensional spatial domain \(\Omega\),2

Footnote 2: Similar techniques can be used for other classes of PDEs and higher dimensions where the method of characteristics applies.

\[a(x,y)\frac{\partial v}{\partial x}+b(x,y)\frac{\partial v}{\partial y}+c(x,y)v =f(x,y) \tag{28}\]

with Dirichlet boundary constraints defined on \(\Gamma\subset\partial\Omega\),

\[v\big{|}_{(x,y)\in\Gamma}=g(x,y). \tag{29}\]

We partition the domain into infinitely many characteristic curves \(\mathcal{C}\), each passing through a point \((x_{0},y_{0})\in\Gamma\). The resulting curve is a parameterized integral curve

\[\mathcal{C}:\begin{cases}x^{\prime}(s)=a(x,y)\\ y^{\prime}(s)=b(x,y)\end{cases}\text{where}\,(\cdot)^{\prime}=\frac{\text{d}}{ \text{d}s}\text{ and }\begin{array}{c}x(0)=x_{0}\\ y(0)=y_{0}.\end{array}\]

For any \((x(s),y(s))\) on \(\mathcal{C}\), functions \((v,a,b,c,f)\) can be viewed as univariate functions of \(s\). By chain rule, there is

\[a(x,y)\frac{\partial v}{\partial x}+b(x,y)\frac{\partial v}{\partial y}=x^{ \prime}(s)\frac{\partial v}{\partial x}+y^{\prime}(s)\frac{\partial v}{\partial y }=v^{\prime}(s).\]

Hence, Eq. 28 is reformulated as an ODE along curve \(\mathcal{C}\),

\[v^{\prime}(s)+c(s)v(s)=f(s)\quad\text{s.t. }v(0)=g(x_{0},y_{0}), \tag{30}\]

where \(v(s)\), \(c(s)\), and \(f(s)\) are shorthand notations for \(v(x(s),y(s))\), \(c(x(s),y(s))\), and \(f(x(s),y(s))\), respectively.

In particular, if \(c(x,y)\neq 0\) for all \((x,y)\in\Omega\), both sides of Eq. 28 can be divided by \(c(x,y)\), resulting in a residual of \(r(x,y)/c(x,y)\) where \(r(x,y)\) is the residual of the original problem. By Eq. 8, a constant error bound on \(\mathcal{C}\) is \(|\eta(s)|\leq\max_{s}|r(s)/c(s)|\). Hence, a (loose) constant error bound \(B\) (see Alg. 5) over the entire domain \(\Omega\) is

\[|\eta(x,y)|\leq B:=\max_{(x,y)\in\Omega}\left|\frac{r(x,y)}{c(x,y)}\right|. \tag{31}\]

**Algorithm 5** Constant Err Bound for Linear 1st-Order PDE

**Input:** Coefficient \(c(x,y)\) in Eq. 28, residual information \(r(x,y)\) and domain of interest \(\Omega\)

**Output:** A constant error bound \(B\in\mathbb{R}^{+}\)

**Require:**\(c(x,y)\neq 0\) for all \((x,y)\in\Omega\)

**Ensure:**\(|\eta(x,y)|\leq B\) for all \((x,y)\in\Omega\)

\(\{(x_{k},y_{k})\}_{k}\leftarrow\) sufficiently dense mesh grid over \(\Omega\)

\(B\leftarrow\max_{k}\left|\frac{r(x_{k},y_{k})}{c(x_{k},y_{k})}\right|\)

**return \(B\)**

**Algorithm 6** Constant Err Bound for Linear 1st-Order PDE

## 7 Relevant Experiments

In this section, we perform experiments on equations with manufactured solutions using the _NeuroDiffEq_ library [Chen et al., 2020], which provides convenient tools for training PINNs.
First, we train networks to solve equations and collect their residual information \(r\). Then, we apply Alg. 1-6 (where applicable) to derive error bounds using only residual information \(r\) and equation structure, characterized by its differential operator \(\mathcal{D}\). Lastly, we show that the absolute error strictly falls within the bounds, regardless of how well the networks are trained.

Throughout this section, we always use networks with two hidden layers, each consisting of 32 hidden units. Depending on whether the problem is an ODE or PDE, a network can have a single input \(t\) or two inputs \((x,y)\), but always have a single output. The activation function is \(\tanh\). Unless otherwise noted, the training domain is \(I=[0,1]\) for ODEs and \(\Omega=[0,1]^{2}\) for PDEs. We use a _PyTorch_ Adam optimizer with default hyperparameters to train networks for 1000 epochs.

Notice that we list these configurations only for the reproducibility of visualizations. Our error-bounding algorithm works under any other configurations.

### Single Linear ODE with Constant Coefficients

Here, we study three equations \(v^{\prime\prime}+3v^{\prime}+2v=f(t)\), \(v^{\prime\prime}+v=g(t)\), and \(v^{\prime\prime}-v^{\prime}=h(t)\), whose characteristic roots are \(\{-1,-2\}\), \(\{\pm i\}\), and \(\{0,1\}\) respectively. By Section 5.1.1, the first two equations can be bounded with either Alg. 1 or Alg. 2, while the last must be bounded with Alg. 2.

They all satisfy initial conditions \(v(0)=v^{\prime}(0)=1\). We pick \(f(t)=2t^{2}+8t+7\), \(g(t)=t^{2}+t+3\), and \(h(t)=1-2t\), so that the manufactured solution is \(v(t)=t^{2}+t+1\) for all three equations. Fig. 1 shows that both \(\mathcal{B}_{loose}\) (Alg. 1) and \(\mathcal{B}_{tight}\) (Alg. 2) strictly bounds the absolute error.

### Linear ODE System with Constant Coefficients

In this subsection, we train \(6\) networks to solve a \(6\)-dimensional linear system of ODEs with constant coefficients, namely, \(\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{v}+A\mathbf{v}=\mathbf{f}\). We pick \(A=PJP^{-1}\) where \(J=\begin{pmatrix}J_{1}&J_{2}&\\ &J_{3}&\\ &&J_{1}=\begin{pmatrix}4&1\\ &4&1\\ &&4\end{pmatrix}\), \(J_{2}=\begin{pmatrix}3&1\\ &3\end{pmatrix}\), \(J_{3}=2\), and \(P\) is a random orthogonal matrix.

We pick the initial conditions to be \(\mathbf{v}(0)=P(0\,0\,1\,0\,1\,1)^{T}\) and the forcing function to be \(\mathbf{f}(t)=P(\cos t+4\sin t+\ln(1+t),\,\frac{1}{1+t}+4\ln(1+t)+(t+1),\,4t+5, \,2t+3t^{2}+e^{t},\,4e^{t},\,2\cos t-\sin t)^{T}\), so that the manufactured exact solution is \(\mathbf{v}(t)=P(\sin t,\,\ln(t+1),\,t+1,\,t^{2},\,e^{t},\,\cos t)^{T}\). After obtaining the residual information \(\mathbf{r}(t)=\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{u}(t)+A\mathbf{u}(t)- \mathbf{f}(t)\), we apply Alg. 3 to obtain componentwise bound and norm bound of \(\boldsymbol{\eta}=\mathbf{u}-\mathbf{v}\). It is shown in Fig. 2 that the bounds hold over the domain.

Figure 1: Loose bound (Alg. 1) and tight bound (Alg. 2) for 3 second-order linear ODE with constant coefficients. Notice that the loose bound cannot be applied to the third equation since it has characteristic roots with positive real part.

Figure 2: _Componentwise_ bound (upper) and _norm_ bound (lower) for linear ODE system with constant coefficients
### Nonlinear ODE - Duffing Equation

In this subsection, we consider a Duffing oscillator, which is characterized by the following \(2\)nd order nonlinear ODE:

\[\frac{\mathrm{d}^{2}v}{\mathrm{d}t^{2}}+3\frac{\mathrm{d}v}{\mathrm{d}t}+2v+ \varepsilon v^{3}=\cos t, \tag{32}\]

under initial conditions \(v(0)=1\) and \(v^{\prime}(0)=1\), where \(\varepsilon\) controls the nonlinearity of the equation. Using Alg. 4, we solve the equation on \(I=[0,2]\) for linspaced \(\varepsilon\in(-0.9,0.9)\) using neural networks and bound the errors. The input \(J\) to Alg. 4 is chosen to be \(6\). Namely, we expand the solution and bound components from degree \(0\) to \(6\).

The analytical solution to Eq. 32 is complicated. Hence, we use the RKF4(5) method to compute numerical solutions that are close enough to exact solutions for visualization purposes. See Fig. 3 for network solutions against RKF4(5) solutions and Fig. 4 for error bounds against absolute error.

### Linear PDE System with Nonconstant Coefficients

#### 7.4.1 PDE Error Bound Evaluation Using Alg. 6

We try to solve the following first-order linear PDE,

\[(-x-y)\frac{\partial v}{\partial x}+(x-y)\frac{\partial v}{\partial y}+v=3x-2y \tag{33}\]

in spatial domain \(\Omega=[-1,1]^{2}\). The boundary constraints are \(v(x,\pm 1)=2x\pm 3\) and \(v(\pm 1,y)=3y\pm 2\). The manufactured solution is given by \(v(x,y)=2x+3y\). The characteristic curves are integral curves \(\mathcal{C}:\begin{cases}x^{\prime}(s)=-x-y\\ y^{\prime}(s)=x-y\end{cases}\), or

\[\mathcal{C}:\begin{cases}x(s)=R_{0}e^{-s}\cos(s+\theta_{0})\\ y(s)=R_{0}e^{-s}\sin(s+\theta_{0})\end{cases}\text{, where }R_{0}=\sqrt{x_{0}^{2}+y_{0} ^{2}}\]

and \(\theta_{0}=\mathrm{atan2}(y_{0},x_{0})\) are constants determined by the starting point \((x_{0},y_{0})\in\Gamma=\partial\Omega\). See Figure 5 for visualization.

Since the analytical expression of the characteristic curves is known, Alg. 6 can be applied to evaluate the bound on each curve. We choose \(16\) characteristic curves with starting points \(A\), \(B\),..., \(P\), equidistantly placed on the boundary (Fig. 4(a)). We plot the absolute error and the computed error bound along these characteristic curves in Fig. 6. It can be seen that absolute error lies strictly within the bounds.

#### 7.4.2 PDE Error Bound Evaluation Using Alg. 5

Consider the following PDE

\[(x^{2}+y^{2}+1)\frac{\partial v}{\partial x}+(x^{2}-y^{2}+2)\frac{\partial v }{\partial y}+(3-2x)v=f \tag{34}\]

over domain \(\Omega=[-1,1]^{2}\), where \(f(x,y)=6-4x\). The boundary constraints are \(v(-1,y)=2\) and \(v(x,1)=2\)

Figure 4: True Error vs. error bound (max-degree 0\(\sim\)6) of neural network solution to Duffing Equation 32 for \(\varepsilon\in(-0.9,0.9)\)

Figure 5: Characteristics curves of Eq. 33 (left) and Eq. 34 (right). The red curves, with staring points \(A\) to \(P\), are selected for visualization of absolute error and error bound in Fig. 6.

Figure 3: RKF45 and Network Solutions (max-degree 0\(\sim\)6) to Duffing equation 32 for \(\varepsilon\in(-0.9,0.9)\)
and the manufactured solution is \(v(x,y)=2\).

The characteristic curves \(\mathcal{C}:\begin{cases}x^{\prime}(s)=x^{2}+y^{2}+1\\ y^{\prime}(s)=x^{2}-y^{2}+2\end{cases}\) are given by a nonlinear ODE, which is hard to solve analytically. (See Fig. 5b for visualization) Therefore, Alg. 6 cannot be applied to evaluate the error bound.

However, the coefficient \((3-2x)\) is nonzero over domain \(\Omega\). Hence, we can use Alg. 5 to compute a constant error bound \(|\eta(x,y)|\leq\mathcal{B}(x,y)\equiv B\) for all \((x,y)\in\Omega\). We visualize the bound and the maximum absolute error \(\max_{(x,y)\in\Omega}|\eta|\) after each training epoch in Fig. 7. As expected, the bound is loose, which is about an order of magnitude larger than the max absolute error. Yet, it consistently holds true for every epoch, even during the early stages of training, when the network performs poorly.

## 8 Conclusion and Future Work

This paper proposes various error-bounding algorithms for any PINN solution to certain classes of ODEs and PDEs. These algorithms only require the residual information \(r(\cdot)\) and the equation structure \(\mathcal{D}v=f\) as input. There are many real-world applications for which the exact solution \(v(\cdot)\) is unknown or hard to compute. However, the residual information \(r(\cdot)\) is usually, if not always, readily available. With our proposed algorithms, PINNs can be trained until the error is gauranteed to fall below a specified tolerance threshold. The mathematical relationship between residual and error bound also sheds light on optimizing PINN solutions for future studies.

The error-bounding algorithms proposed in this paper only apply to certain classes of ODEs and PDEs. However, the insights of this paper can be beneficial to future work that extends to more general classes of ODEs and PDEs, especially nonlinear ones. We also plan to apply these algorithms stochastic differential equations, where the error bound is a probabilistic tail bound.
