Certified machine learning: Rigorous a posteriori error bounds for PDE defined PINNs
B. Hillebrecht acknowledges funding from the International Max Planck Research School for Intelligent Systems (IMPRS-IS). B. Unger acknowledges funding from the DFG under Germany's Excellence Strategy – EXC 2075 – 390740016 and The Ministry of Science, Research and the Arts Baden-Württemberg, 7542.2-9-47.10/36/2. Both authors are thankful for support by the Stuttgart Center for Simulation Science (SimTech).

1,2: Stuttgart Center for Simulation Science, University of Stuttgart, Stuttgart, Germany<EMAIL_ADDRESS>;

ORCID: 0000-0001-5361-0505
ORCID: 0000-0003-4272-1079

Birgit Hillebrecht1
Benjamin Unger2

Prediction error quantification in machine learning has been left out of most methodological investigations of neural networks, for both purely data-driven and physics-informed approaches. Beyond statistical investigations and generic results on the approximation capabilities of neural networks, we present a rigorous upper bound on the prediction error of physics-informed neural networks. This bound can be calculated without the knowledge of the true solution and only with a priori available information about the characteristics of the underlying dynamical system governed by a partial differential equation. We apply this a posteriori error bound exemplarily to four problems: the transport equation, the heat equation, the Navier-Stokes equation and the Klein-Gordon equation.

Physics-informed neural network, machine learning, certification, a posteriori error estimator, Navier-Stokes

§ INTRODUCTION

Physics-informed machine learning is applied to numerous highly complex problems, such as turbulence and climate modeling [44, 7, 45], model predictive control [3, 33], and Hamiltonian system dynamics [16]. The systematic study of physics-informed machine learning as a method, however, remains open [30]. Part of this is the question about the quality of the forecast: how close is the prediction of the physics-informed neural network (PINN) [36] to the actual solution?

Beyond statistical evaluations, such as in [22], first steps to answer this question were taken in [18] by establishing rigorous error bounds for PINNs approximating the solution of ordinary differential equations. We extend these results by deriving guaranteed upper bounds on the prediction error for problems modeled by linear partial differential equations (PDEs). Our methodology is thereby applicable to, e.g., the previously mentioned use cases of PINNs in computational science and engineering. Similar to [18], the error bound can be computed a posteriori without knowing the actual solution.

Since there are extensive studies on the convergence and approximation properties of (physics-informed) neural networks [4, 20, 11, 39], it is clear from a theoretical point of view that a suitable neural network (NN) can be found with the desired accuracy. These a priori results and the associated error estimates [39, 19, 8] are fundamentally different compared to the computable a posteriori error bounds we present here. Moreover, we like to stress that our analysis does not rely on statistical approaches or the Bayesian perspective but provides a rigorous error bound, which serves as a certificate for the NN.

Close in topic is [39] wherein a so called 'a posterior' error estimate is obtained based on norm inequalities. But, as stated in their conclusion, the estimates are not directly computable due to unknown constants. In distinction, the error estimator presented here is computable which we demonstrate with several academic examples. Other earlier works [31, 5] also differ from the following work in fundamental aspects, such as that the error bounds derived in them are either not guaranteed to hold, are restricted to a certain problem type, or require discretization steps.

At this point, we would like to emphasize that the results presented in <Ref> are independent of PINNs as the chosen methodology but can be applied to other surrogate modeling techniques. Nevertheless, since PINNs use the norm of the residual as a regularization term during the training process, they appear as a natural candidate to apply our residual-based error estimator.

We use bold notation for vectors and vector-valued functions and italic letters for operators, especially we use $\mathcal{I}$ to denote the identity operator. The notations $\dot{\mathbf{u}} = \partial_t \mathbf{u} = \frac{\partial \mathbf{u}}{\partial t}$ are used interchangeably to denote partial derivatives w.r.t. time. Similarily, we shorten derivatives w.r.t. spatial variables and use conventional notation for the divergence $\nabla \cdot \phi = \mathrm{div}(\phi)$, the gradient $\nabla \mathbf{u} = \mathrm{grad}(\mathbf{u})$, and the Laplace operator $\Delta \mathbf{u} = \mathrm{div}(\mathrm{grad}(\mathbf{u}))$ w.r.t. to spatial variables only.
For normed spaces $(X, \| \cdot \|_X)$, $(Y, \|\cdot \|_Y)$, and operator $\mathcal{A}\colon D(\mathcal{A}) \subseteq X \rightarrow Y$ we define the induced operator norm $\| \mathcal{A} \|_{X, Y} \coloneqq \mathrm{sup}_{x\in D(\mathcal{A}), \|x\|_X \neq 0} \frac{\|\mathcal{A}x\|_Y}{\|x\|_X}$. In case, $X \subseteq Y$ and $\| \cdot \|_Y = \| \cdot \|_X$, we shorten the notation to $\| \mathcal{A} \|_{X}$, and drop the subscript whenever the norm is clear from the context. For the spatial domain $\Omega \subseteq \R^d$, we denote the boundary by $\partial \Omega$ and define the size of $\Omega$ as $\|\Omega \| = \int_\Omega 1 \mathrm{d}\mathbf{x}$. The considerd temporal domain is denoted by $\timeInt = [0, t_f)$ with $t_f \in \R_+$.
We use the conventional notation for spaces of p-integrable functions $L^p(\Omega)$, and Sobolev spaces $W^{k,p}(\Omega)$ and $H^k(\Omega) = W^{k, 2}(\Omega)$.

§ PROBLEM DESCRIPTION

In the following, we aim at solving the boundary and inital value problem (BIVP): find $\mathbf{u}\colon \timeInt \times \Omega \rightarrow \R^n$ such that

\begin{equation}
\label{eq::pde}
\left\{\quad
\begin{aligned}
\dt \mathbf{u} &= \mathcal{A}\mathbf{u} \qquad &\mathrm{in } \quad  \timeInt \times \Omega , \\
\mathbf{u} &= \mathbf{u_0}  \qquad& \mathrm{in } \quad  \{t=0\} \times \Omega , \\
\mathcal{B} \mathbf{u} &= \mathbf{u_b} \qquad &\mathrm{in } \quad \timeInt \times \partial \Omega .
\end{aligned}\right.
\end{equation}

The initial value $\mathbf{u_0}$ lies in the domain of definition of the linear differential operator $\mathcal{A}\colon D(\mathcal{A}) \subseteq X \rightarrow X$, meaning $\mathbf{u_0} \in D(\mathcal{A})$.
The bounded linear operator $\mathcal{B}\colon D(\mathcal{A}) \rightarrow U$  represents the type of the boundary condition and maps to the Banach space $U$ which contains the boundary condition $\mathbf{u_b}$.
As an example, for the Dirichlet boundary condition $\mathcal{B} = \mathcal{T} \mathcal{I}$ is suitable, wherein $\mathcal{T}$ denotes the trace operator.

With this setup given, the objective can be intuitively formulated as follows.

Given an approximate solution $\hatu$ for the BIVP (<ref>), which is determined, for example, by approximating the system with a NN. Find a computable certificate $\varepsilon \colon \timeInt \rightarrow \R_+$ such that

\begin{equation*}
\| \hat{\mathbf{u}}(t, \cdot)- \mathbf{u}(t, \cdot) \|_X \le \varepsilon(t)
\end{equation*}

without knowing the true solution $\uu$.

In the following, we consider the BIVP (<ref>) in the framework of semigroups  and investigate mild solutions

\begin{equation*}
\mathbf{u}(t, \cdot) = \mathcal{S}(t) \mathbf{u_0}(\cdot)
\end{equation*}

defined by the semigroup of operators $\{\mathcal{S}(t)\}_{t\ge 0}$ generated by $\mathcal{A}$.
To assert that the problem is well-defined and that a unique (mild) solution exists, we make numerous assumptions.
Firstly, for $\mathcal{B}$ to be well-defined, we need to assert the existence of a trace operator $\mathcal{T}$ via trace theorems, e.g. requiring $\Omega$ to be a Lipschitz domain [9] or $\partial \Omega$ to be $C^1$ <cit.>.
By formulating the following assumption, we also include generic boundary conditions [37].

For a reflexive Banach space $X$, the linear (differential) operator $\mathcal{A}' = \mathcal{A}|_{\mathrm{ker} \mathcal{B}}$

generates a strongly continuous one-parameter semigroup of bounded linear operators $\{\mathcal{S}(t)\}_{t\ge 0}$ on $X$. The domain $D(\mathcal{A'}) =D(\mathcal{A}) \cap \mathrm{ker}(\mathcal{B})$ is a linear subspace of X. The boundary operator $\mathcal{B}$ is right invertible with the inverse $\mathcal{B}_0$ and $\mathcal{A}' \mathcal{B}_0$ is a bounded linear operator from $U$ to $X$.

A mild solution can then be defined by an abstract variation-of-constants type formula with an extended solution space. To make the presentation accessible to a large audience, we omit the general functional analytic framework here (see [37] and the references therein for the details). Furthermore, we drop the prime on the operator $\mathcal{A}'$ to simplify the notation.

The results derived below can be generalized directly to systems modeled by inhomogeneous PDEs of type

\begin{equation*}
\partial_t \mathbf{u} = \mathcal{A} \mathbf{u} + \mathbf{f}
\end{equation*}

with $\mathbf{f}(t) \in D(\mathcal{A})$ and $\mathbf{f}$ independent of $\mathbf{u}$. To simplify notation, we continue all investigations without $\mathbf{f}$.

Similarly to [18] for a finite-dimensional setting, the NN approximation of the solution of the BIVP (<ref>) is interpreted as the solution of a perturbed problem. Here we consider two types of perturbed problems, the first has the same boundary condition as the unperturbed problem, but has perturbation terms in the initial condition and the differential equation. The perturbed problem reads

\begin{equation}\label{eq::pde_approx}
\left\{\quad
\begin{aligned}
\dt \hatu &= \mathcal{A}\hatu + \mathbf{R} \qquad &\mathrm{in } \quad  \timeInt \times \Omega, \\
\hatu &= \mathbf{u_0} + \mathbf{R_0} \qquad& \mathrm{in } \quad  \{t=0\} \times \Omega ,\\
\mathcal{B}\hatu &= \mathbf{u_b} \qquad &\mathrm{in } \quad \timeInt \times \partial \Omega ,
\end{aligned}\right.
\end{equation}

with perturbances $\mathbf{R_0} \in D(\mathcal{A})$ and $\mathbf{R}:\; \timeInt \rightarrow X$.
For an approximate solution $\hatu$ given by a NN, we can determine the  terms $\mathbf{R}, \; \mathbf{R_0}$ by computing the residuals

\begin{equation} \label{eq::residuals}
\begin{aligned}
\mathbf{R_0} &\coloneqq \hat{\mathbf{u}}(0, \cdot) -\mathbf{u_0} \quad &\text{in } &\Omega , \\
\mathbf{R} &\coloneqq \dt \hatu - \mathcal{A}\hatu \quad &\text{in }& \timeInt \times \Omega .
\end{aligned}
\end{equation}

This calculation requires $\hatu$ to be differentiable in time, which constrains the choice of the activation function of the NN. Also, we assume the following for the BIVP (<ref>).

The perturbed BIVP (<ref>) is such that $\mathbf{R}$ is Lipschitz continuous in $t\in \timeInt$ and $\mathbf{u_0}+\mathbf{R_0} \in D(\mathcal{A})$.

Regarding the approximation by a NN, for most activation functions, the Lipschitz condition is satisfied, and only the second condition requires a more detailed investigation of the corresponding function spaces.

The second problem involves additionally a perturbation in the boundary condition

\begin{equation}\label{eq::pde_approx_soft_bcs}
\left\{\quad\begin{aligned}
\dt \hatu &= \mathcal{A}\hatu + \mathbf{R} \qquad &\mathrm{in } \quad  \timeInt \times \Omega ,\\
\hatu &= \mathbf{u_0} + \mathbf{R_0} \qquad& \mathrm{in } \quad  \{t=0\} \times \Omega ,\\
\mathcal{B}\hatu &= \mathbf{u_b} + \mathbf{R_b}\qquad &\mathrm{in } \quad \timeInt \times \partial \Omega .
\end{aligned}\right.
\end{equation}

We consider this system as a boundary control system (BCS, c.f. [41]) and use the notion of input-to-state-stable (ISS) systems, see [40], to later derive rigorous error bounds.

System (<ref>) is said to be ISS w.r.t. the boundary term $\mathbf{u_b}+\mathbf{R_b}$, if there exist functions $\beta, \gamma$ and an operator $\mathcal{C}$ such that

\begin{equation}\label{eq::iss}
\begin{aligned}
\| \mathbf{u}(t, \cdot) \| \le & \beta (\|\mathbf{u}(0, \cdot) \|, t) \\ &+ \gamma (\|\mathcal{C} (\mathbf{R_b}(s) + \mathbf{u_b}(s))\|_{L^\infty(0,t; \Omega)}) .
\end{aligned}
\end{equation}

Here, $\gamma\colon \R_+ \rightarrow \R_+$ is in $\mathcal{K}$, which means that it is continuous and strictly increasing with $\gamma(0)=0$. Similarily, $\beta\colon \R_+ \, \times\, \R_+ \rightarrow \R_+$ is in $\mathcal{KL}$, meaning $\beta(\cdot, t) \in \mathcal{K}$ for all $t\ge 0$ and $\beta(s, \cdot)$ is continuous and strictly decreasing to 0 for all $s>0$ [40].

In the following, we consider an ISS BCS, where we set $\mathbf{u_0} = 0$ and $\mathbf{u_b}=0$, such that (<ref>) simplifies to the second term. As an example, for the heat equation (see <Ref>) we use $\gamma(s) = 1/3 s$. For the imposed Dirichlet boundary, the operator $\mathcal{C} = \delta_{\partial\Omega}$ is the Dirac Distribution on the boundary of the domain $\Omega$.

We present reasoning and theorems based on the notion of ISS, however, all results can be transferred easily to integral input-to-state-stable (iISS) systems [40, 21].

The BCS (<ref>) is ISS with respect to the boundary values $\mathbf{u_b}+\mathbf{R_b}$.

With these definitions and assumptions given, we can clarify the intuitive understanding of <Ref> and present the following rigorous formulation.

Assume that the operators $\mathcal{A}$, $\mathcal{B}$ in the BIVP (<ref>) satisfy <Ref> and let $\uu$ denote the unique mild solution of (<ref>). Furthermore, let <Ref> hold true and let $\hatu$ denote the mild solution of (<ref>) or (<ref>). Find a certificate $\varepsilon(t)\colon \timeInt \rightarrow \R_+$ such that

\begin{equation*}
\| \hat{\mathbf{u}}(t, \cdot)- \mathbf{u}(t, \cdot) \|_X \le \varepsilon(t)
\end{equation*}

and can be computed a posteriori without knowledge of the true solution $\mathbf{u}$.

§ ERROR ESTIMATION

We divide this section, which contains the main results of this work, into three subsections, whereof the first two are concerned with the two previously introduced perturbed BIVPs. Here, we adapt the terminology of hard and soft boundary constraints in neural networks [46], meaning that hard boundary constraints denote the exact fulfilledness of the boundary conditions in BIVP (<ref>), and soft boundary constraints correspond to the perturbed BIVP (<ref>). For the networks used to generate these approximate solutions, we require a network to be designed to fulfill the boundary condition exactly for hard boundary constraints or force the network to approximate the boundary condition by a suitable loss term for soft constraints. In the third subsection, we give an auxiliary result to use numerical integration methods to compute the error bounds.

§.§ Hard boundary constraints

A NN approximation with hard boundary constraints can be interpreted as solution to the perturbed BIVP (<ref>) by defining the residuals as in (<ref>).

Let the BIVP (<ref>) be given and suppose that <Ref> is fulfilled. Furthermore, let $\hatu$ be  a mild solution to (<ref>) and presume that the terms $\mathbf{R}$ and $\mathbf{R_0}$ satisfy <Ref>. Then,

\begin{equation} \label{eq::elim}
\varepsilon(t) \coloneqq \| \mathcal{S}(t)\|_X  {\newdelta_0} + \int_0^t \|\mathcal{S}(t-s)\|_X {\newdelta}(s ) \,\mathrm{d}s ,
\end{equation}

wherein $\|\mathbf{R_0} \|_X \le \newdelta_0 \in \mathbb{R}_+$ and $\newdelta \colon \timeInt \rightarrow \mathbb{R}_+$, $\newdelta$ bounded and continuous with $\| \mathbf{R}(t) \|_X \le \newdelta(t)$ is a certificate for $\hatu$.

The dynamics of the error $\mathbf{e}(t,\mathbf{x}) \coloneqq  \hatu(t, \mathbf{x}) - \mathbf{u}(t, \mathbf{x})$ is governed by the BIVP

\begin{equation} \label{eq::pde_diff_system}
\begin{aligned}
\dt \mathbf{e} &= \mathcal{A}\mathbf{e}+ \mathbf{R} \qquad &\mathrm{in } \quad \timeInt \times \Omega,\\
\mathbf{e} &= \mathbf{R_0} \qquad& \mathrm{in } \quad \{t=0\} \times \Omega , \\
\mathcal{B}\mathbf{e} &= 0 \qquad &\mathrm{in } \quad \timeInt \times \partial \Omega ,
\end{aligned}
\end{equation}

with $\mathbf{R}$, $\mathbf{R_0}$ as defined in (<ref>).
Using <Ref>, the unique mild solution <cit.> of BIVP (<ref>) is given by

\begin{equation*}
\mathbf{e}(t,\cdot) = \mathcal{S}(t)\mathbf{R_0}(\cdot) + \int _0^t \mathcal{S}(t-s)\mathbf{R}(s, \cdot) \mathrm{d}s.
\end{equation*}

Controlling the residuals from above by $\newdelta, \; \newdelta_0$, one can limit $\|\mathbf{e}(t, \cdot)\|_X$ from above by $\varepsilon$ as defined in (<ref>).

Let all prerequisites for <Ref> be fulfilled. Then there exist constants $M\geq 1$, $\omega\in\R$ such that

\begin{equation} \label{eq::eps_lim}
\varepsilon(t) \le \tilde{\varepsilon}(t) \vcentcolon= M \mathrm{e}^{\omega t} {\newdelta_0} + \int_0^t M \mathrm{e}^{\omega (t-s)} {\newdelta}(s )\ds .
\end{equation}

If furthermore $\{S(t)\}_{t\ge 0}$ is a

* $\omega$-contraction semigroup, then $M=1$;

* exponentially decaying semigroup, then $\omega < 0$.

By applying general results on strongly continuous semigroups <cit.> to (<ref>), the existence of $M \ge 1$ and $\omega \ge 0$ for (<ref>) is shown.
An $\omega$-contraction semigroup is characterized by $M=1$. By definition of an exponentially decaying semigroup, there exist constants $M \ge 1$ and $\mu > 0$ such that $\|\mathcal{S}(t)\|_X \le M\mexp^{-\mu t}$.

Although <Ref> only establishes the existence of $M$ and $\omega$, we emphasize that these constants can be computed in several applications, see <Ref>. For example, we can use the spectral mapping theorem (e.g., <cit.>) to determine the growth bound $\omega$ to be the real part of the largest eigenvalue of the infinitesimal generator of the semigroup.

§.§ Soft boundary constraints

Enforcing boundary values with soft boundary constraints introduces an additional perturbation in the boundary, which is reflected in the perturbed BIVP (<ref>).

Let the BIVP (<ref>) fulfill <Ref> and assume that it is ISS according to <Ref>. Let the solution $\hatu$ to the perturbed problem (<ref>) satisfy <Ref>, then

\begin{equation*}
\begin{aligned}
\varepsilon(t) \coloneqq& \; M \mexp^{\omega t }{\newdelta_0} + \int_0^t M \mexp^{\omega (t -s)} {\newdelta}(s ) \,\mathrm{d}s  \\
&\qquad+\; \gamma ( \|\mathcal{C}\mathbf{R_b}\|_{L^\infty(0,t;\Omega)}),
\end{aligned}
\end{equation*}

with $\|\mathbf{R_0} \|_X \le \newdelta_0 \in \mathbb{R}_+$ and $\newdelta \colon \timeInt \rightarrow \mathbb{R}_+$, $\newdelta$ bounded and continuous with $\| \mathbf{R}(t) \|_X \le \newdelta(t)$ is a certificate for $\hatu$.

We introduce $\tildeu$ as a mild solution to (<ref>) with the same perturbations $\mathbf{R_0}, \mathbf{R}(t)$ considered in the BIVP (<ref>). Since a mild solution $\uu$ to (<ref>) exists, and <Ref> is valid, there exists a unique mild solution to (<ref>). Hence, the error can be split into

\begin{equation*}
\| \hatu - \uu \|_X \le \| \hatu - \tildeu \|_X + \| \tildeu - \uu\|_X .
\end{equation*}

According to <Ref> we can bound the latter part by (<ref>).
For the first contribution, we proceed similarly to the proof of <Ref> by observing that the dynamics of the error $\tildee(t,\mathbf{x}) \coloneqq \hatu - \tildeu $ is determined by the BIVP

\begin{equation*}
\begin{aligned}
\dt \tildee &= \mathcal{A}\tildee \qquad &\mathrm{in } \quad \timeInt \times \Omega,\\
\tildee &= 0 \qquad& \mathrm{in } \quad \{t=0\} \times \Omega  ,\\
\mathcal{B}\tildee &= \mathbf{R_b}\qquad &\mathrm{in } \quad \timeInt \times \partial \Omega .
\end{aligned}
\end{equation*}

Now, we use the ISS property according to <Ref> and notice that $\tildee(0, \cdot) = 0$ implies $ \beta(\| \tildee(0, \cdot) \|, t) =0$ (with $\beta$ as in <Ref>) and hence conclude

\begin{equation*}
\| \tildee(t, \cdot) \|_X \le \gamma (\|\mathcal{C}\mathbf{R_b} \|_{L^\infty(0, t; \Omega)}).
\end{equation*}

Adding the contributions to the error concludes the proof.

§.§ Integral approximation

As in [18], we give an auxiliary result to apply the previous theorems using numerical integration techniques to compute the integrals in <Ref>, <Ref> and <Ref>, respectively. Here, we use the composite trapezoidal rule to approximate the integral as

\begin{equation*}
\int_0^t M\mexp^{\omega (t-s)} \newdelta(s) \ds \approx M \mexp^{\omega t} \hat{I}_n(0,t,\mathrm{e}^{-\omega s}\newdelta(s))
\end{equation*}

using $n$ equally spaced subintervals $[s_i, s_{i+1}]$, $0 = s_1< ...< s_n= t$, meaning $s_{i+1}-s_{i} = \Delta s$. Therein the right hand side is defined via

\begin{equation}
\begin{aligned} \label{eq::trapz}
\hat{I}_n&(t, \mexp^{-\omega s } \newdelta(s) ) \coloneqq \\
&\sum_{i=1}^{n-1} \frac{\Delta s}{2} \left( \mexp^{-\omega s_i} \newdelta(s_{i}) + \mexp^{-\omega s_{i+1}} \newdelta(s_{i+1}) \right).
\end{aligned}
\end{equation}

Numerical computation of the integral introduces an additional error

\begin{equation} \label{eq::trapz_err}
\Eint \coloneqq M \mexp^{\omega t} \frac{K t^3}{12 n^2},
\end{equation}

assuming that $\newdelta(s)$ is sufficiently smooth and that $K$ fulfills

\begin{equation} \label{eq::K}
\left|\tfrac{\mathrm{d}^2}{\mathrm{d} s^2} (\mexp^{-\omega s} \newdelta(s))\right| \le K .
\end{equation}

Summarizing the previous discussion yields the following result.

Consider the BIVPs (<ref>) and (<ref>) (or (<ref>)) and suppose that <Ref> (and <ref>) are satisfied. Then

\begin{equation*}
\varepsilon(t) \le \;  M \mexp^{\omega t} \newdelta_0 + M \mexp^{\omega t}\hat{I}_{n}(t, M \mexp^{-\omega s}  \newdelta(s)) + \Eint + \Ebc
\end{equation*}

with $\newdelta_0, \,\,\newdelta$ as in <Ref> under the additional assumption that $\newdelta\in \mathrm{C}^3(\timeInt, \mathbb{R})$. Furthermore the numerical integral (<ref>) and the error of the numerical integration (<ref>) is used assuming $K$ can be found according to (<ref>). When considering the BIVP (<ref>), $\Ebc = \gamma (\|\mathcal{C}\mathbf{R_b} \|_{L^\infty(0, t; \Omega)})$, else the term vanishes.

The additional integration error considered in <Ref>lemma::errorlimit_trpz accounts for numerical integration over time only. The numerical error introduced by integrating over space is already included and part of $\newdelta(t)$ and $\newdelta_0$.

The previous derivations also show that an additional regularization of the derivatives (c.f. [46]) of the residuum could prove useful to reduce the integration error.

§ METHODOLOGY: PINNS

We want to apply the results from the previous <Ref> to PINNs as introduced in [36]. Therein the the data-driven loss used during training of feed-forward NNs is supplemented by a contribution reflecting the PDE (the physics-informed contribution), yielding

\begin{equation} \label{eq::loss}
L = \tfrac{1}{1+\kappa} L_\mathrm{eq} +\tfrac{\kappa}{1+\kappa} L_\mathrm{data}  + \rho L_\mathrm{space} .
\end{equation}

The parameter $\kappa \in \R_+$ defines the relation between data-driven and physics-informed contributions, while $\rho\geq 0$ scales the contribution of soft boundary conditions and further spatial properties, see below for further details. The data-driven contribution is encoded in

\begin{equation*}
L_\mathrm{data} = \frac{1}{N_\mathrm{data}} \sum_{i=1}^{N_\mathrm{data}} \|  \hatu(t_{\mathrm{data}, i}, \mathbf{x}_{\mathrm{data}, i}) - \uu_{\mathrm{data}, i}\|^2
\end{equation*}

for a given training dataset

\begin{equation*} (t_{\mathrm{data}, i}, \mathbf{x}_{\mathrm{data}, i}, \uu_{\mathrm{data}, i})_{i=1,...,N_\mathrm{data}}.
\end{equation*}

We need to know the true solution  $\uu_{\mathrm{data}, i} = \uu(t_{\mathrm{data}, i}, \mathbf{x}_{\mathrm{data}, i})$ for this and hence we reduce data-driven contributions to data given by the initial condition

\begin{equation*}
L_\mathrm{data} = \frac{1}{N_\mathrm{data}} \sum_{i=1}^{N_\mathrm{data}} \|  \hatu(0, \mathbf{x}_{\mathrm{data}, i}) - \mathbf{u_0}( \mathbf{x}_{\mathrm{data}, i})\|^2 .
\end{equation*}

The physics-informed loss

\begin{equation*}
L_\mathrm{eq} = \frac{1}{N_\mathrm{eq}} \sum_{i=1}^{N_\mathrm{eq}} \| \partial_t \hatu(t_{\mathrm{eq} , i}, \mathbf{x}_{\mathrm{eq} , i}) - \mathcal{A}\hatu(t_{\mathrm{eq} , i}, \mathbf{x}_{\mathrm{eq} , i})\|^2
\end{equation*}

evaluates the residual of the PDE for a set of $N_\mathrm{eq}$ collocation points $(t_{\mathrm{eq} , i}, \mathbf{x}_{\mathrm{eq} , i})_{i=1,..., N_\mathrm{eq}}$. In particular, this contribution to the loss term aims at minimizing $\mathbf{R}$, which is used in the error bounds derived in <ref>.

The third contribution $L_\mathrm{space}$ accounts for space constraints such as the divergence-freeness of solutions to the Navier-Stokes equations or soft boundary constraints. The parameter $\rho$ scales this loss relative to the other two contributions.

We use NNs with $n_\mathrm{i}\in \mathbb{N}$ input neurons, $n_\mathrm{o}\in \mathbb{N}$ output neurons and $n_\mathrm{hl}\in \mathbb{N}$ hidden layers with $n_\mathrm{hn}\in \mathbb{N}$ neurons each as basis for networks used in the numerical experiments. In the following two subsections, we discuss the modifications in network topology or loss functions which are necessary to realize hard and soft boundary constraints.

In case no further control parameters are given, $n_\mathrm{i}$ corresponds to the dimension of $\timeInt \times \Omega$, i.e., $n_\mathrm{i} = d+1$. Similarly, assuming that no hidden variables are included in the PDE which are not described explicitly in the PDE (like the pressure in the Navier-Stokes equations, c.f. <Ref>), $n_\mathrm{o} $ corresponds to the dimension of $\uu$, i.e., $n_\mathrm{o} = n$.

§.§ Hard boundary constraints

The implementation of hard boundary constraints requires a modification of the network topology. For the following numerical experiments, we need topologies to account for Dirichlet or periodic boundary conditions. For Dirichlet boundary conditions, we extend the purely sequential core NN by a second path as depicted in <Ref>. In the homogeneous case, we use a mask (as presented for arbitrary boundary conditions in [28, 29]  and therein called length factor) for a d-dimensional rectangular domain of width $w = (w_0,..., w_d)$ with center $x_0 = (x_{1,0}, ..., x_{d, 0})$ of the type

\begin{equation*}
\lambda(\mathbf{x}) = \prod_{i=1}^d \left[ (\tfrac{w_i}{2})^2 - (x_i - x_{i,0})^2 \right].
\end{equation*}

For inhomogeneous problems, we add a suitable offset

\begin{equation*}
\boldsymbol{\chi}(t, \mathbf{x}) = \mathbf{u_b}(t, \mathbf{x}) \qquad \mathrm{for }\;\; t\in\timeInt,\quad \mathbf{x}\in \partial \Omega.
\end{equation*}

We can summarize the model to be

\begin{equation} \label{eq::hard_bcs}
\hatu(t, x)  = \lambda(\mathbf{x}) \cdot \mathbf{\hat{\hat{u}}}(t, \mathbf{x}) + \boldsymbol{\chi}(t, \mathbf{x})
\end{equation}

wherein $\mathbf{\hat{\hat{u}}}$ is the output of the sequential core NN.

For periodic boundary condition, we consider again a d-dimensional rectangular domain of width $w = (w_0,..., w_d)$ and follow [10] to construct a periodic base layer for one (the $i$-th) space dimension with $n_\mathrm{p}$ units as

\begin{equation} \label{eq::per_base}
\sigma ( \mathbf{A} \cdot \cos(\boldsymbol{\nu} x + \boldsymbol{\phi}) + \mathbf{B}),
\end{equation}

wherein $\sigma$ is the activation function. Furthermore, $\mathbf{A}, \mathbf{B}, \boldsymbol{\nu}, \boldsymbol{\phi} \in \R^{n_p}$. The frequency $\nu$ may not be trained but needs to be fixed by the period $ w_i$ in the respective space direction as $\boldsymbol{\nu} =2 \pi / w_i  \cdot \mathbf{I}_{n_p}$. Here, $\mathbf{I}_{n_p}$ denotes the vector of all ones in $\R^{n_p}$.

The non-periodic input variables are processed by a dense layer.
Exemplarily, a non-sequential network realizing this is shown in <Ref>. To enhance training performance, we introduce a small dense network before the periodic layers to model the dependence of $\mathbf{A}^i, \mathbf{B}^i, \boldsymbol{\phi}^i$ on $t$ as an extension to the suggested modelling in [10].

[ node distance = 2cm, >=latex,block/.style=draw, fill=white, rectangle, minimum height=3em, minimum width=8em]
[block] (Input) Input;
[node distance=10mm]
[block, below =of Input] (NN)Hidden layers;
[block, below =of NN] (Output)NN Output;

[block, right =of NN] (LF) Mask;

[node distance=28mm]
[block, below =of LF] (RO) $\times$;

[block, right =of LF] (OFS) Offset;

[node distance=46mm]
[block, below =of OFS] (ROO) $+$;

[<-] (Input) – ++ (0,1.5) node[right] $\mathbf{x}$;
[<-, dashed] (Input.120) – ++ (0,1) node[left] $t$;

[->] (Input) – (NN);
[->, dashed] (Input.240) – (NN.120);
[->, dash dot dot] (NN) – ++ (0,-1)  -| node[pos=0.25]  (Output);
[->,dash dot dot] (Output) – ++ (0,-0.8)  -| node[pos=0.05, below] $\mathbf{\hat{\hat{u}}}(t,\mathbf{x})$ (RO.140);

[->] (Input.320) – ++ (0,-0.5)  -| node[pos=0.25, below] $\mathbf{x}$ (LF);
[->] (Input.320) – ++ (0,-0.5)  -| node[above]  (OFS);
[->, dashed] (Input.330) – ++ (0,-0.2)  -| node[pos=0.2, above] $t$  (OFS.50);
[->, dash dot dot] (LF) – ++ (0,-1)  -| node[pos=0.7, right] $\lambda(\mathbf{x})$ (RO);
[->,dash dot dot] (OFS) – ++ (0,-1)  -| node[pos=0.7, right] $\boldsymbol{\chi}(t,\mathbf{x})$ (ROO);

[->,dash dot dot] (RO) – ++ (0,-0.8) -| node[pos=0.2, below] $\lambda(\mathbf{x}) \cdot \mathbf{\hat{\hat{u}}}(t, \mathbf{x})$ (ROO.140);

[->,dash dot dot] (ROO) – ++ (0,-1) node[left] $\hatu(t, \mathbf{x}) = \lambda(\mathbf{x}) \cdot \mathbf{\hat{\hat{u}}}(t, \mathbf{x}) + \boldsymbol{\chi}(t, \mathbf{x})$;

NN realizing hard Dirichlet boundary constraints $\mathbf{u_b}(t, \mathbf{x})$ for $\mathbf{x} \in \partial \Omega$, $t\in \timeInt$. Input parameters are realized as dashed (time) and continuous (space) lines, processed intermediate or output values are shown as dash dotted lines.

[node distance=10mm and 4mm,>=latex,block/.style=draw, fill=white, rectangle, minimum height=3em, minimum width=6em]

[block] (Input) Input;
[<-] (Input.40) – ++ (0,1) node[left] $t$;
[<-, dashed] (Input.90) – ++ (0,1) node[left] $x_1$;
[<-, dotted] (Input.140) – ++ (0,1) node[left] $x_2$;

[block, below =of Input] (A)$\mathbf{A}$;
[block, left =of A] (B)$\mathbf{B}$;
[block, right =of A] (phi)$\boldsymbol{\phi}$;
[block, below =of B] (per)Periodic base layer;
[draw, black!50 , dashed,fit=(A) (B) (phi) (per)] (a12) ;
[right =of  per, black!50, xshift=0.5cm, align=left] (per_all) Enforces periodic
boundary condition;
[->] (Input) – ++ (0,-1)  -| node[pos=0.25]  (phi);
[->] (Input) – ++ (0,-1)  -| node[pos=0.25]  (A);
[->] (Input) – ++ (0,-1)  -| node[pos=0.25, below] $t$ (B);
[->, dotted] (Input) – ++ (-1,0)  -| node[above] $x_2$ (per.155);
[->, dash dot dot] (A) – ++ (0,-1)  -| node[pos=0.25]  (per.50);
[->, dash dot dot] (B) – ++ (0,-1)  -| node[pos=0.25]  (per.50);
[->, dash dot dot] (phi) – ++ (0,-1)  -| node[pos=0.25]  (per.50);

[block, below right =of phi, xshift=1cm] (lin)Dense layer;
[->] (Input) – ++ (1,0)  -| node[pos=0.25, above] $t$ (lin);
[->, dashed] (Input.345) – ++ (1,0)  -| node[pos=0.25, below] $x_1$ (lin.120);

[block, below right =of per] (concat)Hidden layers;

[->, dash dot dot] (concat) – ++ (0,-1) node[left] ...;
[->, dash dot dot] (per) – ++ (0,-1)  -| node[pos=0.25]  (concat);
[->, dash dot dot] (lin) – ++ (0,-1)  -| node[pos=0.25]  (concat);

Extended realization of a non sequential network to enforce the periodic boundary condition in the space variable $x_2$ according to (<ref>).
The periodic base layer is dependent on three parameters $\mathbf{B}, \mathbf{A}$ and $\phi$, which have an explicit dependency on the time $t$. These parameters and the input variable for which the output shall be periodic, $x_2$, serve then as input to the periodic base layer.
The other input variables $x_1$ and $t$ are bypassed around the additional elements through a dense layer.

§.§ Soft boundary constraints and other restrictions

In case, the solution shall meet additional constraints which are not satisfied automatically by the NN, e.g. divergence-freeness, we train the NN by adding $L_\mathrm{space}$ and weighting it with the parameter $\rho$.
Details for the choice of $\rho$ and the implementation of $L_\mathrm{space}$ can be found in the corresponding examples in <Ref>subsec::NSE, <ref> and <ref>.

If soft boundary constraints are used, the boundary error needs to be taken into account according to <Ref>.

§ NUMERICAL RESULTS

In this section, we investigate four examples and illustrate the applicability of the a posteriori error estimator for PDE-defined PINNs. The four examples are the heat equation, the transport equation, the Navier-Stokes equation, and the Klein-Gordon equation.

The code and data used to generate the subsequent results are accessible via doi: 10.5281/zenodo.7156168 under MIT Common License.

Herein, we apply <Ref> and illustrate the different contributions by splitting the predicted error in a part reflecting the propagation of the error in the initial condition $\Einit$ and the integrated error from the misrepresentation of the temporal evolution $\Eeq$ according to

\begin{equation*}
\begin{aligned}
\Einit (t) \coloneqq &\; M \mathrm{e}^{\omega t} \newdelta_0, \\
\Eeq (t) \coloneqq & \;\int_0^t M \mathrm{e}^{\omega (t-s)}  \newdelta(s )\ds .
\end{aligned}
\end{equation*}

We omit the numerical approximation steps as derived in <Ref> for notation simplicity. Furthermore, if we apply <Ref>, we denote the error introduced by the boundary approximation as

\begin{equation*}
\Ebc (t) \coloneqq \;\gamma(\|\mathcal{C}\mathbf{R_b}\|_{L^\infty(0,t; \Omega)} ).
\end{equation*}

The reference error is denoted by

\begin{equation*}
\Eref (t) \coloneqq \; \sqrt{ \int_\Omega \|\hatu (t, \mathbf{x})- \uu (t, \mathbf{x})\|^2\, \mathrm{d}\mathbf{x}},
\end{equation*}

which is computed numerically with composite trapezoidal rule. To apply <Ref>, we estimate the number of required subintervals by requiring

\begin{equation*}
\Eint \le \alpha \Eeqexp   ,
\end{equation*}

wherein $\Eeqexp$ is an initial estimation of the error $\Eeq$. We compute this estimation as

\begin{equation*}
\Eeqexp(t) = M \exp ^ {\omega t} \newdelta_0 + M (\mexp^{\omega t} - 1)\frac{\overline{\newdelta} \cdot \| \Omega \| }{\omega}
\end{equation*}

using the average residual during training

\begin{equation*}
\overline{\newdelta} = \frac{1}{N_\mathrm{eq}} \sum_{i=1}^{N_\mathrm{eq}} \| \partial_t \hatu (t_{\mathrm{eq},i}, \mathbf{x}_{\mathrm{eq},i}) - \mathcal{A}\hatu (t_{\mathrm{eq},i}, \mathbf{x}_{\mathrm{eq},i}) \|_{2, \R^n} .
\end{equation*}

The parameter $\alpha$ scales the numerical integration error relative to the expected error $\Eeqexp$ and is chosen constantly as $\alpha = 0.33$ in the following experiments.

Hence, we can compute the number of required subintervals

\begin{equation*}
N_\mathrm{SI} (t) = \left\lceil \sqrt{\frac{\mexp^{\omega t} K \, t^3}{12 \alpha \Eeqexp(t)}} \,\right\rceil .
\end{equation*}

We setup NNs based on TensorFlow [1] using the hyperbolic tangent as the activation function. Problem specific parameters such as network size are mentioned in the relevant subsections. A set of common parameters and settings is listed here: We train using the L-BFGS optimizer [27] with learning rate $0.1$ unless stated otherwise. The NNs rely on the hyperbolic tangent as activation function and the data-driven and physics-informed contributions are equally weighted with $\kappa = 1$ unless explicitly mentioned.

For testing the NN and our error estimator, we generate a grid of equally spaced data points in the complete spatial domain and evaluate the prediction of the NN and the error estimator at various points in time.

§.§ Heat equation

We consider the scalar BIVP

\begin{equation}\label{eq::heat}
\begin{aligned}
\partial_t u(t,x) &= \tfrac{1}{5} \partial_{xx} u(t,x),& \\
u(0, x) &= \sin(2\pi x ) \qquad &\mathrm{for }&\; x\in \Omega, \\
u(t, 0) &= 0 = u(t, 1) \qquad &\mathrm{for }&\; t \in \mathbb{T} .
\end{aligned}
\end{equation}

for the spatio-temporal domain $\timeInt \times \Omega = [0,0.5] \times (0,1)$ with analytical solution $u(t, x) = \sin(2 \pi x) \mexp^{-\frac{1}{5} (2\pi)^2 t}$; see [14] .

Semigroup properties
The linear operator under consideration is defined by

\begin{equation*}
\begin{aligned}
\mathcal{A}\colon D(\mathcal{A}) = H_0^2(\Omega) \to L^2(\Omega),\qquad
u \mapsto \tfrac{1}{5} \partial_{xx} u.
\end{aligned}
\end{equation*}

The operator $\mathcal{A}$ generates a $0 $-contraction semigroup on $L^2(\Omega)$ <cit.>, which is also exponentially stable on $L^2(\Omega)$ using results from [38].

The operator $\mathcal{A} = \tfrac{1}{5}\partial_{xx}$ generates an exponentially stable semigroup on $L^2(\Omega)$ for $\Omega = (0,1)$ with decay parameter $\mu = \frac{1}{5}\pi^2$.

This is an application of <cit.> to the present problem.

For this problem, we use a NN with 4 hidden layers of 10 neurons each and train it for 3000 epochs. The PDE is enforced by 1000 collocation points in the domain $\timeInt \times \Omega$ and the initial condition is trained with 200 data points.

We perform two experiments, one with hard boundary constraints employing <Ref> and one with soft boundary constraints employing <Ref>. In the former case, we design the network with a one-dimensional mask to enforce Dirichlet boundary condition parameterized by $w_1=1$ and $x_1=\tfrac{1}{2}$ (cf. <Ref>). In the latter case, we train the network to match the boundary condition by using the loss contribution

\begin{equation*}
L_\mathrm{space} = \frac{1}{N_\mathrm{space}} \sum_{i=1}^{N_\mathrm{space}} \left( | u(t_i, x=0) |^2 + | u(t_i, x=1)| ^2 \right)
\end{equation*}

with fixed weighting $\rho = 10$.

Numerical results for hard boundary constraints

The trained network agrees well with the known reference solution, i.e., the deviation of the initial values can be limited from above by $\newdelta_0 = 2.9 \cdot 10^{-3}$ and the average residual over all collocation points is given by $\overline{\newdelta} = 2.3 \cdot 10^{-2}$. As shown in <Ref>, knowing that the semigroup is exponentially decaying improves the error significantly. It outperforms the estimator based on the contraction growth bound by a factor of $2.25$ at the end of the considered time interval.

every tick label/.append style=font=
every axis/.append style=font=
[ymode=log, xmin=0, xmax=0.5, ymin=0.00001, xlabel=Time $t$, ylabel=Error, height=5cm, width=9.2cm,grid=both,
legend cell align=left,
legend pos= south east,
legend columns=4,
legend style=draw=none, nodes=font=]
[no marks, ErrorPred, dashed, line width=] table [x=t, y=E_init, col sep=comma] Figures/heat_equation/exp_run_tanh_test_data_t_error_over_time.csv; $\Einit^\mathrm{exp}$;
[no marks, ErrorPred, dotted, line width=] table [x=t, y=E_PI, col sep=comma] Figures/heat_equation/exp_run_tanh_test_data_t_error_over_time.csv; $\Eeq^\mathrm{exp}$;
[no marks, ErrorPred, line width=] table [x=t, y=E_tot, col sep=comma] Figures/heat_equation/exp_run_tanh_test_data_t_error_over_time.csv; $\Einit^\mathrm{exp} + \Eeq^\mathrm{exp}$;
empty legend

[no marks, ErrorRef2, dashed, line width=] table [x=t, y=E_init, col sep=comma] Figures/heat_equation/con_run_tanh_test_data_t_error_over_time.csv;$\Einit^\mathrm{con}$;
[no marks, ErrorRef2, dotted, line width=] table [x=t, y=E_PI, col sep=comma] Figures/heat_equation/con_run_tanh_test_data_t_error_over_time.csv;$\Eeq^\mathrm{con}$;
[no marks, ErrorRef2, line width=] table [x=t, y=E_tot, col sep=comma] Figures/heat_equation/con_run_tanh_test_data_t_error_over_time.csv;$\Einit^\mathrm{con} + \Eeq^\mathrm{con}$;
[no marks, ErrorRef, line width=] table [x=t, y=E_ref, col sep=comma] Figures/heat_equation/con_run_tanh_test_data_t_error_over_time.csv;$\Eref$;

A posteriori error estimator for the heat equation (<ref>) on $\timeInt=[0,0.5]$ with hard boundary constraints. The a posteriori error estimator $E_\mathrm{PI} + \Einit$ is depicted with the reference error $\Eref$. The error estimator for the semigroup of the heat equation as exponentially decaying semigroup is shown in dark blue denoted by superscript $\mathrm{exp}$, while the results for the contractive semigroup are shown in red marked with the superscript $\mathrm{con}$.

Numerical results for soft boundary constraints
Instead of training the heat equation with hard boundary constraints, we now use the extended theory with ISS as formulated in <Ref> to demonstrate that the error bound including ISS statements is valid for the heat equation.

The error prediction is shown in <Ref>. It is visible that the rise in the true error $\Eref$ towards the end of the time interval $[0,0.5]$ is related to an increased error in the boundary condition. This is properly reflected by the contribution $\Ebc$, which is computed using ISS parameters as analyzed in [21].

every tick label/.append style=font=
every axis/.append style=font=
[ymode=log, xmin=0, xmax=0.5, ymin=0.000008, xlabel=Time $t$, ylabel=Error, height=5cm, width=9.2cm,grid=both,
legend cell align=left,
legend pos= south east,
legend columns=2, transpose legend,
legend style=draw=none, nodes=font=]
[no marks, ErrorPred, dashed, line width=] table [x=t, y=E_init, col sep=comma] Figures/heat_equation_soft_bcs/run_tanh_test_data_t_error_over_time.csv; $\Einit$;
[no marks, ErrorPred, dotted, line width=] table [x=t, y=E_PI, col sep=comma] Figures/heat_equation_soft_bcs/run_tanh_test_data_t_error_over_time.csv; $\Eeq$;
[no marks, ErrorPred, dash dot, line width=] table [x=t, y=E_bc, col sep=comma] Figures/heat_equation_soft_bcs/run_tanh_test_data_t_error_over_time.csv; $\Ebc$;
[no marks, ErrorPred,  line width=] table [x=t, y=E_tot, col sep=comma] Figures/heat_equation_soft_bcs/run_tanh_test_data_t_error_over_time.csv; $\Einit + \Eeq + \Ebc$;
[no marks, ErrorRef, line width=] table [x=t, y=E_ref, col sep=comma] Figures/heat_equation_soft_bcs/run_tanh_test_data_t_error_over_time.csv;$\Eref$;

Heat equation with soft boundary constraints: Prediction error estimation including additional error introduced by the boundary error via ISS properties.

§.§ Homogeneous advection equation

We consider the homogeneous advection equation in two space dimensions with one output value

\begin{equation*}
\begin{aligned}
\partial_t u(t,\mathbf{x}) =\;& (\tfrac{1}{5}, \tfrac{1}{2}) \cdot \nabla u(t,\mathbf{x})\quad (t,\mathbf{x}) \in \timeInt \times \Omega, \\
u(0,\mathbf{x}) =\; & \begin{cases}
\tfrac{1}{2} - \| \mathbf{x} \|_1, \quad\quad\;\; &\mathrm{for }\; \| \mathbf{x} \|_1 \le \tfrac{1}{2}, \\
0, \quad\quad &\mathrm{else},
\end{cases}  \\
u(t, -2, \xtwo ) = \;& u(t, 2, \xtwo) \qquad\quad\;\, t\in\timeInt,\; \xtwo \in [-2,2] ,\\
u(t, \xone, -2) = \;& u(t, \xone, 2) \qquad\quad\;\,  t\in\timeInt,\; \xone \in [-2,2] ,\\
\end{aligned}
\end{equation*}

for $\Omega \coloneqq (-2,2)^2$ and ${\timeInt} \coloneqq [0,8]$.

Semigroup properties
To account for the periodic boundary constraints, we introduce the space

\begin{equation*}
X \coloneqq \left\lbrace f \in L^2_\mathrm{loc}(\R^2)\; \left|\,
\begin{aligned}
&f(\mathbf{x}) = f(\mathbf{x} +\textstyle  \sum_{i=1,2} 4 k_i\mathbf{e_i}) \\ &\mathrm{for }\; k_1, k_2 \in \Z , \text{a.e.}\ \mathbf{x} \in \R^2
\end{aligned} \right.\right\rbrace
\end{equation*}

with norm $\|\cdot \|_{L^2{[-2,2]^2}}$, where $\mathbf{e_i}$ denotes the $i$th cartesian unit vectors in $\R^2$.
The operator

\begin{equation*}
\mathcal{A}\colon D(\mathcal{A}) \to L^2(\Omega), \qquad
u \mapsto (\tfrac{1}{5}, \tfrac{1}{2}) \cdot \nabla u.
\end{equation*}

is bounded on $D(\mathcal{A}) =H^1(\Omega)\cap X $ and generates a strongly continuous semigroup on X. The solution to the advection equation is defined by the shift-semigroup

\begin{equation*}
\mathcal{S}_\mathrm{shift, \mathbf{C}}(s) f(\mathbf{x}) = f(\mathbf{x} + s\mathbf{C})
\end{equation*}

with $\mathbf{C}=(\tfrac{1}{5}, \tfrac{1}{2})^T$, which generates an isometry on $D(\mathcal{A})\subseteq L^2_\mathrm{loc}(\R^2)$ <cit.>. In particular, $\mathcal{S}_\mathrm{shift, \mathbf{C}}$ is a $0$-contraction semigroup.

We set up a NN with 8 hidden layers of 40 neurons each and train it for 10000 epochs with $\kappa=0.3$. For the physics-informed contribution 10000 collocation points in the domain $\tilde{\timeInt} \times \Omega$, with $\tilde{\timeInt} \coloneqq [0,4]$, are used. The initial condition is learned on training points located on a grid of $201 \times 201$ for $t=0$. The periodic boundary condition are enforced according to <Ref>sec::methodology.

Numerical results

In <Ref>, it is clearly visible that the prediction of the error lies in the same magnitude as the actual error and is mostly dominated by the error on the initial condition for the trained time range $\tilde{\timeInt} =[0,4]$. Even for an extended time domain ${\timeInt} = [0,8]$, for which the network has not been trained, the error estimator is a true upper bound on the prediction error. Since the estimator does not depend on NNs as a method or the quality of the approximation, this is an expected result.

The true error is overestimated increasingly over time but remains at the same magnitude of $10^0$. The lower part of <Ref> illustrates that leaving the training domain increases both the error and the overestimation of the error.

every tick label/.append style=font=
every axis/.append style=font=
[ymode=log, xmin=0, xmax=8, ylabel=Error, height=5cm, width=9.2cm,grid=both,
legend cell align=left,
legend pos= south east,
legend style=draw=none,
legend style=draw=none, nodes=font=,
[no marks, ErrorPred, dashed, line width=] table [x=t, y=E_init, col sep=comma] Figures/transport_equation_2dim/run_tanh_input_eval__error_over_time.csv;
[no marks, ErrorPred, dotted, line width=] table [x=t, y=E_PI, col sep=comma] Figures/transport_equation_2dim/run_tanh_input_eval__error_over_time.csv;
[no marks, ErrorPred,  line width=] table [x=t, y=E_tot, col sep=comma] Figures/transport_equation_2dim/run_tanh_input_eval__error_over_time.csv;
[no marks, ErrorRef, line width=] table [x=t, y=E_ref, col sep=comma] Figures/transport_equation_2dim/run_tanh_input_eval__error_over_time.csv;
$\Einit$,$\Eeq$ ,$\Einit+ \Eeq$,$\Eref$

[xmin=0, xmax=8,xlabel=Time $t$, ylabel=Overestimation, height=3.5cm, width=9.2cm, grid=both, yshift=-2.5cm,
legend cell align=left,
legend pos= north west,
legend style=draw=none,
legend style=draw=none, nodes=font=, name=axis2]
[no marks, ErrorPred, line width=] table [x=t, y=E_rel, col sep=comma] Figures/transport_equation_2dim/run_tanh_input_eval__error_over_time.csv;
$\frac{\Einit + \Eeq}{\Eref}$

Top: Error estimation for the two dimensional transport equation for ${\timeInt}=[0,8]$. Bottom: The error estimator overestimates the true error only by an acceptable factor even outside the training range $\tilde{\timeInt} =[0,4]$.

§.§ Navier-Stokes equations

The Navier-Stokes equations (NSEs) are one prominent example of numerical challenges and their solution via PINNs has been investigated previously in [23, 2, 43]. We consider velocity-pressure formulation of the NSEs

\begin{equation*}
\begin{aligned}
\dt \mathbf{u} + \varrho (\mathbf{u} \cdot \nabla) \mathbf{u} &= - \nabla p + \frac{1}{Re} \nabla^2 \mathbf{u} & \mathrm{in }\; \timeInt &\times \Omega,\\
\nabla \cdot \mathbf{u} &= 0  &\mathrm{in }\; \timeInt &\times \Omega,\\
%\mathbf{u} &= \mathbf{u_b} &\mathrm{on }\; \timeInt &\times \partial\Omega,\\
%\mathbf{u} &= \mathbf{u_0} &\mathrm{in }\; \{t=0\}& \times \Omega,
\end{aligned} \label{eq::nses}
\end{equation*}

to address a problem for which the analytical solution is known, the Taylor flow [26]. For the Taylor flow, we consider a two dimensional domain $\Omega \coloneqq (0,\pi)^2 $ with associated velocities $\mathbf{u} = (u,v)^T$, pressure $p$,  $Re = 1$ and $\varrho=1$ on the time domain $\timeInt \coloneqq [0,1]$. The analytical solution is then given by [26]

\begin{equation*}
\begin{aligned}
u(t, \xone, \xtwo) =& - \mathrm{cos}(\xone) \sin (\xtwo) \mexp^{-2t}, \\
v(t, \xone, \xtwo) =& \sin(\xone) \cos (\xtwo) \mexp^{-2t}, \\
p(t, \xone, \xtwo) =& - \frac{1}{4}( \mathrm{cos}(2 \xone) + \cos (2 \xtwo) )\mexp^{-4t} .\\
\end{aligned}
\end{equation*}

Semigroup properties
The semigroup properties for the Navier-Stokes equations can be derived by applying the Helmholtz-Leray projection [6] on (<ref>) and investigating the resulting abstract Cauchy problem. The Helmholtz-Leray projection $\mathcal{P}$ is the projection from $L^2(\Omega) $ into the space of divergence-free vector fields $L_\sigma^2(\Omega ) \coloneqq \mathcal{P} L^2(\Omega)$. Using the Stokes operator

\begin{equation*}
\mathcal{A} = \mathcal{P} \Delta \quad, \; D(\mathcal{A}) = H^2(\Omega )^n \cap H^{1}_0(\Omega )^n \cap L_\sigma ^2(\Omega ) ,
\end{equation*}

the projected NSEs are

\begin{equation*}
\begin{aligned}
\dot{\uu}(t)  = \mathcal{A} \uu (t) - \mathcal{P} (\uu\cdot \nabla)\uu ,
\end{aligned}
\end{equation*}

wherein the rightmost term is considered an inhomogeneity. The Stokes operator generates a bounded analytic semigroup and hence a strongly continuous semigroup in $L^2(\Omega)$ with $\Omega \subset \mathbb{R}^d$ (here, $d=2$) bounded with a smooth boundary and Dirichlet boundary condition [15]. In consequence of <cit.> and [15], the operator $\mathcal{P}\Delta$ generates a 0-contraction semigroup on $L^2_\sigma(\Omega)$.

We train a PINN with the loss function (<ref>), $\kappa = 0.7$, $\rho = 100$ and

\begin{equation*}
L_\mathrm{space} = \sum_{i=1}^{N_\mathrm{eq}} \| \nabla \cdot \mathbf{u} \|^2 .
\end{equation*}

Since the PINN does not automatically satisfy the restriction that the target function must lie in the space of divergence-free functions, we hereby achieve that the violations made by our NN approximation are negligible in comparison to the other approximation errors. Due to this modification in the strategy, our estimator is no longer guaranteed to be a rigorous upper bound.
We use hard boundary conditions as discussed in <ref> with centers $\xonezero = \xtwozero = \pi/2$ and widths $\wone = \wtwo = \pi$. The contribution $\boldsymbol{\chi}(t, \mathbf{x})$ contains the target boundary condition. Its components are given by

\begin{equation*}
\begin{aligned}
\chi_{1,2}(t, \mathbf{x}) = & \left(-\tfrac{x_{1,2}^2}{\pi^2}+1 \right)\left(\mp \sin(x_{2,1})\mexp^{-2t} \right) \\
& + \left(-\tfrac{(x_{1,2}-\pi)^2}{\pi^2}+1 \right)\left(\pm \sin(x_{2,1}) \mexp^{-2t} \right).
\end{aligned}
\end{equation*}

Here, we want to emphasize that the pressure $p$ (or, to be more precise, its spatial derivative) and the non-linear terms are considered only as an additional contribution to the residual (c.f. <Ref>). Especially the pressure contribution to the residual is justified since it is not governed by a distinct evolution equation. For the sake of completeness, it is included in the initial condition, but the pressure is not restricted by any boundary condition.

By using methodology as proposed in [24], one could enforce divergence-free fields in NNs by construction. Since this method relies on Fourier transformation and, therefore, would introduce an extensive bias to a sin/cos-based solution as desired, we decide to approximate the divergence-freeness as a soft constraint.

We train a PINN with 10 hidden layers with 80 neurons in each layer over 20000 epochs with 2000 collocation points. The initial condition is trained on a grid of $31 \times 31$ points over the spatial domain.

Numerical results

The trained network agrees well with the initial condition deviating by $\newdelta_0 = 1.2 \cdot 10^{-3}$, the temporal evolution is also reflected well such that $\overline{\newdelta}= 5 \cdot 10^{-3}$. As visible in <Ref>, the actual error decays slightly over time, which is not reflected by the error estimator. Here, the error estimator overestimates the actual error by approximately two magnitudes.

every tick label/.append style=font=
every axis/.append style=font=
[xmin=0, xmax=1,xlabel=Time $t$, ymode=log, ylabel=Error, height=5cm, width=9.2cm, grid=both,
legend cell align=left,
legend columns = 2,
legend style=at=(0.95,0.7), draw=none, nodes=font=,
[no marks, ErrorPred, dashed, line width=] table [x=t, y=E_init, col sep=comma] Figures/nse/run_tanh_taylor_t_error_over_time.csv;
[no marks, ErrorPred, dotted, line width=] table [x=t, y=E_PI, col sep=comma] Figures/nse/run_tanh_taylor_t_error_over_time.csv;
[no marks, ErrorPred, dash dot, line width=] table [x=t, y=E_tot, col sep=comma] Figures/nse/run_tanh_taylor_t_error_over_time.csv;
[no marks, ErrorRef, line width=] table [x=t, y=E_ref, col sep=comma] Figures/nse/run_tanh_taylor_t_error_over_time.csv;

Predicted error for the PINN simulating the Taylor flow goverened by the NSE. For comparison, the reference error is computed by numerical integration of the difference between the analytical solution and the PINN prediction.

§.§ Klein-Gordon equation with homogeneous Neumann boundary

As a final example, we consider the Klein-Gordon equation with Neumann boundary condition, which is relevant to problems in quantum field theory [35]. The one-dimensional problem can be formulated as

\begin{equation}
\label{eq::kg}
\partial_{tt} u = \partial_{xx} u - \frac{1}{4} u
\end{equation}

in the domain $\timeInt \times \Omega = [0, 0.2] \times (0,1)$. The initial values are given as

\begin{equation*}
\begin{aligned}
u(0,x) =& \cos(2 \pi x) ,\\
\partial_t u(0,x) =& \tfrac{1}{2} \cos( 4 \pi x ) ,
\end{aligned}
\end{equation*}

and we equip (<ref>) with homogeneous Neumann boundary conditions, i.e.

\begin{equation*}
\partial_x u (t, x\in\{0,1\}) =\;0 .
\end{equation*}

The analytical solution is

\begin{equation*}
u(t, x) = \cos(2 \pi x) \cos \left(a_2 t\right) +
\tfrac{1}{2 \cdot a_4} \, \cos(4 \pi x) \sin\left( a_4 t\right)
\end{equation*}

with $a_n = \sqrt{(n \pi)^2 +\tfrac{1}{4}}$; cf. [35].

Semigroup properties
The BIVP for the Klein-Gordon equation can be written as a first order system of the form

\begin{equation*}
\partial_t \begin{pmatrix}
u \\ u_t
\end{pmatrix}
= \begin{pmatrix}
0 & \mathcal{I} \\ \Delta_{xx} - \tfrac{1}{4} & 0
\end{pmatrix}
\begin{pmatrix}
u \\ u_t
\end{pmatrix} .
\end{equation*}

According to [25] the Neumann Laplacian generates a C0-cosine operator function on $L^2(\Omega)$ and applying results <cit.>, the operator

\begin{equation*}
\tilde{\mathcal{A}} =  \begin{pmatrix}
0 & \mathcal{I} \\
\mathcal{A} + \mathcal{B} & 0
\end{pmatrix}
\end{equation*}

for a bounded linear operator $\mathcal{B} \colon D(\mathcal{A}) \rightarrow L^2(\Omega)$ generates a strongly continuous semigroup.

Since $\mathcal{B} = -\tfrac{1}{4} \mathcal{I}$ with $\mathcal{I}$ denoting the identity operator, $\mathcal{B}$ suffices the conditions for <cit.>. Hence, the operator $\tilde{\mathcal{A}}$ generates a strongly continuous semigroup on $X = L^2(\Omega) \times L^2(\Omega)$ for the domain $D(\tilde{\mathcal{A}}) = H^2(\Omega) \times L^2(\Omega)$.

We use the results <cit.> to derive the necessary parameters for error estimation. According to these results, the resolvent set of the new operator $\rho(\tilde{\mathcal{A}})$ is given by all $\lambda \in \C$ whenever $\lambda^2 \in \rho(\mathcal{A}+ \mathcal{B})$. This then transfers to the spectrum, such that $\lambda \in \sigma(\tilde{\mathcal{A}})$ whenever $\lambda ^2 \in \sigma(\mathcal{A}+ \mathcal{B})$.

The spectrum of the Neumann Laplacian consists of the point spectrum only and is given by

\begin{equation*}
\sigma(\Delta_{[0,1]}^N ) = \left\lbrace -k \pi \mid k\in \N_0 \right\rbrace
\end{equation*}

so that the spectrum of $\tilde{\mathcal{A}}$ can be derived to be

\begin{equation*}
\sigma(\tilde{\mathcal{A}}) = \left\lbrace \sqrt{-k \pi - 0.25} \mid k\in \N_0 \right\rbrace.
\end{equation*}

Hence, the largest real part of the eigenvalues is given by $\lambda_\mathrm{max} = 0$ for $ k=0$, which we use to apply <Ref> utilizing the spectral mapping theorem <cit.>.

If we consider the Klein-Gordon-equation on the more restrictive space $X = H^2(\Omega) \times L^2(\Omega)$ with norm

\begin{equation}\label{eq::newnorm_kg}
\| (u, v) \| = \left( \| \nabla u\|^2_{L^2(\Omega)} + 0.25 \|u\|^2_{L^2(\Omega)} + \| v\|^2_{L^2(\Omega)}\right)^{1/2},
\end{equation}

then one can show that the generated semigroup is a 0-contraction semigroup analogous to the proof of <cit.> with minor modifications.

Numerical determination of M
This Klein-Gordon example on $X = L^2(\Omega) \times L^2(\Omega)$ shows nicely the technical difficulties arising when the parameter $M$ is required for a complete understanding of the growth behaviour of the weak solution.
To retrieve a reasonable estimate for $M$, we use the inequality

\begin{equation*}
M \ge \frac{\|\mathcal{S}(t) \mathbf{v}\|}{\mexp^{\omega t} \|\mathbf{v}\| }.
\end{equation*}

Here, we use input data $\mathbf{v}$, which is Gaussian noise with means in the range $[0,20] \times [0,20]$ and variances in the range $[0.05, 0.5] \times [0.05, 0.5]$, and compute $\|\mathcal{S}(t) \mathbf{v}\|$ with the python framework FiPy [17]. Thereby we approximate $M \ge 164.43$.

We use a NN with 6 hidden layers with 10 neurons each. The initial data is given by 201 equally spaced data points in the range $[0,1]$ and 2500 collocation points are selected equally distributed from $\timeInt \times [0,1]$. The Neumann boundary condition is implemented as soft boundary constraint via

\begin{equation*}
L_\mathrm{space} = \frac{1}{N_\mathrm{eq}} \sum_{i=1}^{N_\mathrm{eq}} \left( | u_x(t_{\mathrm{eq},i}, 0) |^2 + | u_x(t_{\mathrm{eq},i}, 1) |^2 \right)
\end{equation*}

with $\rho = 50$. The data-driven loss is more heavily considered by choosing $\kappa=2.5$.
We use multiple optimization algorithms to train this NN properly. Firstly, we run 30000 epochs of the adam optimizer with learning rate $10^{-3}$ to find a suitable starting point for further 50000 epochs of the second order L-BFGS algorithm with learning rate $10^{-2}$.

Numerical results
The error estimator for the Klein-Gordon equation (<ref>) on $X = L^2(\Omega)\times L^2(\Omega)$ overestimates the true error from the beginning significantly by two orders of magnitude (c.f. <Ref>, top). This improves slightly over time to only one order of magnitude, but is most clearly due to the high value of $M$. The fact that $M \neq 1$ is reflected nicely in <Ref> (top), in which the error contribution $\Einit$ is significantly larger than the true error $\Eref$  even for $t=0$.

When evaluating the NN on $X= H^2(\Omega)\times L^2(\Omega)$ with the norm (<ref>) we are considering a 0-contraction semigroup, which is properly reflected in <Ref> (bottom). In consequence, the error estimator lies in the same order of magnitude as the true error and overestimates the error by less than $10\%$.

every tick label/.append style=font=
every axis/.append style=font=
[ymode=log, xmin=0, xmax=0.2, ymin=0.006, ylabel=Error, height=5cm, width=9.2cm,grid=both,
legend cell align=left,
xtick = 0, 0.05, 0.1, 0.15, 0.2,
legend pos= south east,
legend columns=2, transpose legend,
legend style=draw=none, nodes=font=]
[no marks, ErrorPred, dashed, line width=] table [x=t, y=E_init, col sep=comma] Figures/klein_gordon/run_tanh_kleingordon_ref_t_error_over_time.csv; $\Einit$;
[no marks, ErrorPred, dotted, line width=] table [x=t, y=E_PI, col sep=comma] Figures/klein_gordon/run_tanh_kleingordon_ref_t_error_over_time.csv; $\Eeq$;
[no marks, ErrorPred,  line width=] table [x=t, y=E_tot, col sep=comma] Figures/klein_gordon/run_tanh_kleingordon_ref_t_error_over_time.csv; $\Einit + \Eeq$;
[no marks, ErrorRef, line width=] table [x=t, y=E_ref, col sep=comma] Figures/klein_gordon/run_tanh_kleingordon_ref_t_error_over_time.csv;$\Eref$;

every tick label/.append style=font=
every axis/.append style=font=
[ymode=log, xmin=0, xmax=0.2, ymin=0.006, xlabel=Time $t$, ylabel=Error, height=5cm, width=9.2cm,grid=both,
legend cell align=left,
xtick = 0, 0.05, 0.1, 0.15, 0.2,
legend pos= south east,
legend columns=2, transpose legend,
legend style=draw=none, nodes=font=]
[no marks, ErrorPred, dashed, line width=] table [x=t, y=E_init, col sep=comma] Figures/klein_gordon_in_H2L2/run_tanh_kleingordon_ref_t_error_over_time.csv; $\Einit$;
[no marks, ErrorPred, dotted, line width=] table [x=t, y=E_PI, col sep=comma] Figures/klein_gordon_in_H2L2/run_tanh_kleingordon_ref_t_error_over_time.csv; $\Eeq$;
[no marks, ErrorPred,  line width=] table [x=t, y=E_tot, col sep=comma] Figures/klein_gordon_in_H2L2/run_tanh_kleingordon_ref_t_error_over_time.csv; $\Einit + \Eeq$;
[no marks, ErrorRef, line width=] table [x=t, y=E_ref, col sep=comma] Figures/klein_gordon_in_H2L2/run_tanh_kleingordon_ref_t_error_over_time.csv;$\Eref$;

Error predicted via the presented methods for the Klein-Gordon equation with homogeneous Neumann boundary condition on $X= L^2(\Omega)\times L^2(\Omega)$ (top) or $X= H^2(\Omega)\times L^2(\Omega)$ (bottom).

§ DISCUSSION

We derived computable rigorous upper bounds for the prediction error of PINNs and other surrogate modeling techniques that approximate solutions to linear partial differential equations. We illustrated the applicability of the error estimator in several examples. We discuss in the following the most important of the thereby illustrated properties, possibilities, and shortcomings.

Firstly, we demonstrated with the NSEs that the error estimator can also be applied to nonlinear partial differential equations. However, this presumes that the dominant contribution to the time evolution is still linear and leads to a weakening of the guaranteed upper error bound to an approximate error estimation. In the most general case, the previously presented methodology is not directly applicable to nonlinear PDEs, which is subject to further investigation.

Apart from the question of how to thoroughly extend the presented theory to nonlinear PDEs, finding the correct governing semigroup parameters remains a challenge. Even though most real-world problems (such as the Navier-Stokes equations, Maxwell equations, etc.) probably have been discussed extensively as part of many mathematical research areas, the investigations of the parameters necessary for the a posteriori error estimator have been more academic in nature. More precisely, while the existence of these parameters has been proven in various constellations, their explicit calculation is either not directly included in the publications or hidden in lengthy proofs. We have presented for the Klein-Gordon equation how to determine numerically one of the key parameters. This approach could be extended to compute both the growth bound $\omega$ and the scaling factor $M$ by discretizing the problem using finite-element methods and finding the eigenvalues. This problem in general and the latter suggestion need further study.

We have shown for the heat equation that the two different estimators resulting once from the analysis that it is a 0-contractive semigroup and once from the exponential stability analysis, overestimate the true error in different magnitudes. This improvement of the estimator through a more thorough analysis of the semigroup is consistent with the fundamental reason why PINNs are more powerful than purely data-driven NNs: the introduction of a priori knowledge.

Mostly, except for the heat equation and the Klein-Gordon equation, we used hard boundary constraints. This is not always possible, and when it is, it can be laborious to find a suitable representation. Hence, the extension of the theory using input-to-state stability (ISS) is highly relevant for the practical use of this estimator. Since ISS is restricted to exponentially stable semigroups, a similar analysis for, e.g., the Klein-Gordon equation remains open.

Beyond the previously mentioned technical limitations and challenges, the methodology does not provide information about a lower limit to the true error, as it would be required for judging the sharpness of the error estimator [42].

[1]
M. Abadi, A. Agarwal, P. Barham, and E. B. et al.
TensorFlow: Large-scale machine learning on heterogeneous systems,
Software available from tensorflow.org.

[2]
C. Amalinadhi, P. S. Palar, R. Stevenson, and L. Zuhal.
On physics-informed deep learning for solving navier-stokes
In Proc. AIAA SciTech, page 1436, 2022.

[3]
F. Arnold and R. King.
State–space modeling for control based on physics-informed neural
Eng. Appl. Artif. Intell., 101:104195, 2021.

[4]
A. R. Barron.
Approximation and estimation bounds for artificial neural networks.
Mach Learn, 14(1):115–133, 1994.

[5]
S. Berrone, C. Canuto, and M. Pintore.
Solving pdes by variational physics-informed neural networks: an a
posteriori error analysis.
arXiv preprint, arXiv:2205.00786, 2022.

[6]
A. J. Chorin and J. E. Marsden.
A mathematical introduction to fluid mechanics, volume 3 of
Texts Appl. Math.
Springer-Verlag, 1993.

[7]
M. A. Cruz, R. Thompson, L. Sampaio, and R. Bacchi.
The use of the reynolds force vector in a physics informed machine
learning approach for predictive turbulence modeling.
Comput. Fluids, 192:104258, 2019.

[8]
T. De Ryck and S. Mishra.
Generic bounds on the approximation error for physics-informed (and)
operator learning.
arXiv preprint 2205.11393, 2022.

[9]
Z. Ding.
A proof of the trace theorem of sobolev spaces on lipschitz domains.
Proc. Amer. Math. Soc., 124(2):591–600, 1996.

[10]
S. Dong and N. Ni.
A method for representing periodic functions and enforcing exactly
periodic boundary conditions with deep neural networks.
J. Comput. Phys., 435:110242, jun 2021.

[11]
C. Duan, Y. Jiao, Y. Lai, X. Lu, Q. Quan, and J. Z. Yang.
Analysis of deep ritz methods for laplace equations with dirichlet
boundary conditions.
arXiv preprint 2111.02009, 2021.

[12]
K.-J. Engel, R. Nagel, and S. Brendle.
One-parameter semigroups for linear evolution equations, volume
Springer-Verlag, 2000.

[13]
L. C. Evans.
Partial differential equations, volume 19.
Amer. Math. Soc., 2010.

[14]
J. Fourier.
Analytical Theory of Heat.
Cambridge University Press, 1878.

[15]
Y. Giga.
The stokes operator in $ l\_r $ spaces.
Proc. Jpn. Acad. A: Math. Sci., 57(2):85–89, 1981.

[16]
S. Greydanus, M. Dzamba, and J. Yosinski.
Hamiltonian neural networks.
Neural. Inf. Process. Syst., 32, 2019.

[17]
J. E. Guyer, D. Wheeler, and J. A. Warren.
Fipy: Partial differential equations with python.
Comput. Sci. Eng., 11(3):6–15, 2009.

[18]
B. Hillebrecht and B. Unger.
Certified machine learning: A posteriori error estimation for
physics-informed neural networks.
In 2022 Int. Joint Conf. on Neural Networks (IJCNN), pages
1–8, 2022.

[19]
Q. Hong, J. W. Siegel, and J. Xu.
A priori analysis of stable neural network solutions to numerical
arXiv preprint 2104.02903, 2021.

[20]
K. Hornik, M. Stinchcombe, and H. White.
Universal approximation of an unknown mapping and its derivatives
using multilayer feedforward networks.
Neural Netw., 3(5):551–560, 1990.

[21]
B. Jacob, R. Nabiullin, J. R. Partington, and F. L. Schwenninger.
Infinite-dimensional input-to-state stability and orlicz spaces.
SIAM J. Cont. Optim., 56(2):868–889, 2018.

[22]
P. Jin, L. Lu, Y. Tang, and G. E. Karniadakis.
Quantifying the generalization error in deep learning in terms of
data distribution and neural network smoothness.
Neural Netw., 130:85–99, 2020.

[23]
X. Jin, S. Cai, H. Li, and G. Karniadakis.
Nsfnets (navier-stokes flow nets): Physics-informed neural networks
for the incompressible navier-stokes equations.
J. Comput. Phys., 426:109951, 2021.

[24]
K. Kashinath, P. Marcus, et al.
Enforcing physical constraints in cnns through differentiable pde
In ICLR Workshop DeepDiffEq, 2020.

[25]
V. Keyantuo and M. Warma.
The wave equation in lp-spaces.
In Semigroup Forum, volume 71, pages 73–92. Springer-Verlag,

[26]
J. Kim and P. Moin.
Application of a fractional-step method to incompressible
navier-stokes equations.
J. Comput. Phys., 59(2):308–323, 1985.

[27]
D. Liu and J. Nocedal.
On the limited memory BFGS method for large scale optimization.
Math. Program., 45:503–528, 1989.

[28]
K. S. McFall.
An artificial neural network method for solving boundary value
problems with arbitrary irregular boundaries.
PhD thesis, Georgia Institute of Technology, 2006.

[29]
K. McFall and J. Mahan.
Artificial neural network method for solution of boundary value
problems with exact satisfaction of arbitrary boundary conditions.
IEEE Trans. Neural Netw., 20(8):1221–1233, 2009.

[30]
C. Meng, S. Seo, D. Cao, S. Griesemer, and Y. Liu.
When physics meets machine learning: A survey of physics-informed
machine learning, 2022.

[31]
P. Minakowski and T. Richter.
Error estimates for neural network solutions of partial differential
arXiv preprint 2107.11035, 2021.

[32]
D. Mugnolo.
Semigroup methods for evolution equations on networks,
volume 20.
Springer-Verlag, 2014.

[33]
J. Nicodemus, J. Kneifl, J. Fehr, and B. Unger.
Physics-informed neural networks-based model predictive control for
multi-link manipulators.
IFAC-PapersOnLine, 55(20):331–336, 2022.

[34]
A. Pazy.
Semigroups of Linear Operators and Applications to Partial
Differential Equations, volume 44.
Springer-Verlag, 1989.

[35]
A. D. Polyanin.
Handbook of linear partial differential equations for engineers
and scientists.
Chapman and Hall, 2001.

[36]
M. Raissi, P. Perdikaris, and G. Karniadakis.
Physics-informed neural networks: A deep learning framework for
solving forward and inverse problems involving nonlinear partial differential
J. Comput. Phys., 378:686–707, 2019.

[37]
F. L. Schwenninger.
Input-to-state stability for parabolic boundary control: linear and
semilinear systems.
In Control theory of infinite-dimensional systems, pages
83–116. Springer-Verlag, 2020.

[38]
C. Seifert, S. Trostorff, and M. Waurick.
Exponential Stability of Evolutionary Equations.
Springer-Verlag, 2022.

[39]
Y. Shin, Z. Zhang, and G. E. Karniadakis.
Error estimates of residual minimization using neural networks for
linear pdes.
arXiv preprint, arXiv: 2010.08019, 2020.

[40]
E. D. Sontag et al.
Smooth stabilization implies coprime factorization.
IEEE Trans. Automat. Control, 34(4):435–443, 1989.

[41]
M. Tucsnak and G. Weiss.
Observation and control for operator semigroups.
Springer-Verlag, 2009.

[42]
R. Verfürth.
A review of a posteriori error estimation and adaptive
mesh-refinement techniques.
Wiley-Teubner Ser. Adv. Numer. Math., 1996.

[43]
S. Wang, Y. Teng, and P. Perdikaris.
Understanding and mitigating gradient flow pathologies in
physics-informed neural networks.
SIAM J. Sci. Comput., 43(5):A3055–A3081, 2021.

[44]
J.-L. Wu, H. Xiao, and E. Paterson.
Physics-informed machine learning approach for augmenting turbulence
models: A comprehensive framework.
Phys. Rev. Fluids, 3:074602, Jul 2018.

[45]
H. Xiao, J.-L. Wu, J.-X. Wang, and E. Paterson.
Physics-informed machine learning for predictive turbulence modeling:
Progress and perspectives.
Proc. AIAA SciTech, 2017.

[46]
J. Yu, L. Lu, X. Meng, and G. E. Karniadakis.
Gradient-enhanced physics-informed neural networks for forward and
inverse pde problems.
Comput. Meth. Appl. Mech. Eng., 393:114823, 2022.

§.§ Boundary homogenization

In case, the boundary condition of <Ref> are inhomogeneous, meaning $\mathbf{u_b} \not\equiv 0 $, homogenization of the boundary facilitates semigroup reasoning. We aim at finding a mild solution to the BIVP (<ref>) in a Banach space X.

Here, we can assume that there exist $\mathbf{w} \in D(\mathcal{A}) \cap D(\mathcal{B}) \subseteq X$ such that $ \mathcal{B}\mathbf{w} = \mathbf{u_b}$
and $\partial_t \mathbf{w} \in X$ exists. In case this does not exist, there would not be a mild solution to the BIVP (<ref>), and hence not of concern.

We define $\mathbf{v} = \mathbf{u} - \mathbf{w} \in X$ which then fulfills the BIVP

\begin{equation} \label{eq::pde_homogeneous_bc}
\begin{aligned}
\partial_t \mathbf{v} =\,& \mathcal{A} \mathbf{v} + \mathbf{f_w} \quad& \mathrm{in}\; \timeInt \times \Omega, \\
\mathbf{v} =\,& \mathbf{u_0} - \mathbf{w_0} \quad& \mathrm{in } \;\{t = 0 \} \times \Omega ,\\
\mathcal{B} \mathbf{v} =\,& 0 \quad& \mathrm{in }\; \timeInt \times \partial \Omega ,
\end{aligned}
\end{equation}

with $\mathbf{w}_0 = \mathbf{w}(0, \cdot)$ being the state at $t=0$ of the boundary condition enforcing function $\mathbf{w}$.
To compensate for contributions to the temporal evolution made by $\mathbf{w}$ we introduce

\begin{equation*}
\mathbf{f_w}= \mathcal{A}\mathbf{w} - \partial_t \mathbf{w} \in X.
\end{equation*}

Under <Ref>, a mild solution to (<ref>) corresponding to the weak solution according to <cit.> is given by

\begin{equation}
\mathbf{v} = \mathcal{S}(t) (\mathbf{u_0} - \mathbf{w_0}) + \int_0^t \mathcal{S}(t-s) \mathbf{f_w}(s) \ds .
\end{equation}

Since $\mathbf{w} \in X$, a weak solution to <Ref> is then given by

\begin{equation}
\mathbf{u} = \mathbf{w} + \mathcal{S}(t) (\mathbf{u_0} - \mathbf{w_0}) + \int_0^t \mathcal{S}(t-s) \mathbf{f_w}(s) \ds .
\end{equation}

We emphasize that we assume that $\mathbf{w}$ fulfills the boundary condition exactly and is sufficiently smooth $\mathbf{w} \in X$.

§.§ Exponential Stability of the Heat Equation

Here, we address the statement given in <Ref>.

According to [38], the general heat equation

\begin{equation*}
\partial_t u = k \partial_{xx} u
\end{equation*}

with diffusion constant $k$ can be brought into the form

\begin{equation}
\begin{split}
\left( \partial_{t, \nu}  \begin{pmatrix}
M_0 & 0 \\ 0& 0
\end{pmatrix} +
\begin{pmatrix}
0 & 0 \\ 0& M_1(\partial_{t, \nu})
\end{pmatrix}  + \right. \\
\left.
\begin{pmatrix}
0 & \mathrm{div} \\ \mathrm{grad}& 0
\end{pmatrix}
\right) \begin{pmatrix}
u \\ \phi
\end{pmatrix} = \begin{pmatrix}
0 \\ 0
\end{pmatrix}
\end{split}
\end{equation}

with $M_0 = 1$, $M_1(\partial_{t,\nu}) = k^{-1}$. Here, [38] uses exponentially weighted $L_2$-spaces with appropriate definitions of partial derivatives $\partial_{t, \nu}$. Since we already know as stated in <Ref>, that the generated semigroup is strongly continuous in t, the weighted partial derivative $\partial_{t, \nu}$ reduces to the conventional derivative.

According to [38], the semigroup of operators is exponentially stable with decayrate

\begin{equation*}
\rho_0 = \min\left\lbrace \rho_1, \frac{c_1}{\|M_1\|^2_{\infty, Re>-\rho_1} \|M_0\| \|C^{-1}\|^2} \right\rbrace .
\end{equation*}

With known $M_0$, $M_1$, it remains to determine $\|C^{-1}\|$, $c_1$ and $\rho_1$. The latter one is defined via the spectral abscissa of $M_1$, which is defined as

\begin{equation}
\begin{split}
\sigma_b(M_1) := \inf \left\lbrace \nu \in \mathbb{R}; \mathbb{C}_{Re>\nu} \subset D(\mathcal{M}_1) \; \right. \\ \left. \mathrm{and} \; \|M_1\|_{\infty, \mathbb{C}_{Re>\nu} }< \infty \right\rbrace
\end{split}
\end{equation}

\begin{equation*}
\sigma_b(M_1) < -\rho_1 \quad \mathrm{for} \;\; \rho_1>0 .
\end{equation*}

Since $\rho_1$ can be chosen arbitrarily large and it is irrelevant for the determination of $\rho_0$.

The parameter $c_1$ can be determined according to <cit.> to fulfill

\begin{equation*}
\mathrm{Re} M_1(z) \ge c_1,
\end{equation*}

for all $z \in \mathbb{C}_{Re>-\rho_1}$ which means that here $c_1 = k^{-1}$.

To determine $\|C^{-1}\|$, we use the basis $\phi_n, \psi_n$ of $L^2([0,1])$ with $\phi_n$ as above and $\psi_n = \sqrt{2} \mathrm{cos}(n\pi x)$.
Thus, we have

\begin{equation*}
\|C^{-1}\| = \max \left\lbrace \max_{n\in \mathbb{N}} (\mathrm{grad}^{-1} \phi_n)  ,\;\max_{n\in \mathbb{N}} (\mathrm{grad}^{-1} \psi_n) \right\rbrace  = \frac{1}{\pi}  .
\end{equation*}

Consequentially, the semigroup of operators is exponentially stable with decay rate

\begin{equation*}
\rho_0 = \frac{k^{-1}}{k^{-2} \cdot 1 \cdot \|C^{-1}\|^2} = k \pi ^2
\end{equation*}

which agrees with the decay analysis yielded by using the principal eigenvalue of the Laplace operator on that range (c.f. <cit.>)