# Evaluating Error Bound for Physics-Informed Neural Networks on Linear
Dynamical Systems

Shuheng Liu  Institute for Applied Computational Science, Harvard University
Xiyue Huang  Carnegie Mellon University  Pavlos Protopapas  Institute for
Applied Computational Science, Harvard University

(July 2022)

###### Abstract

There have been extensive studies on solving differential equations using
physics-informed neural networks. While this method has proven advantageous in
many cases, a major criticism lies in its lack of analytical error bounds.
Therefore, it is less credible than its traditional counterparts, such as the
finite difference method. This paper shows that one can mathematically derive
explicit error bounds for physics-informed neural networks trained on a class
of linear systems of differential equations. More importantly, evaluating such
error bounds only requires evaluating the differential equation residual
infinity norm over the domain of interest. Our work shows a link between
network residuals, which is known and used as loss function, and the absolute
error of solution, which is generally unknown. Our approach is semi-
phenomonological and independent of knowledge of the actual solution or the
complexity or architecture of the network. Using the method of manufactured
solution on linear ODEs and system of linear ODEs, we empirically verify the
error evaluation algorithm and demonstrate that the actual error strictly lies
within our derived bound.

##  1 Introduction

Differential equations play an essential role in science, engineering, as well
as a wide range of mathematical modeling processes. In many cases, however,
analytical solutions are nonexistent, and therefore many numerical methods
(finite difference, finite volute, finite element, spectral method) have been
studied for these equations. Recently, massive attention has been paid to
solving differential equations with neural networks, which is capable of
learning a continuous solution while occupying relatively less storage. These
networks learn to minimize the squared residual of some differential
equations, which are usually believed to govern some physical systems.

Raissi et al use the term physics-informed neural networks (PINNs) to refer to
networks that are trained on both the differential equation residuals (in an
unsupervised manner) and observed data (in a supervised manner) [1]. However,
in this paper, we only consider networks trained to minimize differential
equation residuals subject to initial/boundary conditions, without any data
from observations. Nonetheless, we still refer to them as physics-informed
neural networks as they are only deprived of observed data while remaining
informed of the underlying physics.

Unlike for tradition numerical methods, little effort has been made to
ascertain the error bound of neural network solutions to differential
equations. Therefore, the reliability of neural network solutions are usually
questionable, making them hardly acceptable for real-world applications.
Additionally, there is little interpretation for the loss function used for
guiding the training process, other than that it should be as close to zero as
possible. However, the link between the loss value (typically an approximation
of the $L_{2}$ norm of ODE residual) and the error bound has not been
formulated.

In this paper, we present a mathematical proof that explicitly bounds the
absolute/relative error of a solution of a class of linear dynamical systems
with neural networks. The method makes no assumption on the network
architecture or whether the network is sufficiently trained. Apart from the
characteristics and structure of the dynamical systems in question, the error
bound we derived ($O(\varepsilon t^{m})$) only depends on time $t$ and the
largest differential equation residual ($\varepsilon$) (infinity norm) over
the domain of interest. We further present that, for strictly stable systems,
one can derive a bound ($O(\varepsilon)$) that is independent of time $t$.
Finally, we present a technique to tighten the error bound by dividing the
time domain into subintervals and evaluating the maximum residual on each one.

Until more straightforward error evaluation approach is identified, our method
can serve as a practical way to estimate the absolute error of neural network,
which is generally unknown, using the network residual, which is easily
computed and used as the loss function.

Our main contributions are threefold. First, we propose an algorithm for
rigorous error bound evaluation. Second, we extend the algorithm to more
complicated linear cases, including linear systems and equation with
nonconstant coefficients. Third, we demonstrated the effectiveness of the
evaluation algorithm using method of manufactured solution.

In Section 2, we introduce the background of physics-informed neural networks
(PINNs) and their application in solving differential equations. In Section 3,
we discuss the previous work in this field. In Section 4, we propose and prove
a method to efficiently evaluate and bound the error. In Section 5, we show
experimental results that verify our proposed method. Finally, we conclusively
summarize our work and provide direction for future work in Sections 6 and 7.

##  2 Background

Neural networks were first introduced as universal function approximators that
learn nonlinear mappings for supervised learning tasks [2]. Lagaris et al.
first proposed solving differential equations using neural networks [3] due to
the smoothness of neural networks with appropriate activation functions.

One main advantage of using neural networks in solving differential equations
is its potential to learn the structure of the solution space and give a
bundle of solutions $\mathbf{u}(t;\theta_{1},\dots,\theta_{n})$ where
$\theta_{1},\dots,\theta_{n}$ parametrize the differential equation and/or
initial condition [4] [5]. For traditional numerical methods, a new solution
must be calculated from scratch whenever any parameter in equation or initial
condition changes. Neural networks, however, can in theory take in
$t,\theta_{1},\dots,\theta_{n}$ as its input and learn the structure of the
solution space parametrized by $\theta_{1},\dots,\theta_{n}$.

Multiple approaches have been proposed to enforce initial/boundary conditions
on network solutions [3] [6] [7] [8] [9]. In this paper, we use a slightly
modified version of Lagaris’ reparametrization technique
$\mathbf{u}(t)=(t-t_{0})\mathbf{NN}(t)+\mathbf{u}_{0}$ to ensure
$\mathbf{u}(t_{0})=\mathbf{u}_{0}$ where $\mathbf{NN}$ is the network. Unless
otherwise specified, we refer to $\mathbf{u}(t)$ as the network solution which
is already reparamtrized to satisfy the initial condition.

To train a network solution $\mathbf{u}(t)$ for a differential equation
$\mathcal{L}\mathbf{u}=\mathbf{f}$, we essentially minimize an approximation
of the $L_{2}$ norm of differential equation residual on a domain $\Omega$

$\int_{I}\left(\mathcal{L}\mathbf{u}-\mathbf{f}\right)^{2}\mathrm{d}t\approx\frac{|I|}{N}\sum_{\begin{subarray}{c}i=1\\\
t_{i}\in
I\end{subarray}}^{N}\left(\mathcal{L}\mathbf{u}(t_{i})-\mathbf{f}(t_{i})\right)^{2}:=\mathrm{Loss}.$

where $\mathcal{L}$ is a (possibly nonlinear) differential operator.

Effort has also been made to redefine the loss function by including energy
constraints or using the weak formulation of equation [10] [11] [12]. By
sampling the training points $\\{t_{i}\\}_{i=1}^{N}\subset I$ in each epoch,
we aim to achieve a residual norm of zero after sufficient training. The
sampling strategy is also an active field of research [13]. However, the
universal approximation theorem relies on the hypothesis of a sufficiently
large network architecture and only asserts the existence of a set of desired
weights. In the context of PINNs, two more relevant questions are 1) how to
interpret the residual and 2) how to bound the error between the network
solution and the exact solution using the residuals.

Little effort has been to study the failure modes and absolute error of
network solutions until recent years [14] [15] [16]. In our work, we
mathematically derive the error bound for a class of linear ODEs, which can be
efficiently computed using only on the infinite norm of ODE residuals on the
domain of interest. We also present a way to tighten the derived error bound
by dividing the domain into smaller subintervals.

##  3 Previous Work

In [de2021error], Ryck and Mishra established a foundation and rationale for
error of PINNs in approximating PDEs. Making Kolmogorov PDEs as an example,
they have shown that there exists PINNs, approximating these PDEs such that
the resulting generalization error and the total error can be made arbitrarily
small. However, the existence of such neural networks does not guarantee
network training converges in practice. A more practical concern is how to
evaluate the error given any network (possibly illy trained) on certain
equations.

##  4 Error Bound Proof and Methodology

Throughout this section, we use $\mathbf{u}:I\to\mathbb{C}^{n}$ to denote the
neural network solution to $\mathcal{L}\mathbf{u}=\mathbf{f}$ where $I$ can be
any of the forms $(t_{0},t_{1}]$, $(t_{0},t_{1})$ or $(t_{0},\infty)$, and
$\mathcal{L}$ is a linear differential operator. In the one-dimensional case,
we use non-bold font $u$ and $f$ instead of $\mathbf{u}$ and $\mathbf{f}$. The
solution residual is defined as
$R\mathbf{u}(t):=\mathcal{L}\mathbf{u}(t)-\mathbf{f}(t)$. The exact solution
$\mathbf{u}^{*}(t)$ satisfies $R\mathbf{u}^{*}(t)\equiv\mathbf{0}$ and the
exact natural response $\mathbf{u}_{n}^{*}(t)$ is defined to satisfy the
associated homogeneous equation $\mathcal{L}\mathbf{u}_{n}^{*}(t)=\mathbf{0}$.
Both $\mathbf{u}^{*}$ and $\mathbf{u}_{n}^{*}$ satisfy the same initial
condition
$\mathbf{u}^{*}(t_{0})=\mathbf{u}_{n}^{*}(t_{0})=\mathbf{u}_{0}^{*}$. We shall
assume that each component of $\mathbf{u^{*}_{0}}$ is nonzero.

###  4.1 First-Order Linear ODE with Constant Coefficients

It is well known that the most general form of first-order linear ODE with
constant coefficients is $u^{\prime}(t)+cu(t)=f(t)$ where $c\in\mathbb{C}$ is
a constant and $u^{\prime}$ is the derivative of $u$.

Proposition If the residual $Ru(t)$ of equation
$u^{\prime}+(\lambda+i\omega)u=f$, where $\lambda,\omega\in\mathbb{R}$, is
bounded by $\varepsilon\geq 0$ on $I$, namely,

$\left|u^{\prime}+(\lambda+i\omega)u-f\right|\leq\varepsilon\quad\forall t\in
I,$ (1)

and the network solution $u$ satisfies initial condition with
$u(t_{0})=u_{0}^{*}\neq 0$, then,

  1. a)

The absolute error is bounded by
$\left|u-u^{*}\right|\leq\cfrac{\varepsilon}{\lambda}\leq O(\varepsilon)$ on
$I$ if the natural response $u_{n}^{*}$ is convergent ($\lambda>0$);

  2. b)

The relative error w.r.t. $u_{n}^{*}$ is bounded by
$\left|\cfrac{u-u^{*}}{u_{n}^{*}}\right|\leq\cfrac{\varepsilon}{-\lambda|u_{0}^{*}|}\leq
O(\varepsilon)$ on $I$ if the natural response $u_{n}^{*}$ is divergent
($\lambda<0$); and

  3. c)

The absolute and relative errors are bounded by $\left|u-u^{*}\right|\leq
O(\varepsilon t)$ and $\left|\cfrac{u-u^{*}}{u_{n}^{*}}\right|\leq
O(\varepsilon t)$ on $I$ if $\lambda=0$.

Proof Multiply the integrating factor $e^{\lambda t+i\omega t}$ on both sides
of Eq. 1 and evaluate the integral on $(t_{0},t)\subseteq I$,

$\left|\int_{t_{0}}^{t}e^{\lambda\tau+i\omega\tau}\Big{(}u^{\prime}(\tau)+(\lambda+\omega
i)u(\tau)-f(\tau)\Big{)}\mathrm{d}\tau\right|\\\
\leq\int_{t_{0}}^{t}\left|e^{\lambda\tau+i\omega\tau}\Big{(}u^{\prime}(\tau)+(\lambda+\omega
i)u(\tau)-f(\tau)\Big{)}\right|\mathrm{d}\tau\leq\int_{t_{0}}^{t}\left|e^{\lambda\tau+i\omega\tau}\right|\varepsilon\mathrm{d}\tau$
(2)

The first part of inequality holds because modulus of integral is smaller than
integral of modulus. The second part holds by multiplying $e^{\lambda
t+i\omega t}$ on both sides of Eq. 1 and taking the integral on $(t_{0},t)$,
both of which preserve inequality property.

$\left|e^{\lambda t+i\omega t}u(t)-e^{\lambda t_{0}+i\omega
t_{0}}u(t_{0})-\int_{t_{0}}^{t}e^{\lambda\tau+i\omega\tau}f(\tau)\mathrm{d}\tau\right|\leq\varepsilon\int_{t_{0}}^{t}e^{\lambda\tau}\mathrm{d}\tau\\\
$ (3)

L.H.S. is reduced using
$\displaystyle{\int_{t_{0}}^{t}e^{\lambda\tau+i\omega\tau}\left(u^{\prime}+(\lambda+i\omega)u\right)\mathrm{d}\tau=\int_{t_{0}}^{t}\mathrm{d}\left(e^{\lambda\tau+i\omega\tau}u(\tau)\right)=\left[e^{\lambda\tau+i\omega\tau}u(\tau)\right]_{t_{0}}^{t}}$
and R.H.S. is reduced using $\left|e^{\lambda\tau+i\omega\tau}\right|\equiv
e^{\lambda\tau}$.

$\left|u(t)-e^{\lambda\left(t_{0}-t\right)+i\omega\left(t_{0}-t\right)}u(t_{0})-e^{-\lambda
t-i\omega
t}\int_{t_{0}}^{t}e^{\lambda\tau+i\omega\tau}f(\tau)\mathrm{d}\tau\right|\leq\varepsilon
e^{-\lambda t}\int_{t_{0}}^{t}e^{\lambda\tau}\mathrm{d}\tau$ (4)

Both sides are divided by $\left|e^{\lambda t+i\omega t}\right|$.

Notice that the analytical solution is given by

$u^{*}(t)=e^{\lambda(t_{0}-t)+i\omega(t_{0}-t)}u^{*}_{0}+e^{-\lambda t-i\omega
t}\int_{t_{0}}^{t}e^{\lambda\tau+i\omega\tau}f(\tau)\mathrm{d}\tau.$

Define the alternative solution to Eq. 1 under perturbed initial condition,
$u(t_{0})$ as

$\widetilde{u}(t):=e^{\lambda(t_{0}-t)+i\omega(t_{0}-t)}u(t_{0})+e^{-\lambda
t-i\omega t}\int_{t_{0}}^{t}e^{\lambda\tau+i\omega\tau}f(\tau)\mathrm{d}\tau.$

With this, Eq. 3 can be rewritten as

$\left|u(t)-\widetilde{u}(t)\right|\leq\varepsilon e^{-\lambda
t}\int_{t_{0}}^{t}e^{\lambda\tau}\mathrm{d}\tau.$ (5)

By the triangle inequality,

$\left|u(t)-u^{*}(t)\right|\leq\left|u(t)-\widetilde{u}(t)\right|+\left|\widetilde{u}(t)-u^{*}(t)\right|\leq\varepsilon
e^{-\lambda
t}\int_{t_{0}}^{t}e^{\lambda\tau}\mathrm{d}\tau+\left|\widetilde{u}(t)-u^{*}(t)\right|.$
(6)

As $\widetilde{u}=u^{*}$ when $u(t_{0})=u_{0}^{*}$, Eq. 6 is reduced to

$\left|u(t)-u^{*}(t)\right|\leq\varepsilon e^{-\lambda
t}\int_{t_{0}}^{t}e^{\lambda\tau}\mathrm{d}\tau.$ (7)

If $\lambda>0$, Eq. 7 gives rise to the absolute error bound

$\left|u(t)-u^{*}(t)\right|\leq\varepsilon\frac{1-e^{\lambda(t_{0}-t)}}{\lambda}\leq\frac{\varepsilon}{\lambda}=O(\varepsilon)\quad(\lambda>0).$
(8)

If $\lambda<0$, dividing Eq. 7 by
$\left|u_{n}^{*}(t)\right|=\left|e^{\lambda(t_{0}-t)+i\omega(t_{0}-t)}u_{0}^{*}\right|=e^{\lambda(t_{0}-t)}|u_{0}^{*}|$
yields the relative error bound

$\left|\frac{u(t)-u^{*}(t)}{u_{n}^{*}(t)}\right|\leq\varepsilon\frac{e^{-|\lambda|(t-t_{0})}-1}{-|\lambda||u_{0}^{*}|}=\varepsilon\frac{1-e^{-|\lambda|(t-t_{0})}}{|\lambda||u_{0}^{*}|}\leq\frac{\varepsilon}{|\lambda||u_{0}^{*}|}=O(\varepsilon)\quad(\lambda<0).$
(9)

If $\lambda=0$, the integral on R.H.S. of Eq. 7 is reduced to $(t-t_{0})$, and
therefore the absolute error bound is

$\left|u(t)-u^{*}(t)\right|\leq\varepsilon(t-t_{0})=O(\varepsilon
t)\quad(\lambda=0).$ (10)

Since the natural response has constant modulus
$\left|u_{n}^{*}(t)\right|=\left|e^{i\omega(t_{0}-t)}u_{0}^{*}\right|\equiv\left|u_{0}^{*}\right|$
when $\lambda=0$, the relative error with respect to the natural response is
bounded by $O(\varepsilon t)$ as well.

###  4.2 Higher-Order Linear ODE with Constant Coefficients

Proposition Let the residual $Ru(t)$ of the higher-order equation
$u^{(n)}+a_{n-1}u^{(n-1)}+\dots+a_{0}u=f$ be bounded by some $\varepsilon\geq
0$ on $I$, where $u^{(n)}$ is the $n$-th order derivative of $u$, namely,

$\left|u^{(n)}+a_{n-1}u^{(n-1)}+\dots+a_{0}u-f\right|\leq\varepsilon\quad\forall
t\in I.$ (11)

Let the network solution $u$ satisfy initial conditions
$u^{(k)}(t_{0})=u_{0}^{*(k)}$ ($k=0,\dots,n-1$). By the fundamental theorem of
algebra, the characteristic polynomial $p_{c}(x)$ can be uniquely factorized
as

$p_{c}(x):=x^{n}+a_{n-1}x^{n-1}+\dots+a_{0}=\prod_{k=0}^{n-1}(x+\lambda_{k}+i\omega_{k}).$
(12)

It is well-known that the exact solution has the form
$u^{*}(t)=u_{p}^{*}(t)\sum\limits_{k=0}^{n-1}c_{k}\exp(\lambda_{k}t+i\omega_{k}t)$,
where $u_{p}^{*}$ is any particular solution to the original equation and
$c_{0},\dots,c_{n-1}$ are constants chosen to satisfy the initial conditions.

Let $m$ be the total number of $k$ in Eq. 12 such that $\lambda_{k}=0$, then
the absolute error is bounded by:

$\left|u-u^{*}\right|\leq O(\varepsilon t^{m})\,\,\,{\rm
if}\,\,\,\lambda_{k}\geq 0\,\,\,{\rm for\,\,all}\,\,\,k.$ (13)

Proof For brevity, we prove the second-order case here to provide an intuition
of the complete proof, which is presented in Appendix B.

In the second-order case, Eq. 11 can be reduced to

$\left|u^{\prime\prime}+\left(\lambda_{1}+i\omega_{1}+\lambda_{2}+i\omega_{2}\right)u^{\prime}+\left(\lambda_{1}+i\omega_{1}\right)\left(\lambda_{2}+i\omega_{2}\right)u-f\right|\leq\varepsilon\quad(\lambda_{1}\geq\lambda_{2}),$
(14)

or, equivalently,

$\left|\Big{(}u^{\prime}+\left(\lambda_{1}+i\omega_{1}\right)u\Big{)}^{\prime}+\left(\lambda_{2}+i\omega_{2}\right)\Big{(}u^{\prime}+\left(\lambda_{1}+i\omega_{1}\right)u\Big{)}-f\right|\leq\varepsilon.$
(15)

Let $v=u^{\prime}+(\lambda_{1}+i\omega_{1})u$, Eq. 15 is then reduced to a
first-order inequality w.r.t. $v$

$\left|v^{\prime}+\left(\lambda_{2}+i\omega_{2}\right)v-f\right|\leq\varepsilon.$
(16)

By Eq. 7,

$\left|v(t)-v^{*}(t)\right|\leq\varepsilon
e^{-\lambda_{2}t}\int_{t_{0}}^{t}e^{\lambda_{2}\tau}\mathrm{d}\tau,$ (17)

where $v^{*}(t)={u^{*}}^{\prime}(t)+(\lambda_{1}+i\omega_{1})u^{*}(t)$.
Substituting $v=u^{\prime}+(\lambda_{1}+i\omega_{1})u$ into Eq. 17 yields

$\left|u^{\prime}(t)+\left(\lambda_{1}+i\omega_{1}\right)u(t)-v^{*}(t)\right|\leq\varepsilon
e^{-\lambda_{2}t}\int_{t_{0}}^{t}e^{\lambda_{2}\tau}\mathrm{d}\tau=\varepsilon\frac{1-e^{\lambda_{2}(t_{0}-t)}}{\lambda_{2}}$
(18)

Multiplying Eq. 18 by $e^{\lambda_{1}t+i\omega_{1}t}$, taking the integral on
$(t_{0},t)\subseteq I$, and dividing by
$\left|e^{\lambda_{1}t+i\omega_{1}t}\right|$, we have

$\left|u(t)-u^{*}(t)\right|\leq\varepsilon\frac{1}{\lambda_{1}\lambda_{2}}\left(1-\frac{\lambda_{1}e^{-\lambda_{2}t}-\lambda_{2}e^{-\lambda_{1}t}}{\lambda_{1}-\lambda_{2}}\right)=:\varepsilon\phi(t;\lambda_{1},\lambda_{2})$
(19)

If $\lambda_{1},\lambda_{2}>0$, it can be verified that
$\phi(t;\lambda_{1},\lambda_{2})$ is strictly increasing on $I$ and is bounded
by $\left[0,\cfrac{1}{\lambda_{1}\lambda_{2}}\right)$. Therefore

$\left|u(t)-u^{*}(t)\right|\leq\cfrac{\varepsilon}{\lambda_{1}\lambda_{2}}=O(\varepsilon)$
(20)

If $\lambda_{1}>\lambda_{2}=0$, taking the limit $\lambda_{2}\to 0$ in Eq. 19,
there is

$\left|u(t)-u^{*}(t)\right|\leq\lim_{\lambda_{2}\to
0}\varepsilon\phi(t;\lambda_{1},\lambda_{2})=\frac{\varepsilon}{\lambda_{1}^{2}}\left(e^{-\lambda_{1}t}+\lambda_{1}t-1\right)\leq\cfrac{\varepsilon
t}{\lambda_{1}}=O(\varepsilon t).$ (21)

If $\lambda_{1}=\lambda_{2}=0$, taking the double limit
$\lambda_{1},\lambda_{2}\to 0$ in Eq. 19, there is

$\left|u(t)-u^{*}(t)\right|\leq\lim_{\lambda_{1},\lambda_{2}\to
0}\varepsilon\phi(t;\lambda_{1},\lambda_{2})=\frac{\varepsilon
t^{2}}{2}=O(\varepsilon t^{2}).$ (22)

A detailed derivation of Eq. 21 and Eq. 22 can be found in Appendex A.

###  4.3 System of First-Order Linear ODEs with Constant Coefficients

Proposition Let the $p$-norm of the residual $\left\|R\mathbf{u}(t)\right\|$
of the linear system $\mathbf{u}^{\prime}+A\mathbf{u}=\mathbf{f}$
($\mathbf{u},\mathbf{f}\in\mathbb{C}^{n}$ and $A\in\mathbb{C}^{n\times n}$) be
bounded by some $\varepsilon\geq 0$ on $I$, namely,

$\left\|\mathbf{u}^{\prime}+A\mathbf{u}-\mathbf{f}\right\|\leq\varepsilon\quad\forall
t\in I,$ (23)

and the network solution satisfy the initial condition
$\mathbf{u}(t_{0})=\mathbf{u}^{*}_{0}$. Denote the Jordan canonical form of
$A$ as

$J=M^{-1}AM=\begin{pmatrix}J_{1}&&&\\\ &J_{2}&&\\\ &&\ddots&\\\
&&&J_{m}\end{pmatrix}\text{ where
}J_{k}=\begin{pmatrix}\lambda_{k}+i\omega_{k}&1&&\\\ &\ddots&\ddots&\\\
&&\lambda_{k}+i\omega_{k}&1\\\ &&&\lambda_{k}+i\omega_{k}\end{pmatrix}\quad
k=1,\dots,m$

(24)

where $M$ is composed of generalized eigenvectors and $J_{k}$ ($1\leq k\leq
m\leq n$) is a $n_{k}\times n_{k}$ Jordan block ($n_{1}+\dots+n_{m}=n$). Then,
the absolute error is bounded by $\|\mathbf{u}-\mathbf{u}^{*}\|\leq
O(\varepsilon)$ if $\lambda_{k}>0$ for all $k$.

Proof With the substitution $\mathbf{v}:=M^{-1}\mathbf{u}$,
$\mathbf{g}:=M^{-1}\mathbf{f}$, Eq. 23 can be transformed into

$\left\|\mathbf{v}^{\prime}+J\mathbf{v}-\mathbf{g}\right\|=\left\|M^{-1}\mathbf{u}^{\prime}+M^{-1}A\mathbf{u}-M^{-1}\mathbf{f}\right\|\leq\left\|M^{-1}\right\|\left\|\mathbf{u}^{\prime}+A\mathbf{u}-\mathbf{f}\right\|\leq\left\|M^{-1}\right\|\varepsilon$
(25)

where $\left\|M^{-1}\right\|$ is the induced $p$-norm of $M^{-1}$. Each entry
in $\left(\mathbf{v}^{\prime}+J\mathbf{v}-\mathbf{g}\right)$ must be no
greater than $\left\|M^{-1}\right\|\varepsilon$ in order for Eq. 25 to hold.
To bound the error for each Jordan chain, we first define two auxiliary
sequence of functions $\\{h_{k}\\}$ and $\\{H_{k}\\}$, which will be useful in
following derivations.

$h_{k}(t;\lambda):=\frac{1}{\lambda^{k}}\left(1-\sum_{j=0}^{k-1}\frac{\lambda^{j}(t-t_{0})^{j}}{j!}e^{\lambda(t_{0}-t)}\right)\quad\text{and}\quad
H_{k}(t;\lambda):=\sum_{j=1}^{k}h_{k}(t;\lambda).$ (26)

Notice the property that, if $\lambda>0$

$0\leq h_{k}(t;\lambda)<\frac{1}{\lambda^{k}}\quad 0\leq
H_{k}(t;\lambda)<\sum_{j=1}^{k}\frac{1}{\lambda^{j}}\quad\forall t\in I.$

Now, consider the first Jordan chain,

$\displaystyle\left|v^{\prime}_{1}+(\lambda_{1}+i\omega_{1})v_{1}+v_{2}-g_{1}\right|$
$\displaystyle\leq\left\|M^{-1}\right\|\varepsilon$ (27) $\displaystyle\vdots$
$\displaystyle\left|v^{\prime}_{n_{1}-1}+(\lambda_{1}+i\omega_{1})v_{n_{1}-1}+v_{n_{1}}-g_{n_{1}-1}\right|$
$\displaystyle\leq\left\|M^{-1}\right\|\varepsilon$ (28)
$\displaystyle\left|v^{\prime}_{n_{1}}+(\lambda_{1}+i\omega_{1})v_{n_{1}}-g_{n_{1}}\right|$
$\displaystyle\leq\left\|M^{-1}\right\|\varepsilon$ (29)

If $\lambda_{1}>0$, Eq. 29 implies (by section 4.1) the absolute error bound
on $v_{n_{1}}$

$\left|v_{n_{1}}-v^{*}_{n_{1}}\right|\leq\left\|M^{-1}\right\|\varepsilon\frac{1-e^{\lambda_{1}(t_{0}-t)}}{\lambda_{1}}=H_{1}(t;\lambda_{1})\left\|M^{-1}\right\|\varepsilon$
(30)

Plugging Eq. 28 and Eq. 30 into the following triangle inequality yields

$\displaystyle\left|v^{\prime}_{n_{1}-1}+(\lambda_{1}+i\omega_{1})v_{n_{1}-1}+v^{*}_{n_{1}}-g_{n_{1}-1}\right|$
$\displaystyle\leq\left|v^{\prime}_{n_{1}-1}+(\lambda_{1}+i\omega_{1})v_{n_{1}-1}+v_{n_{1}}-g_{n_{1}-1}\right|+\left|v^{*}_{n_{1}}-v_{n_{1}}\right|$
$\displaystyle\leq\left\|M^{-1}\right\|\varepsilon+H_{1}(t;\lambda_{1})\left\|M^{-1}\right\|\varepsilon$
(31)

Apply the integrating factor technique again, there is

$\left|v_{n_{1}-1}-v^{*}_{n_{1}-1}\right|\leq
H_{2}(t;\lambda)\left\|M^{-1}\right\|\varepsilon$ (32)

Repeating the above procedure, there is

$\left|v_{1}-v^{*}_{1}\right|\leq
H_{n_{1}}(t;\lambda_{1})\left\|M^{-1}\right\|\varepsilon,\quad\left|v_{2}-v^{*}_{2}\right|\leq
H_{n_{1}-1}(t;\lambda_{1})\left\|M^{-1}\right\|\varepsilon,\quad\dots\quad,\left|v_{n_{1}}-v^{*}_{n_{1}}\right|\leq
H_{1}(t;\lambda_{1})\left\|M^{-1}\right\|\varepsilon$

(33)

If $\lambda_{1}=0$, it can be proven (see Appendix C) that

$\left|v_{1}-v^{*}_{1}\right|\leq\left\|M^{-1}\right\|\varepsilon\sum\limits_{j=1}^{n_{1}}\frac{(t-t_{0})^{j}}{j!},\quad\left|v_{2}-v^{*}_{2}\right|\leq\left\|M^{-1}\right\|\varepsilon\sum\limits_{j=1}^{n_{1}-1}\frac{(t-t_{0})^{j}}{j!},\quad\dots\quad,\left|v_{n_{1}}-v^{*}_{n_{1}}\right|\leq\left\|M^{-1}\right\|\varepsilon(t-t_{0})$

(34)

Similarly, if $\lambda_{k}>0$ for the $k$-th Jordan chain, then

$\displaystyle\left|v_{n_{1}+\dots+n_{k-1}+1}-v^{*}_{n_{1}+\dots+n_{k-1}+1}\right|$
$\displaystyle\leq H_{n_{k}}(t;\lambda)\left\|M^{-1}\right\|\varepsilon$
$\displaystyle\left|v_{n_{1}+\dots+n_{k-1}+2}-v^{*}_{n_{1}+\dots+n_{k-1}+2}\right|$
$\displaystyle\leq H_{n_{k}-1}(t;\lambda)\left\|M^{-1}\right\|\varepsilon$
$\displaystyle\vdots$
$\displaystyle\left|v_{n_{1}+\dots+n_{k-1}+n_{k}}-v^{*}_{n_{1}+\dots+n_{k-1}+n_{k}}\right|$
$\displaystyle\leq H_{1}(t;\lambda)\left\|M^{-1}\right\|\varepsilon$

It can be shown that, if $\lambda_{k}>0$ for all $k$, then

$\left\|\mathbf{v}-\mathbf{v}^{*}\right\|\leq\sqrt[p]{n}\left(\max_{k}\sum_{j=1}^{n_{k}}\frac{1}{\lambda_{k}^{j}}\right)\left\|M^{-1}\right\|\varepsilon.$
(35)

Substituting $\mathbf{u}=M\mathbf{v}$ into Eq. 35, we have the absolute error
bound on $\mathbf{u}$,

$\left\|\mathbf{u}-\mathbf{u}^{*}\right\|=\left\|M\mathbf{v}-M\mathbf{v}^{*}\right\|\leq\|M\|\left\|\mathbf{v}-\mathbf{v}^{*}\right\|\leq\sqrt[p]{n}\left(\max_{k}\sum_{j=1}^{n_{k}}\frac{1}{\lambda_{k}^{j}}\right)\mathrm{cond}(M)\varepsilon=O(\varepsilon)$
(36)

where $\mathrm{cond}(M)=\left\|M\right\|\left\|M^{-1}\right\|$ is the
condition number of $M$. Note that the matrix of generalized eigenvectors,
$M$, can be replaced with $MD$ where $D\in\mathbb{C}^{n\times n}$ is a
diagonal matrix. The infimum of condition number under right multiplication

$\mathrm{cond^{R}}(M):=\inf_{D\text{ diagonal}}\mathrm{cond}(MD)=\inf_{D\text{
diagonal}}\|MD\|\left\|D^{-1}M^{-1}\right\|$

has been studied for induced $1$-norm, $2$-norm, and $\infty$-norm in [17],
[18], and [19].

###  4.4 First-Order Linear ODE with Nonconstant Coefficients

Proposition Let the residual $\left|Ru(t)\right|$ of
$u^{\prime}+\left(p(t)+iq(t)\right)u=f(t)$
($p,q:I\to\mathbb{R},f:I\to\mathbb{C}$) be bounded by some $\varepsilon\geq 0$
on $I$, namely,

$\left|u^{\prime}+\left(p(t)+iq(t)\right)u-f(t)\right|\leq\varepsilon\quad\forall
t\in(t_{0},\infty),$ (37)

and the network satisfy the initial condition $u(t_{0})=u^{*}_{0}$, then the
absolute error is bounded by

$\left|u-u^{*}\right|\leq O(\varepsilon t)$ (38)

if $p(t)\geq 0$ for sufficiently large $t$ on $I$.

Proof Denote the antiderivatives of $p(t)$ and $q(t)$ as

$P(t)=\int_{t_{0}}^{t}p(\tau)\mathrm{d}\tau\quad
Q(t)=\int_{t_{0}}^{t}q(\tau)\mathrm{d}\tau.$

Applying the integrating factor technique again, there is

$\left|\int_{t_{0}}^{t}e^{P(\tau)+iQ(\tau)}\Big{(}u^{\prime}(\tau)+\big{(}p(\tau)+iq(\tau)\big{)}u(\tau)-f(\tau)\Big{)}\mathrm{d}\tau\right|\\\
\leq\int_{t_{0}}^{t}\left|e^{P(\tau)+iQ(\tau)}\right|\left|u^{\prime}(\tau)+\big{(}p(\tau)+iq(\tau)\big{)}u(\tau)-f(\tau)\right|\mathrm{d}\tau$
$\displaystyle\left|e^{P(t)+iQ(t)}u(t)-u(t_{0})-\int_{t_{0}}^{t}e^{P(\tau)+iQ(\tau)}f(\tau)\mathrm{d}\tau\right|$
$\displaystyle\leq\varepsilon\int_{t_{0}}^{t}e^{P(\tau)}\mathrm{d}\tau$
$\displaystyle\left|u(t)-e^{-P(t)-iQ(t)}u^{*}_{0}-e^{-P(t)-iQ(t)}\int_{t_{0}}^{t}e^{P(\tau)+iQ(\tau)}f(\tau)\mathrm{d}\tau\right|$
$\displaystyle\leq\varepsilon
e^{-P(t)}\int_{t_{0}}^{t}e^{P(\tau)}\mathrm{d}\tau$
$\displaystyle\left|u(t)-u^{*}(t)\right|$ $\displaystyle\leq\varepsilon
e^{-P(t)}\int_{t_{0}}^{t}e^{P(\tau)}\mathrm{d}\tau.$ (39)

Rewriting the R.H.S. of of Eq. 39, there is

$\left|u(t)-u^{*}(t)\right|\leq\varepsilon
t\left(1+\cfrac{\phi(t)}{te^{P(t)}}\right),$ (40)

where

$\phi(t)=\int_{t_{0}}^{t}e^{P(\tau)}\mathrm{d}\tau-
te^{P(t)}=\int_{t_{0}}^{t}\left(e^{P(\tau)}-e^{P(t)}\right)\mathrm{d}\tau.$
(41)

Let $p(t)\geq 0$ for $t>t^{\prime}$. Subsequently, $P(t)$ is nondecreasing for
$t>t^{\prime}$. Therefore,

$\phi(t)=\int_{t_{0}}^{t^{\prime}}\left(e^{P(\tau)}-e^{P(t)}\right)\mathrm{d}\tau+\int_{t^{\prime}}^{t}\left(e^{P(\tau)}-e^{P(t)}\right)\mathrm{d}\tau\leq\int_{t_{0}}^{t^{\prime}}\left(e^{P(\tau)}-e^{P(t)}\right)\mathrm{d}\tau=\phi(t^{\prime})\quad
t>t^{\prime}.$ (42)

Consequently,

$\frac{\phi(t)}{te^{P(t)}}\leq\max_{\tau\in[t_{0},t^{\prime}]}\left[\frac{\phi(\tau)}{\tau
e^{P(\tau)}}\right]=:M,$ (43)

and finally,

$\left|u(t)-u^{*}(t)\right|\leq\varepsilon t\left(1+M\right)=O(\varepsilon
t).$ (44)

###  4.5 Dividing the Intervals for a Tightened Error Bound

In Sections 4.1 to 4.4, we only consider the global maximum residual norm
$\varepsilon$ on $I$. However, one can also partition $I$ into subintervals
$I=I_{1}\cup I_{2}\cup\dots$ and consider the local maximum residual norm
$\varepsilon_{k}$ on $I_{k}$. This leads to an even tighter error bound since
$\varepsilon_{k}\leq\varepsilon$ for all $k$.

For instance, in the case for first-order linear ODE with constant
coefficients, the bound in Eq. 8 becomes

$\left|u-u^{*}\right|\leq e^{-\lambda
t}\int_{t_{0}}^{t}e^{\lambda\tau}\left|Ru(\tau)\right|\mathrm{d}\tau$ (45)

as $\max\limits_{k}\rho(I_{k})\to 0$, where $\rho(I_{k})$ is the diameter of
interval $I_{k}$.

##  5 Experimental Results

For any ODE (or system of ODEs) discussed in Section 4, we are able to bound
the error of any network by simply evaluating its infinite norm (maximum
residual). This is true for any neural network solution, regardless of how
well it is trained or trained at all. In addition, the error bound can be
tightened by dividing the domain into subintervals and evaluating the infinite
norm on each one.

We run the following experiments with the NeuroDiffEq library [20], which
provides a convenient and flexible framework for training neural networks to
solve differential equations. Unless otherwise specified, we use an Adam
optimizer with a learning rate of $1.0\times 10^{-3}$ and
$(\beta_{1},\beta_{2})=(0.9,\,0.999)$ for training the networks. The neural
networks are simple fully-connected neural networks, with two 32-unit hidden
layers and $\tanh$ activation function. The loss function we use is the
$L_{2}$-norm of the ODE residuals at sampled points in the domain. The
solution we choose is the one from the epoch with the lowest validation loss.
We apply the following reparametrization to enforce the initial conditions
$\mathbf{u}(t_{0})=\mathbf{u}_{0}$ and, where required,
$\cfrac{\mathrm{d}}{\mathrm{d}t}\mathbf{u}(t_{0})=\mathbf{u}^{\prime}_{0}$:

$\displaystyle\mathbf{u}(t)$
$\displaystyle=\mathbf{u}_{0}+\left(1-e^{-(t-t_{0})}\right)\mathrm{ANN}(t),$
(46) $\displaystyle\text{or}\quad\mathbf{u}(t)$
$\displaystyle=\mathbf{u}_{0}+\left(t-t_{0}\right)\mathbf{u}^{\prime}_{0}+\left(1-e^{-\left(t-t_{0}\right)^{2}}\right)\mathrm{ANN}(t).$
(47)

###  5.1 First-Order Linear ODE with Constant Coefficients

We consider the real-valued ODE with constant coefficients, under the initial
condition $u(0)=2$,

$u^{\prime}+3u=f(t)\quad t\in I=[0,3]$

for $f(t)=3t^{2}+5t+4$, $f(t)=6\cos 3t$, $f(t)=4e^{t}$, and
$f(t)=-9\ln(1+t)-(1-t)^{-2}$, where the exact solutions are
$u(t)=e^{-3t}+t^{2}+t+1$, $u(t)=e^{-3t}+\sin 3t+\cos 3t$,
$u(t)=e^{-3t}+e^{t}$, and $u(t)=e^{-3t}-3\ln(t+1)+(1+t)^{-1}$. By Eq. 8, the
error bound for a single interval is given by $\varepsilon(1-e^{-3t})/3$ where
$\varepsilon$ is the largest absolute residual over the interval.

Here we adopt the method of manufactured solutions for ease of computing
absolute error. Namely, we pick specific forcing terms $f(t)$ so that the
exact solutions are known. We then compute the absolute errors of neural
networks using the known exact solutions.

For each epoch, we uniformly sample 1024 points from the domain $I=[0,3]$ for
training. After $100$ and $1000$ epoch of training, we evaluate the residual
of the ANN solution over $I$ and derive the corresponding error bounds using
the propositions in Sections 4.1 and 4.5. In Figure 1, we evaluate the error
bounds by linearly dividing the domain $I$ into $1$, $10$, $100$ sub-intervals
and the corresponding error bounds are tighter as more sub-intervals are used.
The absolute error computed using exact solution consistently fall within our
derived bound, which verifies our approach for evaluating error bound.

Figure 1: ODE Residuals and Absolute Error of ANN solution of Stable First
Order System Under Forcing. The absolute error is calculated using the exact
solution.

###  5.2 Higher-Order Linear ODE with Constant Coefficients

Here we consider two types of second-order differential equation,

$u^{\prime\prime}+u=f\quad\text{and}\quad u^{\prime\prime}+4u^{\prime}+3u=f$

where the solution space of the the associated homogeneous solution has basis
$\\{\sin t,\cos t\\}$ and $\\{e^{-t},e^{-3t}\\}$ respectively. By Eq. 22 and
19 the error bounds for a single interval are $\varepsilon t^{2}/2$ and
$\varepsilon\left(2+e^{-3t}-3e^{-t}\right)/6\leq\varepsilon/3$ respectively,
where $\varepsilon$ is the largest absolute residual over the interval.

We pick the forcing terms and initial conditions as described in Table 1.

Equation | Forcing $f(t)$ | $u(0)$ | $u^{\prime}(0)$ | Exact Solution $u(t)$
---|---|---|---|---
$u^{\prime\prime}+u=f$ | $2e^{t}$ | $2.0$ | $2.0$ | $\sin t+\cos t+e^{t}$
$u^{\prime\prime}+u=f$ | $t^{2}+t+3$ | $2.0$ | $2.0$ | $\sin t+\cos t+t^{2}+t+1$
$u^{\prime\prime}+u=f$ | $\ln(t+1)-(t+1)^{-2}$ | $1.0$ | $2.0$ | $\sin t+\cos t+\ln(t+1)$
$u^{\prime\prime}+u=f$ | $2\cos t^{2}+(1-4t^{2})\sin t^{2}$ | $1.0$ | $1.0$ | $\sin t+\cos t+\sin t^{2}$
$u^{\prime\prime}+3u^{\prime}+4u=f$ | $8e^{t}$ | $3.0$ | $-3.0$ | $e^{-t}+e^{-3t}+e^{t}$
$u^{\prime\prime}+3u^{\prime}+4u=f$ | $3t^{2}+11t+9$ | $3.0$ | $-3.0$ | $e^{-t}+e^{-3t}+t^{2}+t+1$
$u^{\prime\prime}+3u^{\prime}+4u=f$ | $3\ln(t+1)+4(t+1)^{-1}-(t+1)^{-2}$ | $2.0$ | $-3.0$ | $e^{-t}+e^{-3t}+\ln(t+1)$
$u^{\prime\prime}+3u^{\prime}+4u=f$ | $6\cos t-2\sin t$ | $3.0$ | $-3.0$ | $e^{-t}+e^{-3t}+\sin t+\cos t$

Table 1: Experiment Setup for Section 5.2

We train the network on $I=[0,3]$ for $100$ and $1000$ epochs. The ODE
residual and error bound with $n=1,10,100$ subintervals are plotted in Figures
2 and 3.

Figure 2: Residuals and Error Bounds of Harmonic Oscillator Under Various
Forcing

Figure 3: Residuals and Error Bounds of Stable 2nd Order Equation Under
Various Forcing

###  5.3 System of First-Order Linear ODEs with Constant Coefficients

We consider the following system of linear ODEs

$\mathbf{u}^{\prime}+MJM^{-1}\mathbf{u}=\mathbf{f}$

under the initial condition
$\mathbf{u}(t_{0})=M\begin{pmatrix}1&1&1&1&1&1\end{pmatrix}^{T}$ where
$J\in\mathbb{R}^{6\times 6}$ is the Jordan canonical form
$\begin{pmatrix}J_{1}\\\ &J_{2}\\\ &&J_{3}\end{pmatrix}$ with Jordan blocks
$J_{1}=\begin{pmatrix}4&1\\\ &4&1\\\ &&4\end{pmatrix}$,
$J_{2}=\begin{pmatrix}3&1\\\ &3\end{pmatrix}$, and $J_{3}=2$, and the forcing
is determined by

$\mathbf{f}(t)=M\begin{pmatrix}\cos t+4\sin
t+e^{t}-1&5e^{t}-4+t^{2}&4t^{2}+2t&3t^{3}+3t^{2}+e^{2t}-1&5e^{2t}-3&\frac{1}{t+1}+2\ln(t+1)\end{pmatrix}^{T}$

We randomly sample orthogonal $6\times 6$ matrices $M=M^{-T}$ to ensure
$\mathrm{cond}(M)=1$. By Eq. 33 and 36, the error bound is
$\varepsilon\sqrt{H^{2}_{3}(t;4)+H^{2}_{2}(t;4)+H^{2}_{1}(t;4)+H^{2}_{2}(t;3)+H^{2}_{1}(t;3)+H^{2}_{1}(t;2)}$
$\leq\sqrt{6}\varepsilon/2$ for a single interval where $\varepsilon$ is the
largest residual norm over the interval. The system is solved for $t\in
I=[0,3]$ for 1000 epochs with 1024 uniformly resampled points from the
expanded domain $[-1,4]$ at each epoch. We use networks with two 512-unit
(instead of 32-unit) hidden layers due to the coupled nature of the system.
However, it should be pointed out that, the error bound holds true regardless
of the network size or how well the network is trained. Again, we divide $I$
into $n=1,10,100$ subintervals for increasingly tighter bounds. Figure 4 shows
the system residual norm, network solution, as well as error bounds.

Figure 4: ODE Residuals and Absolute Error of ANN solution (System of Linear
ODEs)

###  5.4 First-Order Linear ODE with Nonconstant Coefficients

In this section, we consider linear ODE with time-dependent coefficients, with
$p(t)$, $f(t)$, initial condition and derived error bound for a single
interval according to Eq. 39 tabulated in Table 2.

$u^{\prime}+p(t)u=f(t)\quad t\in I=[0,3]$ (48)

Forcing $f(t)$ | Coefficient $p(t)$ | IC $u(0)$ | Exact Solution $u(t)$ | Bound
---|---|---|---|---
$(2t+1)(t+1)^{-1}\cos t-t\sin t$ | $(t+1)^{-1}$ | $1.0$ | $(t+1)^{-1}+t\cos t$ | $\varepsilon\frac{t^{2}/2+t}{t+1}$
$(t+1)^{2}\left(t^{2}+1\right)^{-1}e^{t}$ | $2t\left(t^{2}+1\right)^{-1}$ | $2.0$ | $(t^{2}+1)^{-1}+e^{t}$ | $\varepsilon\frac{t^{3}/3+t}{t^{2}+1}$
$2t+t^{2}\cos t\left(1+\sin t\right)^{-1}$ | $\cos t\left(1+\sin t\right)^{-1}$ | $1.0$ | $\left(1+\sin t\right)^{-1}+t^{2}$ | $\varepsilon\frac{t-\cos t}{\sin t+1}$
$\left(1+(t+2)\ln(t+1)\right)\left(t+1\right)^{-1}$ | $(t+2)(t+1)^{-1}$ | $1.0$ | $e^{-t}\left(t+1\right)^{-1}+\ln(t+1)$ | $\varepsilon\frac{e^{t}}{e^{t}+1}$

Table 2: Experiment Setup $u^{\prime}+p(t)u=f(t)$ for Section 5.3, where
$\varepsilon$ in derived bounds is the largest absolute residual over a single
interval

We train the network for 100 and 1000 epoch with 1024 points uniformly
resampled from $I=[0,3]$. According to Section 4.4, the absolute error bound
is $O(\varepsilon t)$. By evenly dividing $I$ into $n=1,10,100$ subintervals,
we obtain the error bounds in Figure 5.

Figure 5: ODE Residuals and Absolute Error of ANN solution (Nonconstant
Coefficients)

##  6 Conclusion

In this work, we have ascertained the link between ODE residuals and error
bound. We have proven, that for stable ODE systems discussed above, the bound
only depends on characteristics of the ODE (or system of ODEs) and the maximum
residual norm. There are efficient ways to evaluate the bound over the
interval, as we did in Section 5.

In our experiments, we have shown that while these bounds are sometimes too
loose using only the global maximum residual norm, they are usually
asymptotically bounded by some constant. One can further tighten the bound by
dividing the domain into smaller subintervals. In our experiments, the
subintervals are linearly divided, but one can also use an adaptive quadrature
for this task [21].

##  7 Future Work

In this work, we tie linear ODEs residuals to the absolute error and showed
the error can be bounded by a function of residuals. A subsequent research
topic is what strategy can be used to ensure a low residual, which in turn
guarantees a low absolute error.

As another future extension, our proposed method may be generalizable to
system of linear ODEs with time-dependent coefficients. It is also interesting
to study if this method is applicable to local linear approximation of
nonlinear ODEs. Furthermore, spatial or spatiotemporal PDEs with Dirichlet
boundary conditions is also worth exploring.

## References

  * [1] Maziar Raissi, Paris Perdikaris and George E Karniadakis  “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations”  In _Journal of Computational Physics_ 378 Elsevier, 2019, pp. 686–707
  * [2] Kurt Hornik, Maxwell Stinchcombe and Halbert White  “Multilayer feedforward networks are universal approximators”  In _Neural networks_ 2.5 Elsevier, 1989, pp. 359–366
  * [3] Isaac E Lagaris, Aristidis Likas and Dimitrios I Fotiadis  “Artificial neural networks for solving ordinary and partial differential equations”  In _IEEE transactions on neural networks_ 9.5 IEEE, 1998, pp. 987–1000
  * [4] Cedric Flamant, Pavlos Protopapas and David Sondak  “Solving Differential Equations Using Neural Network Solution Bundles”  In _arXiv preprint arXiv:2006.14372_ , 2020
  * [5] Shaan Desai et al.  “One-Shot Transfer Learning of Physics-Informed Neural Networks”, 2021
  * [6] Isaac E Lagaris, Aristidis C Likas and Dimitris G Papageorgiou  “Neural-network methods for boundary value problems with irregular boundaries”  In _IEEE Transactions on Neural Networks_ 11.5 IEEE, 2000, pp. 1041–1049
  * [7] Kevin Stanley McFall and James Robert Mahan  “Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions”  In _IEEE Transactions on Neural Networks_ 20.8 IEEE, 2009, pp. 1221–1233
  * [8] Pola Lydia Lagari, Lefteri H Tsoukalas, Salar Safarkhani and Isaac E Lagaris  “Systematic construction of neural forms for solving partial differential equations inside rectangular domains, subject to initial, boundary and interface conditions”  In _International Journal on Artificial Intelligence Tools_ 29.05 World Scientific, 2020, pp. 2050009
  * [9] N Sukumar and Ankit Srivastava  “Exact imposition of boundary conditions with distance functions in physics-informed deep neural networks”  In _arXiv preprint arXiv:2104.08426_ , 2021
  * [10] Marios Mattheakis, David Sondak, Akshunna S Dogra and Pavlos Protopapas  “Hamiltonian neural networks for solving differential equations”  In _arXiv preprint arXiv:2001.11107_ , 2020
  * [11] Bing Yu  “The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems”  In _arXiv preprint arXiv:1710.00211_ , 2017
  * [12] Justin Sirignano and Konstantinos Spiliopoulos  “DGM: A deep learning algorithm for solving partial differential equations”  In _Journal of computational physics_ 375 Elsevier, 2018, pp. 1339–1364
  * [13] Kshitij Parwani and Pavlos Protopapas  “Adversarial Sampling for Solving Differential Equations with Neural Networks”  In _arXiv preprint arXiv:2111.12024_ , 2021
  * [14] Olga Graf, Pablo Flores, Pavlos Protopapas and Karim Pichara  “Uncertainty Quantification in Neural Differential Equations”  In _arXiv preprint arXiv:2111.04207_ , 2021
  * [15] Mengwu Guo and Ehsan Haghighat  “An energy-based error bound of physics-informed neural network solutions in elasticity”  In _arXiv preprint arXiv:2010.09088_ , 2020
  * [16] Aditi Krishnapriyan et al.  “Characterizing possible failure modes in physics-informed neural networks”  In _Advances in Neural Information Processing Systems_ 34, 2021
  * [17] Friedrich L Bauer and Alston S Householder  “Absolute norms and characteristic roots”  In _Numerische Mathematik_ 3.1 Springer, 1961, pp. 241–246
  * [18] Friedrich L Bauer  “Optimally scaled matrices”  In _Numerische Mathematik_ 5.1 Springer, 1963, pp. 73–87
  * [19] Richard D Braatz and Manfred Morari  “Minimizing the Euclidean condition number”  In _SIAM Journal on Control and Optimization_ 32.6 SIAM, 1994, pp. 1763–1768
  * [20] Feiyu Chen et al.  “NeuroDiffEq: A Python package for solving differential equations with neural networks”  In _Journal of Open Source Software_ 5.46, 2020, pp. 1931
  * [21] William Marshall McKeeman  “Algorithm 145: Adaptive numerical integration by Simpson’s rule”  In _Communications of the ACM_ 5.12 ACM New York, NY, USA, 1962, pp. 604

de2021error

##  Appendix A Derivation of Eq. 21 and Eq. 22

Consider the case when $\lambda_{1}>0$, $\lambda_{2}\to 0$, we have the
following limit

$\displaystyle\left|u-u^{*}\right|$ $\displaystyle\leq\lim_{\lambda_{2}\to
0}\frac{\varepsilon}{\lambda_{1}\lambda_{2}}\left(1-\frac{\lambda_{1}e^{-\lambda_{2}}-\lambda_{2}e^{-\lambda_{1}}}{\lambda_{1}-\lambda_{2}}\right)$
$\displaystyle=\lim_{\lambda_{2}\to
0}\frac{\varepsilon}{\lambda_{1}\lambda_{2}}\frac{\lambda_{1}-\lambda_{2}-\left(\lambda_{1}e^{-\lambda_{2}}-\lambda_{2}e^{-\lambda_{1}}\right)}{\lambda_{1}-\lambda_{2}}$
$\displaystyle=\lim_{\lambda_{2}\to
0}\frac{\varepsilon}{\lambda_{1}\lambda_{2}}\frac{\lambda_{1}\left(1-e^{-\lambda_{2}t}\right)-\lambda_{2}\left(1-e^{-\lambda_{1}t}\right)}{\lambda_{1}-\lambda_{2}}$
$\displaystyle=\lim_{\lambda_{2}\to
0}\frac{\varepsilon}{\lambda_{1}\left(\lambda_{1}-\lambda_{2}\right)}\frac{\lambda_{1}\left(1-e^{-\lambda_{2}t}\right)-\lambda_{2}\left(1-e^{-\lambda_{1}t}\right)}{\lambda_{2}}$
$\displaystyle=\lim_{\lambda_{2}\to
0}\frac{\varepsilon}{\lambda_{1}\left(\lambda_{1}-\lambda_{2}\right)}\left(\frac{\lambda_{1}\left(1-e^{-\lambda_{2}t}\right)}{\lambda_{2}}-\left(1-e^{-\lambda_{1}t}\right)\right)$
$\displaystyle=\frac{\varepsilon}{\lambda_{1}^{2}}\lim_{\lambda_{2}\to
0}\left(\frac{\lambda_{1}\left(1-e^{-\lambda_{2}t}\right)}{\lambda_{2}}-\left(1-e^{-\lambda_{1}t}\right)\right)$
$\displaystyle=\frac{\varepsilon}{\lambda_{1}^{2}}\left(\lambda_{1}t-1+e^{-\lambda_{1}t}\right)$
(49)
$\displaystyle\leq\frac{\varepsilon}{\lambda_{1}^{2}}\left(\lambda_{1}t\right)=\frac{\varepsilon
t}{\lambda_{1}}$

If we take the limit $\lambda_{1}\to 0$ on top of $\lambda_{2}\to 0$, step 49
can be simplified using Taylor expansion,

$\displaystyle\left|u-u^{*}\right|$ $\displaystyle\leq\lim_{\lambda_{1}\to
0}\frac{\varepsilon}{\lambda_{1}^{2}}\left(\lambda_{1}t-1+e^{-\lambda_{1}t}\right)$
$\displaystyle=\lim_{\lambda_{1}\to
0}\frac{\varepsilon}{\lambda_{1}^{2}}\left(\lambda_{1}t-1+1-\lambda_{1}t+\cfrac{1}{2}\lambda_{1}^{2}t^{2}+t^{3}O\left(\lambda_{1}^{3}\right)\right)$
$\displaystyle=\cfrac{\varepsilon t^{2}}{2}$

##  Appendix B General Case Proof of Section 4.2

Define the following sequence of auxiliary functions
$\\{\phi_{n}\\}_{n=1}^{\infty}$ on $I$,

$\phi_{n}(t;\lambda_{1:n})=\frac{1}{\prod\limits_{j=1}^{n}\lambda_{j}}-\sum_{k=1}^{n}\frac{e^{-\lambda_{k}(t-t_{0})}}{\lambda_{k}\prod\limits_{j=1,j\neq
k}^{n}\left(\lambda_{j}-\lambda_{k}\right)},$

where $\lambda_{1:n}$ is a tuple
$\left(\lambda_{1},\lambda_{2},\dots,\lambda_{n}\right)$. Note that with
$\phi_{0}(t)=1$, it can be demonstrated that $\\{\phi_{n}\\}_{n=1}^{\infty}$
satisfies the recurrence relation

$\phi_{n+1}(t;\lambda_{1:n+1})=e^{-\lambda_{n+1}t}\int_{t_{0}}^{t}e^{\lambda_{n+1}\tau}\phi_{n}(\tau;\lambda_{1:n+1})\mathrm{d}\tau\quad\text{for
}n\geq 0.$ (50)

It can also be proven that $\phi_{n}(t;\lambda_{1:n})$ is monotonically
increasing on $I$ if $\lambda_{1},\dots,\lambda_{n}\geq 0$ because

$\frac{\mathrm{d}}{\mathrm{d}t}\phi_{n}(t,\lambda_{1:n})=\sum_{k=1}^{n}\frac{e^{-\lambda_{k}(t-t_{0})}}{\prod\limits_{j=1,j\neq
k}^{n}\left(\lambda_{j}-\lambda_{k}\right)}\geq 0.$ (51)

Also, if $\lambda_{1},\dots,\lambda_{n}>0$, there is
$\lim\limits_{t\to\infty}\phi_{n}(t,\lambda_{1:n})=\prod\limits_{j=1}^{n}\lambda_{j}^{-1}$.

Let $u_{0}(t):=u(t)$,
$u_{1}(t):=u^{\prime}_{0}(t)+(\lambda_{n}+i\omega_{n})u_{0}(t)$,
$u_{2}(t):=u^{\prime}_{1}(t)+(\lambda_{n-1}+i\omega_{n-1})u_{1}(t)$, …,
$u_{n-1}(t)=u^{\prime}_{n-2}(t)+(\lambda_{2}+i\omega_{2})u_{n-2}(t)$, Eq. 11
can be written as

$\left|u_{n-1}^{\prime}+(\lambda_{1}+i\omega_{1})u_{n-1}-f\right|\leq\varepsilon,$
(52)

which is a first-order inequality in terms of $u_{n-1}$ as discussed in
Section 4.1. By Eq. 7,

$\left|u_{n-1}-u^{*}_{n-1}\right|\leq\varepsilon
e^{-\lambda_{1}t}\int_{t_{0}}^{t}e^{\lambda_{1}\tau}\mathrm{d}\tau=\varepsilon\phi_{1}(t;\lambda_{1}).$
(53)

Substitute
$u_{n-1}(t)=u^{\prime}_{n-2}(t)+(\lambda_{2}+i\omega_{2})u_{n-2}(t)$ back into
Eq. 53, we have

$\left|u_{n-2}+(\lambda_{2}+i\omega_{2})-u^{*}_{n-1}\right|\leq\varepsilon\phi_{1}(t;\lambda_{1}),$

which is a first order inequality in terms of $u_{n-2}$. Applying the
integrating factor trick again, we have

$\left|u_{n-2}-u^{*}_{n-2}\right|\leq\varepsilon
e^{-\lambda_{2}t}\int_{t_{0}}^{t}e^{\lambda_{2}t}\phi_{1}(\tau,\lambda_{1})\mathrm{d}\tau=\varepsilon\phi_{2}(t;\lambda_{1:2}).$
(54)

Repeating the above process yields

$\left|u-u^{*}\right|=\left|u_{0}-u^{*}_{0}\right|\leq\varepsilon\phi_{n}(t;\lambda_{1:n})$
(55)

##  Appendix C Proof of Equation 34

Take the limit $\lambda\to 0$ in Eq. 26, and applying Taylor expansions where
necessary, we have

$\displaystyle h_{k}(t;0)$ $\displaystyle=\lim_{\lambda\to
0}\frac{1}{\lambda^{k}}\left(1-\sum_{j=0}^{k-1}\frac{\lambda^{j}(t-t_{0})^{j}}{j!}e^{\lambda(t_{0}-t)}\right)$
$\displaystyle=\lim_{\lambda\to
0}\cfrac{e^{\lambda(t_{0}-t)}}{\lambda^{k}}\left(e^{\lambda(t-t_{0})}-\sum_{j=0}^{k-1}\frac{\lambda^{j}(t-t_{0})^{j}}{j!}\right)$
$\displaystyle=\lim_{\lambda\to
0}\cfrac{e^{\lambda(t_{0}-t)}}{\lambda^{k}}\sum_{j=k}^{\infty}\frac{\lambda^{j}(t-t_{0})^{j}}{j!}$
$\displaystyle=\lim_{\lambda\to
0}\cfrac{1}{\lambda^{k}}\left(\sum_{l=0}^{\infty}\frac{\lambda^{l}(t_{0}-t)^{l}}{l!}\right)\left(\sum_{j=k}^{\infty}\frac{\lambda^{j}(t-t_{0})^{j}}{j!}\right)$

Notice the lowest order term w.r.t. $\lambda$ in
$\displaystyle{\left(\sum_{l=0}^{\infty}\frac{\lambda^{l}(t_{0}-t)^{l}}{l!}\right)\left(\sum_{j=k}^{\infty}\frac{\lambda^{j}(t-t_{0})^{j}}{j!}\right)}$
is $\lambda_{k}$, which is attained only when $l=0$ and $j=k$. The coefficient
for the $\lambda^{k}$ term is given by

$\cfrac{(t_{0}-t)^{0}}{0!}\cdot\cfrac{(t-t_{0})^{k}}{k!}=\cfrac{(t-t_{0})^{k}}{k!}.$

Consequently,

$\displaystyle h_{k}(t;0)$ $\displaystyle=\lim_{\lambda\to
0}\cfrac{1}{\lambda^{k}}\left(\cfrac{(t-t_{0})^{k}}{k!}\lambda^{k}+O\left(\lambda^{k+1}\right)\right)=\cfrac{(t-t_{0})^{k}}{k!}$
$\displaystyle H_{k}(t;0)$
$\displaystyle=\sum_{j=1}^{k}h_{k}(t;0)=\sum_{j=1}^{k}\cfrac{(t-t_{0})^{j}}{j!}$

Eq. 34 is attained by plugging the above equality into Eq. 33.

##  Appendix D Examples of dividing domain into subintervals

In Section 4.5, we show that the error bound on $I=[0,t]$ can be further
tightened by evaluating the maximum absolute residuals on a sequence of
subintervals $I_{i}=[t_{i-1},t_{i}]$. We apply this technique for the
experiments in Section 5.

###  D.1 Second Order Linear Equation with Constant Coefficients

Consider a second-order linear equation with constant coefficients (assuming
$\lambda_{1},\lambda_{2}\geq 0$,

$u^{\prime\prime}(t)+(\lambda_{1}+i\omega_{1}+\lambda_{2}+i\omega_{2})u^{\prime}(t)+(\lambda_{1}+i\omega_{1})(\lambda_{2}+i\omega_{2})u(t)=f(t)$
(56)

An approximated solution yielded by a neural network does not exactly satisfy
the Eq. 56, but instead incurs what we call a residual term $r(t)$

$u^{\prime\prime}(t)+(\lambda_{1}+i\omega_{1}+\lambda_{2}+i\omega_{2})u^{\prime}(t)+(\lambda_{1}+i\omega_{1})(\lambda_{2}+i\omega_{2})u(t)=f(t)+r(t)$
(57)

Solutions of Eq. 57 and 56 differ by

$\displaystyle\Delta(t)$
$\displaystyle=e^{-{\lambda_{1}}t}\int_{s=0}^{s=t}e^{{\lambda_{1}}s}e^{-{\lambda_{2}}s}\left(\int_{\tau=0}^{\tau=s}e^{{\lambda_{2}}\tau}r(\tau)\mathrm{d}\tau\right)\mathrm{d}s$
$\displaystyle=e^{-{\lambda_{1}}t}\int_{s=0}^{s=t}e^{({\lambda_{1}}-{\lambda_{2}})s}\left(\int_{\tau=0}^{\tau=s}e^{{\lambda_{2}}\tau}r(\tau)\mathrm{d}\tau\right)\mathrm{d}s$
$\displaystyle=e^{-{\lambda_{1}}t}\int_{s=0}^{s=t}\int_{\tau=0}^{\tau=s}e^{({\lambda_{1}}-{\lambda_{2}})s}e^{{\lambda_{2}}\tau}r(\tau)\mathrm{d}\tau\mathrm{d}s$

Therefore

$\displaystyle|\Delta(t)|$ $\displaystyle\leq
e^{-{\lambda_{1}}t}\int_{s=0}^{s=t}\int_{\tau=0}^{\tau=s}e^{({\lambda_{1}}-{\lambda_{2}})s}e^{{\lambda_{2}}\tau}|r(\tau)|\mathrm{d}\tau\mathrm{d}s$
$\displaystyle=e^{-{\lambda_{1}}t}\int_{\tau=0}^{\tau=t}\int_{s=\tau}^{s=t}e^{({\lambda_{1}}-{\lambda_{2}})s}e^{{\lambda_{2}}\tau}|r(\tau)|\mathrm{d}s\mathrm{d}\tau$
$\displaystyle=e^{-{\lambda_{1}}t}\int_{\tau=0}^{\tau=t}e^{{\lambda_{2}}\tau}|r(\tau)|\left(\int_{s=\tau}^{s=t}e^{({\lambda_{1}}-{\lambda_{2}})s}\mathrm{d}s\right)\mathrm{d}\tau$
$\displaystyle=e^{-{\lambda_{1}}t}\int_{\tau=0}^{\tau=t}e^{{\lambda_{2}}\tau}|r(\tau)|\frac{e^{({\lambda_{1}}-{\lambda_{2}})t}-e^{({\lambda_{1}}-{\lambda_{2}})\tau}}{{\lambda_{1}}-{\lambda_{2}}}\mathrm{d}\tau$
$\displaystyle=\int_{\tau=0}^{\tau=t}|r(\tau)|\frac{e^{{\lambda_{2}}(\tau-t)}-e^{{\lambda_{1}}(\tau-t)}}{{\lambda_{1}}-{\lambda_{2}}}\mathrm{d}\tau$

Notice that
$\cfrac{e^{{\lambda_{2}}(\tau-t)}-e^{{\lambda_{1}}(\tau-t)}}{{\lambda_{1}}-{\lambda_{2}}}\geq
0$ for $\tau<t$. Let $M(a,b)=\max\limits_{a\leq\tau\leq
b}\left|r(\tau)\right|$, we have

$\displaystyle|\Delta(t)|$
$\displaystyle\leq\sum_{i=1}^{n}M(t_{i-1},t_{i})\int_{\tau=t_{i-1}}^{\tau=t_{i}}\frac{e^{{\lambda_{2}}(\tau-t)}-e^{{\lambda_{1}}(\tau-t)}}{{\lambda_{1}}-{\lambda_{2}}}\mathrm{d}\tau$
(58) $\displaystyle\leq
M_{(0,t)}\int_{\tau=0}^{\tau=t}\frac{e^{{\lambda_{2}}(\tau-t)}-e^{{\lambda_{1}}(\tau-t)}}{{\lambda_{1}}-{\lambda_{2}}}\mathrm{d}\tau.$
(59)

where $0=t_{0}<t_{1}<\dots<t_{n}=t$.

Eq. 58 sheds light on how to evaluate the error bound by subdividing interval
$[0,t]$ into $n$ subintervals. Namely, we first evaluate the maximum absolute
residual $M(t_{i-1},t_{i})$ on $[t_{i-1},t_{i}]$ as well as the integral
$\displaystyle{\int_{\tau=0}^{\tau=t}\frac{e^{{\lambda_{2}}(\tau-t)}-e^{{\lambda_{1}}(\tau-t)}}{{\lambda_{1}}-{\lambda_{2}}}\mathrm{d}\tau}$,
which always has a closed-form expression depending $\lambda_{1}$ and
$\lambda_{2}$. The absolute error at any $t$ is then bounded by the sum of the
products. In particular, Eq. 59 is the special case where we do not divide
$[0,t]$ into subintervals ($n=1$), which is discussed in Section 4.2.

In the special case where $\max\limits_{1\leq i\leq n}(t_{i}-t_{i-1})\to 0$,
there is

$\int_{\tau=0}^{\tau=t}|r(\tau)|\frac{e^{{\lambda_{2}}(\tau-t)}-e^{{\lambda_{1}}(\tau-t)}}{{\lambda_{1}}-{\lambda_{2}}}\mathrm{d}\tau=\sum_{i=1}^{n}M(t_{i-1},t_{i})\int_{\tau=t_{i-1}}^{\tau=t_{i}}\frac{e^{{\lambda_{2}}(\tau-t)}-e^{{\lambda_{1}}(\tau-t)}}{{\lambda_{1}}-{\lambda_{2}}}\mathrm{d}\tau.$
(60)

###  D.2 System of ODEs

Consider a Jordan chain of length $3$ and eigenvalue $(\lambda+i\omega)$.

$\displaystyle u^{\prime}_{1}(t)+(\lambda+i\omega)u_{1}(t)+u_{2}(t)$
$\displaystyle=f_{1}(t)$ $\displaystyle
u^{\prime}_{2}(t)+(\lambda+i\omega)u_{2}(t)+u_{3}(t)$ $\displaystyle=f_{2}(t)$
$\displaystyle u^{\prime}_{3}(t)+(\lambda+i\omega)u_{3}(t)$
$\displaystyle=f_{3}(t)$

The approximated solution given by the neural network incurs residuals
$r_{1}(t),r_{2}(t),r_{3}(t)$, namely,

$\displaystyle u^{\prime}_{1}(t)+(\lambda+i\omega)u_{1}(t)+u_{2}(t)$
$\displaystyle=f_{1}(t)+r_{1}(t)$ (61) $\displaystyle
u^{\prime}_{2}(t)+(\lambda+i\omega)u_{2}(t)+u_{3}(t)$
$\displaystyle=f_{2}(t)+r_{2}(t)$ (62) $\displaystyle
u^{\prime}_{3}(t)+(\lambda+i\omega)u_{3}(t)$ $\displaystyle=f_{3}(t)+r_{3}(t)$
(63)

Eq. 63 implies that

$\displaystyle|u_{3}-u^{*}_{3}|\leq e^{-\lambda
t}\int_{\tau=0}^{\tau=t}e^{\lambda\tau}|r_{3}(\tau)|\mathrm{d}\tau$

By triangle inequality, Eq. 62 becomes

$\displaystyle\left|u^{\prime}_{2}+\lambda u_{2}+u^{*}_{3}\right|$
$\displaystyle\leq\left|u^{\prime}_{2}+\lambda
u_{2}+u_{3}\right|+|u_{3}-u^{*}_{3}|\leq|r_{2}(t)|+e^{-\lambda
t}\int_{\tau=0}^{\tau=t}e^{\lambda t}|r_{3}(\tau)|\mathrm{d}\tau$
$\displaystyle\left|u_{2}-u^{*}_{2}\right|$ $\displaystyle\leq e^{-\lambda
t}\int_{\tau=0}^{\tau=t}e^{\lambda\tau}|r_{2}(\tau)|\mathrm{d}\tau+e^{-\lambda
t}\int_{s=0}^{s=t}\left(\int_{\tau=0}^{\tau=s}e^{\lambda\tau}|r_{3}(\tau)|\mathrm{d}\tau\right)\mathrm{d}s$
$\displaystyle=e^{-\lambda
t}\int_{\tau=0}^{\tau=t}e^{\lambda\tau}|r_{2}(\tau)|\mathrm{d}\tau+e^{-\lambda
t}\int_{\tau=0}^{\tau=t}\left(\int_{s=\tau}^{s=t}e^{\lambda\tau}|r_{3}(\tau)|\mathrm{d}s\right)\mathrm{d}\tau$
$\displaystyle=e^{-\lambda
t}\int_{\tau=0}^{\tau=t}e^{\lambda\tau}|r_{2}(\tau)|\mathrm{d}\tau+e^{-\lambda
t}\int_{\tau=0}^{\tau=t}(t-\tau)e^{\lambda\tau}|r_{3}(\tau)|\mathrm{d}\tau$

Apply the same procedure for Eq. 61, there is

$\displaystyle\left|u^{\prime}_{1}+\lambda u_{1}+u^{*}_{2}\right|$
$\displaystyle\leq\left|u^{\prime}_{1}+\lambda
u_{1}+u_{2}\right|+|u_{2}-u^{*}_{2}|$ $\displaystyle\leq|r_{1}(t)|+e^{-\lambda
t}\int_{\tau=0}^{\tau=t}e^{\lambda\tau}|r_{2}(\tau)|\mathrm{d}\tau+e^{-\lambda
t}\int_{\tau=0}^{\tau=t}(t-\tau)e^{\lambda\tau}|r_{3}(\tau)|\mathrm{d}\tau$
$\displaystyle|u_{1}-u^{*}_{1}|$ $\displaystyle\leq e^{-\lambda
t}\int_{\tau=0}^{\tau=t}e^{\lambda\tau}|r_{1}(\tau)|\mathrm{d}\tau+e^{-\lambda
t}\int_{s=0}^{s=\tau}\left(\int_{\tau=0}^{\tau=s}e^{\lambda\tau}|r_{2}(\tau)|\mathrm{d}\tau\right)\mathrm{d}s$
$\displaystyle\quad+e^{-\lambda
t}\int_{s=0}^{s=\tau}\left(\int_{\tau=0}^{\tau=s}(s-\tau)e^{\lambda\tau}|r_{3}(\tau)|\mathrm{d}\tau\right)\mathrm{d}s$
$\displaystyle=e^{-\lambda
t}\int_{\tau=0}^{\tau=t}e^{\lambda\tau}|r_{1}(\tau)|\mathrm{d}\tau+e^{-\lambda
t}\int_{\tau=0}^{\tau=t}\left(\int_{s=\tau}^{s=t}e^{\lambda\tau}|r_{2}(\tau)|\mathrm{d}s\right)\mathrm{d}\tau$
$\displaystyle\quad+e^{-\lambda
t}\int_{\tau=0}^{\tau=t}\left(\int_{s=\tau}^{s=t}(s-\tau)e^{\lambda\tau}|r_{3}(\tau)|\mathrm{d}s\right)\mathrm{d}\tau$
$\displaystyle=e^{-\lambda
t}\int_{\tau=0}^{\tau=t}e^{\lambda\tau}|r_{1}(\tau)|\mathrm{d}\tau+e^{-\lambda
t}\int_{\tau=0}^{\tau=t}(t-\tau)e^{\lambda\tau}|r_{2}(\tau)|\mathrm{d}\tau+e^{-\lambda
t}\int_{\tau=0}^{\tau=t}\frac{(t-\tau)^{2}}{2}e^{\lambda\tau}|r_{3}(\tau)|\mathrm{d}\tau$
$\displaystyle=\int_{\tau=0}^{\tau=t}e^{\lambda(\tau-t)}|r_{1}(\tau)|\mathrm{d}\tau+\int_{\tau=0}^{\tau=t}(t-\tau)e^{\lambda(\tau-t)}|r_{2}(\tau)|\mathrm{d}\tau+\int_{\tau=0}^{\tau=t}\frac{(t-\tau)^{2}}{2}e^{\lambda(\tau-t)}|r_{3}(\tau)|\mathrm{d}\tau$

Note that, with $M_{k}(a,b)=\max\limits_{a\leq\tau\leq b}|r_{k}(\tau)|$
($k=1,2,3$) and $0=t_{0}\leq t_{1}\leq\dots\leq t_{n}=t$,

$\displaystyle
0\leq\int_{0}^{t}\frac{(t-\tau)^{k}}{k!}e^{\lambda(\tau-t)}|r(\tau)|\mathrm{d}\tau$
$\displaystyle\leq\sum_{i=1}^{n}M_{k}(t_{i-1},t_{i})\int_{t_{i-1}}^{t_{i}}\frac{(t-\tau)^{k}}{k!}e^{\lambda(\tau-t)}\mathrm{d}\tau$
(64) $\displaystyle\leq
M_{k}(0,t)\int_{0}^{t}\frac{(t-\tau)^{k}}{k!}e^{\lambda(\tau-t)}\mathrm{d}\tau$
(65)

Again, Eq. 64 shows one can evaluate the absolute error bound by dividing $n$
subintervals. For each interval, one evaluates the maximum residual as well as
the integral (which has a closed-form expression). Eq. 65 is the special case
as discussed in Section 4.3, where subintervals are not used ($n=1$).