# On the Compatibility between Neural Networks and Partial Differential
Equations for Physics-informed Learning

Kuangdai Leng  Jeyan Thiyagalingam

###### Abstract

We shed light on a pitfall and an opportunity in physics-informed neural
networks (PINNs). We prove that a multilayer perceptron (MLP) only with ReLU
(Rectified Linear Unit) or ReLU-like Lipschitz activation functions will
always lead to a vanished Hessian. Such a network-imposed constraint
contradicts any second- or higher-order partial differential equations (PDEs).
Therefore, a ReLU-based MLP cannot form a permissible function space for the
approximation of their solutions. Inspired by this pitfall, we prove that a
linear PDE up to the $n$-th order can be strictly satisfied by an MLP with
$C^{n}$ activation functions when the weights of its output layer lie on a
certain hyperplane, as called the out-layer-hyperplane. An MLP equipped with
the out-layer-hyperplane becomes “physics-enforced”, no longer requiring a
loss function for the PDE itself (but only those for the initial and boundary
conditions). Such a hyperplane exists not only for MLPs but for any network
architecture tailed by a fully-connected hidden layer. To our knowledge, this
should be the first PINN architecture that enforces point-wise correctness of
PDEs. We show a closed-form expression of the out-layer-hyperplane for second-
order linear PDEs, which can be generalised to higher-order nonlinear PDEs.

###### keywords:

partial differential equation , deep learning , physics-informed neural
network , ReLU activation function

††journal: MASKED

[inst1]organization=Scientific Computing Department, STFC,
addressline=Rutherford Appleton Laboratory, city=Didcot, postcode=OX11 0QX,
country=UK

[height=4.5, nodespacing=14mm, nodesize=25pt, layerspacing=23mm] [count=3,
bias=false, title=Input, text=$x^{]}_{\hiddenlayer}$[count=3, bias=false,
title=First
hidden
layer, text=$x^{]}_{\linklayers}$[count=4, bias=false, title=Last
hidden
layer, text=$x^{]}_{\linklayers}$[count=1, title=Output
layer, text=$\tilde{u}$[thick, ->] (7.2,-2.9) – (8.2,-1.3); [align=center] at
(9.1, -1.0) Data loss ;

[thick, ->] (7.5,-3.3) – (8.1,-3.3); [align=center] at (9.1, -4.9) AutoDiff;
[fill=blue!35, rounded corners=.4cm, align=center, inner sep=.2cm] at (9.1,
-3.3) $\dfrac{\partial\tilde{u}}{\partial\mathbf{x}^{0}}$,

$\dfrac{\partial}{\partial\mathbf{x}^{0}}\dfrac{\partial\tilde{u}}{\partial\mathbf{x}^{0}}$,

$\cdots$ ; [thick, ->] (10.1,-3.3) – (10.7,-3.3);

[from layer=2, from node=1, to layer=3, to node=1, label=$w^{2}_{\link}$rom
layer=2, from node=2, to layer=3, to node=1, label=$w^{2}_{\link}$rom layer=2,
from node=3, to layer=3, to node=1, label=$w^{2}_{\link}$rom layer=2, from
node=4, to layer=3, to node=1, label=$w^{2}_{\node}$lign=center] at (11.5,
-3.3) ~~PDE loss~~
IC loss
BC loss ;

[align=center] at (2.3,-6.2) $\underbrace{\quad\quad\quad}$; [fill=pink,
rounded rectangle, rounded rectangle east arc=none, align=center, anchor=east]
at (2.35,-6.9) $\mathbf{W}^{1}$; [fill=blue!35, rounded rectangle, rounded
rectangle west arc=none, align=center, anchor=west] at (2.3,-6.9)
$\,\mathbf{z}^{1}$; [align=center] at (2.3,-7.5) $\downarrow$;

[align=center] at (4.6,-6.2) $\underbrace{\quad\quad\quad}$; [fill=pink,
rounded rectangle, rounded rectangle east arc=none, align=center, anchor=east]
at (4.65,-6.9) $\mathbf{W}^{2}$; [fill=blue!35, rounded rectangle, rounded
rectangle west arc=none, align=center, anchor=west] at (4.6,-6.9)
$\,\mathbf{z}^{2}$; [align=center] at (4.6,-7.5) $\downarrow$;

[fill=green!40, rounded rectangle, align=center, inner sep=0.2cm] at
(-.1,-8.2) PDE
coefficients; [align=center] at (1.2,-8.2) $\rightarrow$; [draw, align=center,
inner sep=.2cm] at (3.5,-8.2) Steps 1$\sim$5 in Algorithm 1; [align=center] at
(5.8,-8.2) $\rightarrow$;

[fill=blue!35, rounded rectangle, align=center] at (7.6,-8.2)
$\mathbf{g}^{[1]},\mathbf{g}^{[2]},\mathbf{g}^{[3]},\mathbf{g}^{[4]}$;

[fill=pink, rounded rectangle, align=center, inner sep=.2cm] at (7.8,-6.9)
Trainables
$\lambda^{[2]},\lambda^{[3]},\lambda^{[4]}$;

[align=center] at (9.4,-8.1) $\nearrow$;

[align=center] at (9.4,-7.1) $\searrow$;

[fill=orange!50, rounded rectangle,inner sep=.2cm, align=center] at
(10.3,-7.6) $\mathbf{w}^{2}$;

[->, thick, orange] (10.3,-7.2) to [out=70,in=-70] (5.8,-4.8);

We prove that a multilayer perceptron (MLP) with only ReLU-like activation
functions will always lead to a vanished Hessian, and thus cannot provide a
feasible solution space for any second- or higher-order partial differential
equations (PDEs).

We prove that an MLP with $C^{n}$ activation functions can _always_ (i.e.,
regardless of data) satisfy a $n$-th-order linear PDE _exactly_ (i.e., zero
generalisation error) if the weights of its output layer lie on a hyperplane
decided by the given PDE.

We give a closed-form expression of this hyperplane for second-order linear
PDEs along with an implementation. To the best of our knowledge, this should
be the first network architecture that enforces point-wise correctness of PDEs
(instead of their initial or boundary conditions).

##  1 Introduction

Simulation and inversion of partial differential equations (PDEs) play a
crucial role in applied physics and scientific computing. With conventional
numerical methods being computationally expensive, neural networks and deep
learning technologies, in particular the physics-informed neural networks
(PINNs), have made a new and promising route for simulating and inverting
PDEs. Originally formalised by Raissi et al. [1], PINNs have found numerous
applications across various domains [2, 3], featuring remarkable technical
advancements from many aspects, such as the extension to integro-differential
equations [4], graph neural operators [5], identification of nonlinear
operators [6], multi-fidelity models [7], scalable learning in subdomains [8]
and variational inference [9]. Different from data-driven surrogate models
[10] aimed only at minimising misfit to data, PINNs are trained also to
minimise the physics-based loss functions for the governing equation and the
initial and boundary conditions. These PDE-based loss functions are expected
to aid the neural networks in understanding the underlying physics whereby to
achieve a lower generalisation error with a reduced amount of training data. A
vanilla PINN using a multilayer perceptron (MLP) architecture is shown in Fig.
1, where the key is to utilise the utmost strength of deep neural networks in
automatic differentiation (AutoDiff).

[height=4.5, nodespacing=14mm, nodesize=25pt, layerspacing=23mm] [count=3,
bias=false, title=Input, text=$x^{]}_{\hiddenlayer}$[count=3, bias=false,
title=First
hidden
layer, text=$x^{]}_{\linklayers}$[count=4, bias=false, title=Last
hidden
layer, text=$x^{]}_{\linklayers}$[count=1, title=Output
layer, text=$\tilde{u}$[thick, ->] (7.2,-2.9) – (8.2,-1.3); [align=center] at
(9.1, -1.0) Data loss ;

[thick, ->] (7.5,-3.3) – (8.1,-3.3); [align=center] at (9.1, -4.9) AutoDiff;
[fill=blue!35, rounded corners=.4cm, align=center, inner sep=.2cm] at (9.1,
-3.3) $\dfrac{\partial\tilde{u}}{\partial\mathbf{x}^{0}}$,

$\dfrac{\partial}{\partial\mathbf{x}^{0}}\dfrac{\partial\tilde{u}}{\partial\mathbf{x}^{0}}$,

$\cdots$ ; [thick, ->] (10.1,-3.3) – (10.7,-3.3);

[align=center] at (11.5, -3.3) PDE loss
IC loss
BC loss ;

Figure 1: A vanilla PINN with an MLP architecture. The input of the MLP or
$\mathbf{x}^{0}$ contains the independent variables of the PDE, and the
output, $\tilde{u}(\mathbf{x}^{0})$, is the wanted approximate solution to the
PDE. The acronyms “IC” and “BC” are used to denote initial and boundary
conditions, respectively.

Despite their initial success across many domains, PINNs are still far from
being qualified as a routine for solving or inverting PDEs. Many previous
studies have reported or been motivated by the underperformance of a vanilla
PINN (i.e., an end-to-end MLP) from different perspectives. Examples relevant
to our study here include: training in the Fourier space for better
consistency between a PINN and a target PDE in terms of their
parameterisations of solutions [11, 12], using multi-scale or multi-frequency
neural networks [13, 14] or domain decomposition techniques [8, 15] to enhance
the learnability at high frequencies [16], using adaptive activation functions
[17] or adaptive data sampling [18, 19] to address unbalanced gradients of the
loss functions [20, 21], and explicitly embedding relevant laws or conditions
of physics into a network architecture, such as boundary conditions [22, 23,
24], physical symmetry [25] and invariance [26]. These studies more or less
boil down to one fundamental question: _given a data regime, to what extent a
neural network can be compatible with the target PDE system_.

For clarity, let us consider the following function spaces. Let
$\mathcal{F}_{\mathrm{PDE}}$, $\mathcal{F}_{\mathrm{IC}}$,
$\mathcal{F}_{\mathrm{BC}}$ be the function spaces that respectively satisfy
the governing equation, the initial conditions and the boundary conditions.
Given a well-posed PDE system, their intersection
$\mathcal{F}_{\mathrm{PDE}}\cap\mathcal{F}_{\mathrm{IC}}\cap\mathcal{F}_{\mathrm{BC}}$
must be a skeleton set, i.e., it must contain one and only one element, as
denoted by $u$, which is the true solution. A neural network also form a
function space parameterised by its weights, as denoted by
$\mathcal{F}_{\mathrm{NN}}$. Clearly, a neural network can be a good candidate
for physics-informed learning only if $u\in\mathcal{F}_{\mathrm{NN}}$ within
the data regime, of which the necessary conditions are
$\mathcal{F}_{\mathrm{NN}}\cap\mathcal{F}_{\mathrm{PDE}}\neq\emptyset$,
$\mathcal{F}_{\mathrm{NN}}\cap\mathcal{F}_{\mathrm{IC}}\neq\emptyset$ and
$\mathcal{F}_{\mathrm{NN}}\cap\mathcal{F}_{\mathrm{BC}}\neq\emptyset$.
Meanwhile, for a neural network to work efficiently (in terms of data quantity
and convergence rate), the following differences,
$\mathcal{F}_{\mathrm{NN}}-\mathcal{F}_{\mathrm{PDE}}$,
$\mathcal{F}_{\mathrm{NN}}-\mathcal{F}_{\mathrm{IC}}$ and
$\mathcal{F}_{\mathrm{NN}}-\mathcal{F}_{\mathrm{BC}}$, should be as small as
possible. The above mentioned studies can be well associated with these two
purposes; for example, [13, 14, 8, 15, 17] can be understood as increasing the
expressiveness of $\mathcal{F}_{\mathrm{NN}}$ to capture localised features in
$u$ at high frequencies, [11, 12, 25, 26] as decreasing
$\mathcal{F}_{\mathrm{NN}}-\mathcal{F}_{\mathrm{PDE}}$ by changing basis
functions or imposing energy conservation, and [22, 23, 24] as enforcing
$\mathcal{F}_{\mathrm{NN}}-\mathcal{F}_{\mathrm{BC}}=\emptyset$ or
$\mathcal{F}_{\mathrm{NN}}\in\mathcal{F}_{\mathrm{BC}}$.

In this paper, we will show a case of incompatibility between
$\mathcal{F}_{\mathrm{NN}}$ and $\mathcal{F}_{\mathrm{PDE}}$
($\mathcal{F}_{\mathrm{NN}}\cap\mathcal{F}_{\mathrm{PDE}}=\emptyset$) and
provide a simple and novel approach to enforce full compatibility between them
($\mathcal{F}_{\mathrm{NN}}\in\mathcal{F}_{\mathrm{PDE}}$), both regardless of
data. We will prove that an MLP with only Rectified Linear Unit (ReLU) or
ReLU-like Lipschitz activation functions will always cause a vanished Hessian,
which contradicts any second- or higher-order PDEs. In other words, any
approximate solution yielded by a ReLU-based MLP is non-permissible by such
PDEs. While this incompatibility has been explicitly mentioned and avoided in
some applications [27, 28], we do see quite a few applications having
(probably) falling into this trap. Notably, He et al. [29] has shown that
ReLU-based MLPs are formally equivalent to the piecewise linear interpolation
used in the finite element methods [30], implying that a ReLU-based MLP may
still be available for second- or higher-order PDEs. We argue that such
equivalence is only apparent, and the two methods are essentially different in
terms of the locality of function parameterisation; to take advantage of such
formal equivalence, certain measures (e.g., learning in subdomains [8, 15])
must be taken to bridge the gap between the two. We will revisit this issue in
latter part of this paper.

As inspired by the above ReLU-induced incompatibility, we will prove that a
_linear_ PDE up to the $n$-th order can be _strictly_ satisfied by an MLP with
$C^{n}$ activation functions when the weights of its output layer lie on a
hyperplane decided by the PDE, which we refer to as the _out-layer-
hyperplane_. Such a hyperplane exists for any network architecture tailed by a
fully-connected hidden layer. A PINN equipped with the out-layer-hyperplane
becomes “physics-enforced”, or
$\mathcal{F}_{\mathrm{NN}}\in\mathcal{F}_{\mathrm{PDE}}$, no longer requiring
a loss function for the governing equation. To the best of our knowledge, this
should be the first network architecture that enforces point-wise correctness
of the governing equation. We will provide a closed-form expression of the
out-layer-hyperplane for second-order linear PDEs approximated by an MLP and
provide an implementation. We will also discuss its generalisation to other
network architectures and higher-order nonlinear PDEs.

The rest of this paper is organised as follows. In Section 2, we derive the
Hessian of an MLP-approximated solution, followed by Section 3 where we
discuss the incompatibility caused by ReLU. Next, we propose the concept of
out-layer-hyperplane in Section 4, along with its proof and realisation. After
a brief summary of supportive codes in Section 5, we conclude this paper in
Section 6. As for notations adopted in this paper, we use superscripts to
denote layer indices (not powers, except those on $\mathbb{R}$) and subscripts
to denote element indices. Einstein summation notation is assumed for
subscript indices except those put in parentheses. We also omit the dot symbol
for inner product. For example, a matrix multiplication
$\mathbf{A}\cdot\mathbf{B}$ will be written as
$\mathbf{A}\mathbf{B}=A_{ij}B_{jk}=\sum_{j}A_{i(j)}B_{(j)k}$.

##  2 Hessian of an MLP approximation

Let $u(\mathbf{x}^{0})\in\mathbb{R}$ denote the true solution of a PDE, where
$\mathbf{x}^{0}\in\mathbb{R}^{d^{0}}$ are the independent variables such as
the spatial and temporal coordinates. Taking a 2D homogeneous wave equation on
a membrane for example, we have $\mathbf{x}^{0}=\\{x,y,t\\}$ with $d^{0}=3$,
and $u(\mathbf{x}^{0})$ satisfies the governing equation
$u_{xx}+u_{yy}-u_{tt}/c^{2}=0$ (with $c$ being the wave velocity) and certain
initial and boundary conditions. For simplicity, here we assume scalar-valued
PDEs, but our conclusions hold for vector-valued ones (e.g.,
$\mathbf{u}(\mathbf{x}^{0})\in\mathbb{R}^{3}$).

Let $\tilde{u}(\mathbf{x}^{0})$ be a solution approximated by an MLP that has
$L$ layers, the $k$-th layer of which has $d^{k}$ neurons with weights
$\mathbf{W}^{k}\in\mathbb{R}^{d^{k}\times d^{k-1}}$ and biases
$\mathbf{b}^{k}\in\mathbb{R}^{d^{k}}$. This approximate solution can be
formulated as the following composition of functions:

$\tilde{u}(\mathbf{x}^{0})=\mathcal{L}^{L}\circ\sigma^{L-1}\circ\mathcal{L}^{L-1}\circ\cdots\circ\sigma^{2}\circ\mathcal{L}^{2}\circ\sigma^{1}\circ\mathcal{L}^{1}(\mathbf{x}^{0}),$
(1)

where $\mathcal{L}^{k}$ is the affine transformation endowed by the $k$-th
layer:

$\mathbf{z}^{k}:=\mathcal{L}^{k}(\mathbf{x}^{k-1})=\mathbf{W}^{k}\mathbf{x}^{k-1}+\mathbf{b}^{k},$
(2)

for $\mathbf{x}^{k-1}\in\mathbb{R}^{d^{k-1}}$ and
$\mathbf{z}^{k}\in\mathbb{R}^{d^{k}}$, and $\sigma^{k}$ is the activation
function of the $k$-th layer,

$x^{k}_{i}=\sigma^{k}(z^{k}_{i}),\quad i\in\\{1,2,\cdots,d^{k}\\}.$ (3)

Because we assume $d^{L}=1$, we write $\mathbf{W}^{L}$ as $\mathbf{w}^{L}$ for
a consistent notation.

Based on the chain rule, the Jacobian of $\tilde{u}$, i.e.,
$\nabla\tilde{u}\in\mathbb{R}^{d^{0}}$, is given by

$\nabla\tilde{u}=\mathbf{w}^{L}\mathbf{F}^{L-1}\mathbf{W}^{L-1}\cdots\mathbf{F}^{2}\mathbf{W}^{2}\mathbf{F}^{1}\mathbf{W}^{1},$
(4)

where $\mathbf{F}^{k}\in\mathbb{R}^{d^{k}\times d^{k}}$ is a diagonal matrix
(it is diagonal because $\sigma^{k}$ is an element-wise operation) containing
the first derivatives of $\sigma^{k}(z)$:

$\mathbf{F}^{k}:=\frac{\partial\mathbf{x}^{k}}{\partial\mathbf{z}^{k}}=\mathrm{diag}\
\Big{\\{}\left.\frac{d\sigma^{k}}{dz}\right|_{z=z^{k}_{1}},\left.\frac{d\sigma^{k}}{dz}\right|_{z=z^{k}_{2}},\cdots,\left.\frac{d\sigma^{k}}{dz}\right|_{z=z^{k}_{3}}\Big{\\}}.$
(5)

For the purpose of generalising our results to other network architectures and
PDE classes, we rewrite $\nabla\tilde{u}$ in the following top-down manner:

$\nabla\tilde{u}=\mathbf{w}^{L}\mathbf{U}=\mathbf{w}^{L}\mathbf{P}^{k}\mathbf{F}^{k}\mathbf{Q}^{k},\quad\forall
k\in\\{1,2,\cdots,L-1\\},$ (6)

for which we define

$\displaystyle\mathbf{P}^{k}:=$ $\displaystyle\
\frac{\partial\mathbf{x}^{L-1}}{\partial\mathbf{x}^{k}}=\mathbf{F}^{L-1}\mathbf{W}^{L-1}\cdots\mathbf{F}^{k+1}\mathbf{W}^{k+1}\in\mathbb{R}^{d^{L-1}\times
d^{k}},$ (7) $\displaystyle\mathbf{Q}^{k}:=$ $\displaystyle\
\frac{\partial\mathbf{z}^{k}}{\partial\mathbf{x}^{0}}=\mathbf{W}^{k}\mathbf{F}^{k-1}\mathbf{W}^{k-1}\cdots\mathbf{F}^{2}\mathbf{W}^{2}\mathbf{F}^{1}\mathbf{W}^{1}\in\mathbb{R}^{d^{k}\times
d^{0}},$ (8) $\displaystyle\mathbf{U}:=$ $\displaystyle\
\frac{\partial\mathbf{x}^{L-1}}{\partial\mathbf{x}^{0}}=\mathbf{P}^{k}\mathbf{F}^{k}\mathbf{Q}^{k}\in\mathbb{R}^{d^{L-1}\times
d^{0}},\quad\forall k\in\\{1,2,\cdots,L-1\\}.$ (9)

To show the Hessian of $\tilde{u}$, i.e.,
$\nabla\nabla\tilde{u}\in\mathbb{R}^{d^{0}\times d^{0}}$, we first express eq.
(6) in index notation,

$\tilde{u}_{,m}=w_{i}^{L}P_{ij}^{k}F_{jl}^{k}Q_{lm}^{k},\quad\forall
k\in\\{1,2,\cdots,L-1\\},$ (10)

where $()_{,i}:=\frac{\partial}{\partial x_{i}^{0}}()$. The Hessian of
$\tilde{u}$ can be shown by the following steps:

$\displaystyle\tilde{u}_{,mn}=\frac{\partial\tilde{u}_{,m}}{\partial
x_{n}^{0}}=$ $\displaystyle\ w_{i}^{L}\sum_{k=1}^{L-1}P_{ij}^{k}\frac{\partial
F_{jl}^{k}}{\partial x^{0}_{n}}Q_{lm}^{k}$ (11a) $\displaystyle=$
$\displaystyle\ w_{i}^{L}\sum_{k=1}^{L-1}P_{ij}^{k}\frac{\partial
F_{jl}^{k}}{\partial z^{k}_{p}}\frac{\partial z^{k}_{p}}{\partial
x^{0}_{n}}Q_{lm}^{k}$ (11b) $\displaystyle=$ $\displaystyle\
w_{i}^{L}\sum_{k=1}^{L-1}P_{ij}^{k}\frac{\partial F_{jl}^{k}}{\partial
z^{k}_{p}}Q_{pn}^{k}Q_{lm}^{k}$ (11c) $\displaystyle=$ $\displaystyle\
w_{i}^{L}\sum_{k=1}^{L-1}\sum_{j=1}^{d^{k}}P_{i(j)}^{k}s_{(j)}^{k}Q_{(j)m}^{k}Q_{(j)n}^{k}.$
(11d)

In the above derivation, eq. (11a) takes the derivative of _each_
$\mathbf{F}^{k}$ in eq. (4) with respect to $\mathbf{x}^{0}$ and sums them up,
showing the result in the spirit of eq. (10); eq. (11b) expands
$\frac{\partial\mathbf{F}^{k}}{\partial\mathbf{x}^{0}}$ by the chain rule,
i.e.,
$\frac{\partial\mathbf{F}^{k}}{\partial\mathbf{x}^{0}}=\frac{\partial\mathbf{F}^{k}}{\partial\mathbf{z}^{k}}\frac{\partial\mathbf{z}^{k}}{\partial\mathbf{x}^{0}}$;
eq. (11c) substitutes $\frac{\partial\mathbf{z}^{k}}{\partial\mathbf{x}^{0}}$
with $\mathbf{Q}^{k}$ as per its definition; and, finally, eq. (11d) takes
into account that $\frac{\partial F_{jl}^{k}}{\partial z^{k}_{p}}$ vanishes
unless $j=l=p$, as $\sigma^{k}$ is element-wise, whereas the non-zero terms
$s_{j}^{k}$ are the second derivatives of $\sigma^{k}(z)$:

$s_{j}^{k}:=\left.\frac{d}{dz}\frac{d\sigma^{k}(z)}{dz}\right|_{z=z^{k}_{j}}.$
(12)

For its later generalisation, we rewrite eq. (11) as

$\nabla\nabla\tilde{u}=\mathbf{w}^{L}\mathbf{V},$ (13)

where $\mathbf{V}$ is a third-order tensor defined by

$\displaystyle\mathbf{V}:=$ $\displaystyle\
\frac{\partial}{\partial\mathbf{x}^{0}}\frac{\partial\mathbf{x}^{L-1}}{\partial\mathbf{x}^{0}}\in\mathbb{R}^{d^{L-1}\times
d^{0}\times d^{0}},$ (14a) $\displaystyle V_{imn}=$ $\displaystyle\
\sum_{k=1}^{L-1}\sum_{j=1}^{d^{k}}P_{i(j)}^{k}s_{(j)}^{k}Q_{(j)m}^{k}Q_{(j)n}^{k}.$
(14b)

##  3 Incompatibility by ReLU

For a $\sigma^{k}$ of Lipschitz continuity like ReLU, we have
$\mathbf{s}^{k}\equiv\mathbf{0}$. It is then straightforward to see from eq.
(13) that the Hessian $\nabla\nabla\tilde{u}\equiv\mathbf{0}$ if the $(L-1)$
activation functions are _all_ ReLU-like. We provide a simple code to verify
this; see Section 5. Such a vanished Hessian as imposed by the neural network
is unwanted, which limits $\tilde{u}$ to a function space that excludes the
true solution $u$, or
$\mathcal{F}_{\mathrm{NN}}\cap\mathcal{F}_{\mathrm{PDE}}=\emptyset$,
regardless of the training data. From the perspective of Probably
Approximately Correct (PAC) learning [31],
$\mathcal{F}_{\mathrm{NN}}\cap\mathcal{F}_{\mathrm{PDE}}=\emptyset$ means the
breakdown of the “realisability assumption”, which fundamentally weakens a
model’s learnability. For better understanding, let us exam the 2D wave
equation on a membrane, considering the following three cases:

  1. 1.

if the wave equation is inhomogeneous, i.e., $u_{xx}+u_{yy}-u_{tt}/c^{2}=f$
with $f=f(x,y,t)$ being the source term, it can never be satisfied by a ReLU-
based $\tilde{u}$ because we have proved that
$\tilde{u}_{xx}=\tilde{u}_{yy}=\tilde{u}_{tt}\equiv 0$;

  2. 2.

if the wave equation is homogeneous, i.e., $u_{xx}+u_{yy}-u_{tt}/c^{2}=0$, it
will always be satisfied by a ReLU-based $\tilde{u}$; however, this illusory
benefit is invalid because the network-imposed conditions (not only
$\tilde{u}_{xx}=\tilde{u}_{yy}=\tilde{u}_{tt}\equiv 0$ but also
$\tilde{u}_{xy}=\tilde{u}_{xt}=\tilde{u}_{yt}\equiv 0$) are much stronger than
the wave equation itself;

  3. 3.

if the wave equation is both forced and damped, i.e,
$u_{xx}+u_{yy}-u_{tt}/c^{2}-\mu u_{t}/c^{2}=f$ (with $\mu$ being the
coefficient of friction), a ReLU-based MLP will effectively be trained to
minimise $-\mu\tilde{u}_{t}/c^{2}-f$, subject to
$\nabla\nabla\tilde{u}\equiv\mathbf{0}$, rather than to minimise
$\tilde{u}_{xx}+\tilde{u}_{yy}-\tilde{u}_{tt}/c^{2}-\mu\tilde{u}_{t}/c^{2}-f$
as wanted.

The vanished Hessian can be easily avoided by using a smooth activation
function, as many applications have done, with some being alert to the
downside of ReLU for PINNs [27, 28]. Another way is to use non-affine layer
operations, or non-affine $\mathcal{L}$’s in eq. (1), because
$(\sigma\mathcal{L})^{\prime\prime}$ will contain at least one non-zero term
${\sigma}^{\prime}\mathcal{L}^{\prime}$; an example of this is the Fourier
neural operators tailed by ReLU [11]. Note that convolutional filters are also
affine (with a sparse $\mathbf{W}^{k}$), so a ReLU-based convolutional neural
network (CNN) will also end up with a vanished Hessian.

He et al. [29] have studied the similarity between ReLU-based MLPs and the
finite element methods (FEMs) with a linear shape function [30]. They show
that both of them present a piecewise linear function space for the
approximation of the PDE solution, which seems to hint that a ReLU-based MLP
may achieve the same accuracy as a linear FEM does. We argue that such a
similarity or equivalence is only apparent because of the following two major
differences. First, in an FEM, the spatial domain is discretised by a mesh
whereby linear interpolation is localised in each (small) element with fixed
anchor points (i.e., the nodes). Such a local parameterisation with compact
support can fit any complicated functions within in the mesh resolution. In
contrast, a ReLU-based MLP presents a global piecewise linear parameterisation
whose anchor points are stochastically located, due to the stochastic nature
of neural weights, mostly uncontrollable a priori and difficult to learn from
data. Such a difference is visualised in Fig. 2. Second, an FEM is based on a
_weak form_ (or variational form) of the PDE [15] in which the second-order
derivatives have disappeared. However, such a weak form is not considered by a
vanilla PINN, which still computes the Hessian by AutoDiff as part of the loss
function (and it will only get zeros if ReLU is used). The above two
differences are most relevant to the idea of PINNs with domain decomposition
[8, 15], that is, many smaller neural networks are trained in parallel, each
approximating the solution in a smaller subdomain, either using a weak form
[15] or a strong form [8] of the PDE. Domain decomposition can be understood
as to achieve controlled localisation of basis functions by using many neural
networks, similar to an FEM by using many elements. In addition to the weak
form, an alternative way to avoid computing the Hessian is to reformulate a
second-order PDE as a system of first-order PDEs [32], similar to the
staggered-grid finite difference methods [30].

$x$$u(x)$True solutionFinite-element solutionMesh-allowed
solutions$x_{1}$$x_{2}$$x_{3}$$x_{4}$$x_{5}$

(a) FEM

$x$$u(x)$True solutionNetwork-allowed solutions$x_{1}$$x_{5}$

(b) ReLU-based MLP

Figure 2: Sketch of 1D piecewise linear parameterisations (a) by an FEM and
(b) by a ReLU-based MLP. The main difference between (a) and (b) is that the
anchor points are predefined in (a), i.e., $x_{1},x_{2},\cdots,x_{5}$, based
on the required mesh resolution, but are stochastic and uncontrollable in (b).
A simple code is provided to produce the sketched patterns in (b) with
randomly initialised MLPs; see Section 5.

##  4 The out-layer-hyperplane

In the previous section, we have shown that a ReLU-based MLP is incompatible
with a second-order PDE, or
$\mathcal{F}_{\mathrm{NN}}\cap\mathcal{F}_{\mathrm{PDE}}=\emptyset$. The above
process has delightfully inspired a simple and elegant way to enforce an MLP
to satisfy any second-order linear PDEs, i.e.,
$\mathcal{F}_{\mathrm{NN}}\in\mathcal{F}_{\mathrm{PDE}}$, without imposing any
non-physical constraints. Now we assume that the activation functions are of
$C^{2}$ or above so that $\mathbf{V}$ in eq. (14) does not vanish.

Consider a general second-order, scalar-valued linear PDE as follow:

$\Gamma_{mn}u_{,mn}+\gamma_{m}u_{,m}+\beta u+\alpha=0,$ (15)

where $\alpha\in\mathbb{R}$, $\beta\in\mathbb{R}$,
$\bm{\upgamma}\in\mathbb{R}^{d^{0}}$ and
$\mathbf{\Gamma}\in\mathbb{R}^{d^{0}\times d^{0}}$ are known PDE coefficients.
Let $u(\mathbf{x}^{0})$ be approximated by $\tilde{u}(\mathbf{x}^{0})$.
Substituting $u$ by $z^{L}$ in eq. (2), $u_{,m}$ by eq. (6) and $u_{,mn}$ by
eq. (13), we obtain the following linear system that enforces the MLP to
satisfy the above PDE:

$\psi_{i}w_{i}^{L}+\beta b^{L}+\alpha=0,$ (16)

where $\bm{\uppsi}\in\mathbb{R}^{d^{L-1}}$ are defined by

$\psi_{i}:=\Gamma_{mn}V_{imn}+\gamma_{m}U_{im}+\beta x_{i}^{L-1}.$ (17)

The LHS of eq. (16) is a function of all the weights and biases of the MLP and
the input $\mathbf{x}^{0}$. However, from the angle of network design, we can
regard the weights and biases of the hidden layers ($\mathbf{w}^{k}$ and
$b^{k}$, $k\in\\{1,2,\cdots,L-1\\}$) as free variables or trainable parameters
while constraining those of the output layer ($\mathbf{w}^{L}$ and $b^{L}$) so
that eq. (16) can always hold. Obviously,
$\mathbf{w}^{L}\in\mathbb{R}^{d^{L-1}}$ and $b^{L}\in\mathbb{R}$ lie on a
hyperplane in $\mathbb{R}^{d^{L-1}+1}$, as defined by eq. (16), which we call
the _out-layer-hyperplane_.

We realise the out-layer-hyperplane using its parametric equation. Without
loss of generality, we ignore $b^{L}$ by assuming $\beta=0$ for a compact
notation. Known from basic linear algebra, the parametric equation of the
hyperplane is given by

$\mathbf{w}^{L}=\mathbf{g}^{[1]}+\sum_{p=2}^{d^{L-1}}\lambda^{[p]}\mathbf{g}^{[p]},$
(18)

where the basis vectors $\mathbf{g}^{[p]}\in\mathbb{R}^{d^{L-1}}$,
$p\in\\{1,2,\cdots,d^{L-1}\\}$, can be chosen as

$\left[\begin{array}[]{c}\mathbf{g}^{[1]}\\\
\hdashline[2pt/2pt]\mathbf{g}^{[2]}\\\ \mathbf{g}^{[3]}\\\ \vdots\\\
\mathbf{g}^{[d^{L-1}]}\end{array}\right]=\left[\begin{array}[]{c;{2pt/2pt}cccccc}\frac{-\alpha\psi_{1}}{|\bm{\uppsi}|^{2}}&\frac{-\alpha\psi_{2}}{|\bm{\uppsi}|^{2}}&\frac{-\alpha\psi_{3}}{|\bm{\uppsi}|^{2}}&\cdots&\frac{-\alpha\psi_{d^{L-1}}}{|\bm{\uppsi}|^{2}}\\\
\hdashline[2pt/2pt]-\psi_{2}&\psi_1&0&\cdots&0\\\
-\psi_{3}&0&\psi_1&\cdots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
-\psi_{d^{L-1}}&0&0&\cdots&\psi_{1}\end{array}\right],$ (19)

and $\lambda^{[p]}\in\mathbb{R}$, $p\in\\{2,3,\cdots,d^{L-1}\\}$ (note that
$p$ starts from 2 here), are the new network parameters taking the place of
$\mathbf{w}^{L}$. Here we put the superscript indices for $\lambda$ and
$\mathbf{g}$ in $[\cdot]$ to differentiate them with the layer indices. Note
that $\mathbf{g}^{[1]}$ represents a point on the hyperplane, and
$\mathbf{g}^{[2]},\mathbf{g}^{[3]},\cdots,\mathbf{g}^{[d^{L-1}]}$ are linearly
independent vectors parallel to the hyperplane. The choice of these basis
vectors is non-unique, such as
$\mathbf{g}^{[1]}\leftarrow\\{-\alpha/\psi_{1},0,0,\cdots,0\\}$, but our
choice in eq. (19) avoids division by elements of $\bm{\uppsi}$ (which may be
close to zero), while the division by $|\bm{\uppsi}|^{2}$ is always safe and
stable. Also note that the total number of the $\lambda$’s is one less than
the size of $\mathbf{w}^{L}$, as $\mathbf{w}^{L}$ must lie on a hyperplane,
which certifies that we do not impose any extra conditions other than the PDE
itself. The forward pass through an MLP equipped with the out-layer-hyperplane
is sketched in Fig. 3 and elaborated in Algorithm 1. A code implementation is
also provided; see Section 5.

Thus far, we have shown the existence and a realisation of the out-layer-
hyperplane for second-order linear PDEs approximated by an MLP. What makes
this idea more attractive is that it can be generalised to other network
architectures and PDE classes, most of which require a trivial effort, as
detailed below.

  * 1.

Network architecture. Equation (16) holds for any network architecture tailed
by a fully-connected hidden layer as long as we relax the definitions of
$\mathbf{U}$ and $\mathbf{V}$ in eqs. (9) and (14) respectively to
$\mathbf{U}:=\frac{\partial\mathbf{x}^{L-1}}{\partial\mathbf{x}^{0}}$ and
$\mathbf{V}:=\frac{\partial}{\partial\mathbf{x}^{0}}\frac{\partial\mathbf{x}^{L-1}}{\partial\mathbf{x}^{0}}$
(i.e., abandon eq. (14b)). For an MLP, we have given their closed-form
expressions in eqs. (9) and (14), corresponding to Steps 2$\sim$4 in Algorithm
1. For a general network architecture, we can simply compute $\mathbf{U}$ and
$\mathbf{V}$ by AutoDiff using the first $(L-1)$ layers based on their
generalised definitions.

  * 2.

Higher-order linear PDEs. It is straightforward to show that the out-layer-
hyperplane exists for an $n$-th order linear PDE given $C^{n}$ activation
functions. For example, a general third-order linear PDE will require an extra
term $\Omega_{mnr}u_{,mnr}$ in eq. (15), with $\Omega_{mnr}$ being known PDE
coefficients; this term will add $\Omega_{mnr}T_{imnr}$ to the hyperplane
coefficients ${\psi_{i}}$ in eq. (17), where
$\mathbf{T}:=\frac{\partial}{\partial\mathbf{x}^{0}}\frac{\partial}{\partial\mathbf{x}^{0}}\frac{\partial\mathbf{x}^{L-1}}{\partial\mathbf{x}^{0}}$,
and $\mathbf{T}\neq\mathbf{0}$ with $C^{3}$ or above activation functions.

  * 3.

Vector-valued linear PDEs. Let us consider a general vector-valued, second-
order linear PDE with solution $\mathbf{u}(\mathbf{x}^{0})\in\mathbb{R}^{3}$,
formulated as

$\Theta_{ijmn}u_{j,mn}+\Pi_{ijm}u_{j,m}+\Lambda_{ij}u_{j}+\alpha_{i}=0,\quad
i\in\\{1,2,3\\},$ (20)

which is the $\mathbb{R}^{3}$ version of eq. (15). Consequently, one can show
that eq. (16) will be generalised to

$\Psi_{ijk}W_{jk}^{L}+\Lambda_{ij}b^{L}_{j}+\alpha_{i}=0,\quad
i\in\\{1,2,3\\},$ (21)

where
$\Psi_{ijk}=\Theta_{ijmn}V_{kmn}+\Pi_{ijm}U_{km}+\Lambda_{ij}x_{k}^{L-1}$,
with $\mathbf{U}$ and $\mathbf{V}$ remaining the same as in eqs. (9) and (14).
Evidently, the above linear system states that the PDE can be enforced if the
weights and biases of the output layer, $\mathbf{W}^{L}\in\mathbb{R}^{3\times
d^{L-1}}$ and $\mathbf{b}^{L}\in\mathbb{R}^{3}$, lie on three coupled
hyperplanes, which can be realised by parametric equations similar to eq.
(18).

  * 4.

Nonlinear PDEs. It is not difficult to guess that, to satisfy a nonlinear PDE,
the weights of the output layer must lie on a hypersurface. This is true. For
example, consider $(u_{,i}u)$, a common nonlinear term that appears in many
PDEs such as Navier–Stokes and Bateman-Burgers [33]. It will introduce many
quadratic terms about $\mathbf{w}^{L}$ into eq. (16), namely,
$A_{ijk}w_{j}^{L}w_{k}^{L}$, where $A_{ijk}=U_{ij}x_{k}^{L-1}$. For a general
nonlinear PDE, there is no universal way to parameterise its out-layer-
hypersurface. For certain PDEs, however, a closed-form parameterisation
similar to eq. (18) may be found. At least, concerning second- and third-order
nonlinear PDEs, the parametric equations for the generalised quadratic and
cubic hypersurfaces have been well established [34].

[height=4.5, nodespacing=14mm, nodesize=25pt, layerspacing=23mm] [count=3,
bias=false, title=Input, text=$x^{]}_{\hiddenlayer}$[count=3, bias=false,
title=First
hidden
layer, text=$x^{]}_{\linklayers}$[count=4, bias=false, title=Last
hidden
layer, text=$x^{]}_{\linklayers}$[count=1, title=Output
layer, text=$\tilde{u}$[thick, ->] (7.2,-2.9) – (8.2,-1.3); [align=center] at
(9.1, -1.0) Data loss ;

[thick, ->] (7.5,-3.3) – (8.1,-3.3); [align=center] at (9.1, -4.9) AutoDiff;
[fill=blue!35, rounded corners=.4cm, align=center, inner sep=.2cm] at (9.1,
-3.3) $\dfrac{\partial\tilde{u}}{\partial\mathbf{x}^{0}}$,

$\dfrac{\partial}{\partial\mathbf{x}^{0}}\dfrac{\partial\tilde{u}}{\partial\mathbf{x}^{0}}$,

$\cdots$ ; [thick, ->] (10.1,-3.3) – (10.7,-3.3);

[from layer=2, from node=1, to layer=3, to node=1, label=$w^{2}_{\link}$rom
layer=2, from node=2, to layer=3, to node=1, label=$w^{2}_{\link}$rom layer=2,
from node=3, to layer=3, to node=1, label=$w^{2}_{\link}$rom layer=2, from
node=4, to layer=3, to node=1, label=$w^{2}_{\node}$lign=center] at (11.5,
-3.3) ~~PDE loss~~
IC loss
BC loss ;

[align=center] at (2.3,-6.2) $\underbrace{\quad\quad\quad}$; [fill=pink,
rounded rectangle, rounded rectangle east arc=none, align=center, anchor=east]
at (2.35,-6.9) $\mathbf{W}^{1}$; [fill=blue!35, rounded rectangle, rounded
rectangle west arc=none, align=center, anchor=west] at (2.3,-6.9)
$\,\mathbf{z}^{1}$; [align=center] at (2.3,-7.5) $\downarrow$;

[align=center] at (4.6,-6.2) $\underbrace{\quad\quad\quad}$; [fill=pink,
rounded rectangle, rounded rectangle east arc=none, align=center, anchor=east]
at (4.65,-6.9) $\mathbf{W}^{2}$; [fill=blue!35, rounded rectangle, rounded
rectangle west arc=none, align=center, anchor=west] at (4.6,-6.9)
$\,\mathbf{z}^{2}$; [align=center] at (4.6,-7.5) $\downarrow$;

[fill=green!40, rounded rectangle, align=center, inner sep=0.2cm] at
(-.1,-8.2) PDE
coefficients; [align=center] at (1.2,-8.2) $\rightarrow$; [draw, align=center,
inner sep=.2cm] at (3.5,-8.2) Steps 1$\sim$5 in Algorithm 1; [align=center] at
(5.8,-8.2) $\rightarrow$;

[fill=blue!35, rounded rectangle, align=center] at (7.6,-8.2)
$\mathbf{g}^{[1]},\mathbf{g}^{[2]},\mathbf{g}^{[3]},\mathbf{g}^{[4]}$;

[fill=pink, rounded rectangle, align=center, inner sep=.2cm] at (7.8,-6.9)
Trainables
$\lambda^{[2]},\lambda^{[3]},\lambda^{[4]}$;

[align=center] at (9.4,-8.1) $\nearrow$;

[align=center] at (9.4,-7.1) $\searrow$;

[fill=orange!50, rounded rectangle,inner sep=.2cm, align=center] at
(10.3,-7.6) $\mathbf{w}^{2}$;

[->, thick, orange] (10.3,-7.2) to [out=70,in=-70] (5.8,-4.8);

Figure 3: Forward pass through an MLP ($L=3$) equipped with the out-layer-
hyperplane. In an ordinary MLP, the four elements of $\mathbf{w}^{2}$
(coloured in orange) are trained. Using the out-layer-hyperplane, the three
$\lambda$’s (coloured in pink) become the new trainable parameters, which
always yield a $\mathbf{w}^{2}$ lying on the out-layer-hyperplane so that
$\tilde{u}$ always satisfies the target linear PDE.  Algorithm 1 Forward pass
through an MLP equipped with the out-layer-hyperplane.

Input: $\mathbf{x}^{0}\in\mathbb{R}^{d^{0}}$.
Output: $\tilde{u}\in\mathbb{R}$.
Network parameters (for backprop): $\mathbf{W}^{k}\in\mathbb{R}^{d^{k}\times
d^{k-1}}$ and $\mathbf{b}^{k}\in\mathbb{R}^{d^{k}}$,
$k\in\\{1,2,\cdots,L-1\\}$; $\lambda^{[p]}\in\mathbb{R}$,
$p\in\\{2,3,\cdots,d^{L-1}\\}$.
PDE coefficients: $\alpha\in\mathbb{R}$, $\beta\in\mathbb{R}$,
$\bm{\upgamma}\in\mathbb{R}^{d^{0}}$ and
$\mathbf{\Gamma}\in\mathbb{R}^{d^{0}\times d^{0}}$.

1:Perform ordinary forward pass until the output layer, obtaining
$\mathbf{x}^{L-1}\in\mathbb{R}^{d^{L-1}}$ and
$\mathbf{z}^{k}\in\mathbb{R}^{d^{k}}$, $k\in\\{1,2,\cdots,L-1\\}$; #
$\mathbf{z}^{k}$’s are usually unsaved in an ordinary MLP, but we need them to
compute the first and second derivatives of the activation functions;

2:Compute $\mathbf{F}^{k}\in\mathbb{R}^{d^{k}\times d^{k}}$ and
$\mathbf{s}^{k}\in\mathbb{R}^{d^{k}}$, $k\in\\{1,2,\cdots,L-1\\}$ by eqs. (5)
and (12);

3:Compute $\mathbf{P}^{k}\in\mathbb{R}^{d^{L-1}\times d^{k}}$ and
$\mathbf{Q}^{k}\in\mathbb{R}^{d^{k}\times d^{0}}$ $k\in\\{1,2,\cdots,L-1\\}$
by eqs. (7) and (8);

4:Compute $\mathbf{U}\in\mathbb{R}^{d^{L-1}\times d^{0}}$ and
$\mathbf{V}\in\mathbb{R}^{d^{L-1}\times d^{0}\times d^{0}}$ by eqs. (9) and
(14);

5:Compute $\bm{\uppsi}\in\mathbb{R}^{d^{L-1}}$ by eq.(17) and then
$\mathbf{g}^{[p]}\in\mathbb{R}^{d^{L-1}}$, $p\in\\{1,2,\cdots,d^{L-1}\\}$ by
(19); # this is where the PDE coefficients are used;

6:Reconstruct $\mathbf{w}^{L}\in\mathbb{R}^{d^{L-1}}$ by (18);  # this is
where the network parameters $\lambda^{[p]}$’s are used;

7:Return $\tilde{u}=\mathbf{w}^{L}\mathbf{x}^{L-1}+b^{L}$.

##  5 Validation by code

To support some of the key statements in this paper, we provide a lightweight
code repository (https://github.com/stfc-sciml/PINN-PDE-compatibility) from
which the following three standalone PyTorch [35] scripts can be found:

  * 1.

relu_causes_zero_hessian.py verifies that a ReLU-based MLP will always lead to
a vanished Hessian;

  * 2.

piecewise_linear_by_relu.py plots the piecewise linear functions generated by
some random ReLU-based MLPs, as a support to our sketch in Figure 2(b);

  * 3.

out_layer_hyperplane.py implements the out-layer-hyperplane in an MLP for a 2D
inhomogeneous wave equation on a membrane:
$u_{xx}+u_{yy}-u_{tt}-u_{t}=\sin(x+y-t)$, following Algorithm 1.

##  6 Conclusions

In physics-informed learning, it is crucial that the function space granted by
a neural network, $\mathcal{F}_{\mathrm{NN}}$, is compatible with that
required by the target PDE, $\mathcal{F}_{\mathrm{PDE}}$. We have proved that
an MLP with only ReLU-like activation functions will always lead to a vanished
Hessian, so it is incompatible with any second- or higher-order PDEs, namely,
$\mathcal{F}_{\mathrm{NN}}\cap\mathcal{F}_{\mathrm{PDE}}=\emptyset$.
Therefore, we proclaim that ReLU should not be used in a vanilla PINN. This
has led us to a simple method to make $\mathcal{F}_{\mathrm{NN}}$ fully
compatible with $\mathcal{F}_{\mathrm{PDE}}$, namely,
$\mathcal{F}_{\mathrm{NN}}\in\mathcal{F}_{\mathrm{PDE}}$. We have proved that,
using sufficiently smooth activation functions, a neural network tailed by a
fully-connected hidden layer can strictly satisfy any high-order linear PDE
when the weights of its output layer lie on a certain hyperplane, as called
the out-layer-hyperplane. A closed-form expression of this hyperplane has been
given for second-order linear PDEs. To the best of our knowledge, this should
be the first PINN architecture that enforces point-wise correctness of PDEs.
Our future work will be focused on the efficacy and performance of the out-
layer-hyperplane for different types of architectures and PDEs.

## Acknowledgement

This work is supported by the EPSRC grant, Blueprinting for AI for Science at
Exascale (BASE-II, EP/X019918/1), which is Phase II of the Benchmarking for AI
for Science at Exascale (BASE) grant.

## References

  * [1] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics 378 (2019) 686–707.
  * [2] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, Physics-informed machine learning, Nature Reviews Physics 3 (6) (2021) 422–440.
  * [3] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, F. Piccialli, Scientific machine learning through physics-informed neural networks: Where we are and what’s next, Journal of Scientific Computing 92 (2022) 88.
  * [4] L. Lu, X. Meng, Z. Mao, G. E. Karniadakis, Deepxde: A deep learning library for solving differential equations, SIAM Review 63 (1) (2021) 208–228.
  * [5] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, A. Stuart, K. Bhattacharya, A. Anandkumar, Multipole graph neural operator for parametric partial differential equations, Advances in Neural Information Processing Systems 33 (2020) 6755–6766.
  * [6] L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis, Learning nonlinear operators via deeponet based on the universal approximation theorem of operators, Nature Machine Intelligence 3 (3) (2021) 218–229.
  * [7] X. Meng, G. E. Karniadakis, A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse pde problems, Journal of Computational Physics 401 (2020) 109020.
  * [8] B. Moseley, A. Markham, T. Nissen-Meyer, Finite basis physics-informed neural networks (fbpinns): a scalable domain decomposition approach for solving differential equations, arXiv preprint arXiv:2107.07871 (2021).
  * [9] L. Yang, X. Meng, G. E. Karniadakis, B-pinns: Bayesian physics-informed neural networks for forward and inverse pde problems with noisy data, Journal of Computational Physics 425 (2021) 109913.
  * [10] M. Kasim, D. Watson-Parris, L. Deaconu, S. Oliver, P. Hatfield, D. Froula, G. Gregori, M. Jarvis, S. Khatiwala, J. Korenaga, et al., Building high accuracy emulators for scientific simulations with deep neural architecture search, Machine Learning: Science and Technology 3 (1) (2021) 015013.
  * [11] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, A. Anandkumar, Fourier neural operator for parametric partial differential equations, arXiv preprint arXiv:2010.08895 (2020).
  * [12] M. Tancik, P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, J. Barron, R. Ng, Fourier features let networks learn high frequency functions in low dimensional domains, Advances in Neural Information Processing Systems 33 (2020) 7537–7547.
  * [13] Ziqi, Liu, , 9594, , Z. Liu, Wei, Cai, , 9595, , W. Cai, Zhi-Qin, J. Xu, , 9596, , Z.-Q. J. Xu, Multi-scale deep neural network (mscalednn) for solving poisson-boltzmann equation in complex domains, Communications in Computational Physics 28 (5) (2020) 1970–2001.
  * [14] W. Cai, X. Li, L. Liu, A phase shift deep neural network for high frequency approximation and wave problems, SIAM Journal on Scientific Computing 42 (5) (2020) A3285–A3312.
  * [15] E. Kharazmi, Z. Zhang, G. E. Karniadakis, hp-vpinns: Variational physics-informed neural networks with domain decomposition, Computer Methods in Applied Mechanics and Engineering 374 (2021) 113547.
  * [16] B. Ronen, D. Jacobs, Y. Kasten, S. Kritchman, The convergence rate of neural networks for learned functions of different frequencies, Advances in Neural Information Processing Systems 32 (2019).
  * [17] A. D. Jagtap, K. Kawaguchi, G. Em Karniadakis, Locally adaptive activation functions with slope recovery for deep and physics-informed neural networks, Proceedings of the Royal Society A 476 (2239) (2020) 20200334.
  * [18] Z. Mao, A. D. Jagtap, G. E. Karniadakis, Physics-informed neural networks for high-speed flows, Computer Methods in Applied Mechanics and Engineering 360 (2020) 112789.
  * [19] J. Yu, L. Lu, X. Meng, G. E. Karniadakis, Gradient-enhanced physics-informed neural networks for forward and inverse pde problems, Computer Methods in Applied Mechanics and Engineering 393 (2022) 114823.
  * [20] S. Wang, X. Yu, P. Perdikaris, When and why pinns fail to train: A neural tangent kernel perspective, Journal of Computational Physics 449 (2022) 110768\.
  * [21] S. Mishra, R. Molinaro, Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for pdes, IMA Journal of Numerical Analysis 42 (2) (2022) 981–1022.
  * [22] H. Sheng, C. Yang, Pfnn: A penalty-free neural network method for solving a class of second-order boundary-value problems on complex geometries, Journal of Computational Physics 428 (2021) 110085.
  * [23] P. L. Lagari, L. H. Tsoukalas, S. Safarkhani, I. E. Lagaris, Systematic construction of neural forms for solving partial differential equations inside rectangular domains, subject to initial, boundary and interface conditions, International Journal on Artificial Intelligence Tools 29 (05) (2020) 2050009.
  * [24] S. Dong, N. Ni, A method for representing periodic functions and enforcing exactly periodic boundary conditions with deep neural networks, Journal of Computational Physics 435 (2021) 110242.
  * [25] M. Mattheakis, P. Protopapas, D. Sondak, M. Di Giovanni, E. Kaxiras, Physical symmetries embedded in neural networks, arXiv preprint arXiv:1904.08991 (2019).
  * [26] J. Ling, A. Kurzawski, J. Templeton, Reynolds averaged turbulence modelling using deep neural networks with embedded invariance, Journal of Fluid Mechanics 807 (2016) 155–166.
  * [27] S. Markidis, The old and the new: Can physics-informed deep-learning replace traditional linear solvers?, Frontiers in big Data (2021) 92.
  * [28] B. Moseley, A. Markham, T. Nissen-Meyer, Solving the wave equation with physics-informed deep learning, arXiv preprint arXiv:2006.11894 (2020).
  * [29] J. He, L. Li, J. Xu, C. Zheng, Relu deep neural networks and linear finite elements, Journal of Computational Mathematics 38 (3) (2020) 502–527.
  * [30] H. Igel, Computational seismology: a practical introduction, Oxford University Press, 2017.
  * [31] S. Shalev-Shwartz, S. Ben-David, Understanding machine learning: From theory to algorithms, Cambridge university press, 2014.
  * [32] R. J. Gladstone, M. A. Nabian, H. Meidani, Fo-pinns: A first-order formulation for physics informed neural networks, arXiv preprint arXiv:2210.14320 (2022).
  * [33] L. Debnath, L. Debnath, Nonlinear partial differential equations for scientists and engineers, Springer, 2005.
  * [34] C. Bajaj, Quadric and cubic hypersurface parameterization, Tech. rep., Purdue University (1989).
  * [35] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al., Pytorch: An imperative style, high-performance deep learning library, Advances in neural information processing systems 32 (2019).