# Arithmetic properties of the Genocchi numbers and their generalization

Bakir FARHI
Laboratoire de Mathématiques appliquées
Faculté des Sciences Exactes
Université de Bejaia, 06000 Bejaia, Algeria
<EMAIL_ADDRESS>
http://farhi.bakir.free.fr/

###### Abstract

This note is devoted to establish some new arithmetic properties of the
generalized Genocchi numbers $G_{n,a}$ ($n\in\mathbb{N}$, $a\geq 2$). The
resulting properties for the usual Genocchi numbers $G_{n}=G_{n,2}$ are then
derived. We show for example that for any even positive integer $n$, the
Genocchi number $G_{n}$ is a multiple of the odd part of $n$.

MSC 2010: Primary 11B68, 13F25, 11A07.
Keywords: Genocchi numbers, generalized Genocchi numbers, Bernoulli numbers,
formal power series, IDC-series, congruences.

##  1 Introduction and Notations

Throughout this note, we let ${\mathbb{N}}^{*}$ denote the set of positive
integers. For a given prime number $p$, we let $\vartheta_{p}$ denote the
usual $p$-adic valuation. For a given rational number $r$, we let
$\mathrm{num}(r)$ and $\mathrm{den}(r)$ respectively denote the numerator and
the denominator of $r$; precisely, if $d$ is the smallest positive integer
satisfying $dr\in{\mathbb{Z}}$ then we have: $\mathrm{den}(r):=d$ and
$\mathrm{num}(r):=dr$. Next, for given positive integers $a$ and $n$, we let
$\pi_{a}(n)$ denote the greatest positive divisor of $n$ which is coprime with
$a$. Then, it is immediate that:

$\pi_{a}(n)~{}=~{}\prod_{\begin{subarray}{c}p\text{ prime}\\\ p\nmid
a\end{subarray}}p^{\vartheta_{p}(n)}.$

In the particular case $a=2$, $\pi_{a}(n)$ is called _the odd part_ of $n$; it
is the greatest odd positive divisor of $n$.

Additionally, we need to extend the congruences in ${\mathbb{Z}}$ to
${\mathbb{Q}}$ as follows:

> For $a,b\in{\mathbb{Q}}$ and $n\in{\mathbb{N}}^{*}$, we write $a\equiv b\
> (\mathrm{mod}\ n)$ if the numerator of the rational number $(a-b)$ is a
> multiple of $n$.

In the same context, an equivalent meaning of the congruence $a\equiv b\
(\mathrm{mod}\ n)$ consists to say that for any prime divisor $p$ of $n$, we
have $\vartheta_{p}(a-b)\geq\vartheta_{p}(n)$. We can easily check that some
properties (but not all) of the congruences in ${\mathbb{Z}}$ become valid in
${\mathbb{Q}}$. For example, we can sum side to side several congruences in
${\mathbb{Q}}$ with a same modulus. Also, if $a,b,c\in{\mathbb{Q}}$ and
$n\in{\mathbb{N}}^{*}$ such that $\mathrm{den}(c)$ is coprime with $n$ then we
have:

$a\equiv b\ (\mathrm{mod}\ n)~{}\Longrightarrow~{}ac\equiv bc\ (\mathrm{mod}\
n).$

However, we cannot multiply (side to side) several congruences in
${\mathbb{Q}}$ with a same modulus (contrary to the congruences in
${\mathbb{Z}}$).

Further, we denote by ${\mathbb{Q}}[[t]]$ the ring of formal power series in
$t$ with coefficients in ${\mathbb{Q}}$. If an element $f$ of
${\mathbb{Q}}[[t]]$ is represented as

$f(t)~{}=~{}\sum_{n=0}^{+\infty}a_{n}\frac{t^{n}}{n!}$

(where $a_{n}\in{\mathbb{Q}}$, $\forall n\in{\mathbb{N}}$), we call the
$a_{n}$’s _the differential coefficients_ of $f$ (because each $a_{n}$ is the
$n$th derivative of $f$ at $0$). If the $a_{n}$’s are all integers, we say
that $f$ is an _IDC-series_ (IDC abbreviates the expression “with Integral
Differential Coefficients”). Many usual functions are IDC-series; we can cite
for example the functions $x\mapsto e^{x}$, $x\mapsto\sin{x}$,
$x\mapsto\cos{x}$, $x\mapsto\ln(1+x)$, and so on. We can easily see that the
sum and the product of two IDC-series become an IDC-series. Other important
properties of the IDC-series are given by the author [3] who used them to
establish his generalization of the Genocchi theorem (see below). One of those
properties is given by the following proposition:

######  Proposition 1.1 ([3, Corollary 2.2]).

Let $f(t)=\sum_{n=0}^{+\infty}a_{n}\frac{t^{n}}{n!}$ be an IDC-series with
$a_{0}\neq 0$. Then the formal power series $\frac{a_{0}}{f(a_{0}t)}$ is also
an IDC-series.

Showing that a given function is an IDC-series is, in general, not easy. For
example, the statement that the function $t\mapsto\frac{2t}{e^{t}+1}$ is an
IDC-series constitute a profound theorem due to Genocchi [5]. The differential
coefficients of the last function (which thus are all integers) are called
_the Genocchi numbers_ and denoted by $G_{n}$ ($n\in{\mathbb{N}}$); so we
have:

$\frac{2t}{e^{t}+1}~{}=~{}\sum_{n=0}^{+\infty}G_{n}\frac{t^{n}}{n!}.$

For a study of the Genocchi numbers, the reader can consult [1, 2, 3, 4, 5,
8].

Very recently, the author [3] has generalized the Genocchi numbers by
considering, for a given integer $a\geq 2$, the function
$t\mapsto\frac{at}{e^{(a-1)t}+e^{(a-2)t}+\dots+e^{t}+1}$ and its expansion
into a power series:

$\frac{at}{e^{(a-1)t}+e^{(a-2)t}+\dots+e^{t}+1}~{}=~{}\sum_{n=0}^{+\infty}G_{n,a}\frac{t^{n}}{n!}.$

So, for $a=2$, we simply obtain the usual Genocchi numbers; that is
$G_{n,2}=G_{n}$ ($\forall n\in{\mathbb{N}}$). The main result of [3] then
states that the $G_{n,a}$’s are all integers, which generalizes the Genocchi
theorem. The Genocchi numbers as well as their generalization are ultimately
related to the Bernoulli numbers ${(B_{n})}_{n\in{\mathbb{N}}}$, which can be
defined by the exponential generating function:

$\frac{t}{e^{t}-1}~{}=~{}\sum_{n=0}^{+\infty}B_{n}\frac{t^{n}}{n!}$

(see e.g., [1, 7]). Among others, we have the following connection between the
generalized Genocchi numbers and the Bernoulli numbers:

######  Proposition 1.2 ([3, Proposition 4.1]).

For any positive integers $a$ and $n$, with $a\geq 2$, we have:

$G_{n,a}~{}=~{}\sum_{k=0}^{n-1}\binom{n}{k}B_{k}a^{k}.$

Furthermore, the Bernoulli numbers have many arithmetic properties. The most
famous is certainly that given by the von Staudt-Clausen theorem (see e.g.,
[6, 7]), recalled below:

######  Theorem 1.3 (von Staudt-Clausen).

For any even positive integer $n$, we have that:

$B_{n}+\sum_{\begin{subarray}{c}p\text{ prime}\\\ (p-1)\mid
n\end{subarray}}\frac{1}{p}~{}\in~{}{\mathbb{Z}}.$

Because it is known that $B_{0}=1$, $B_{1}=-\frac{1}{2}$ and $B_{n}=0$ for any
odd integer $n\geq 3$, we can derive from the von Staudt-Clausen theorem the
following immediate corollary which is useful for our purpose:

######  Corollary 1.4.

For any natural number $n$ and any prime number $p$, we have that:

$\vartheta_{p}(B_{n})\geq-1.$

The purpose of this note consists to establish arithmetic properties of the
numbers $G_{n,a}$ (in addition to their integrality, previously proven in
[3]). Precisely, we obtain for the $G_{n,a}$’s a nontrivial simple divisor and
a nontrivial congruence modulo $a$.

##  2 The results and the proofs

We begin by establishing nontrivial divisors for the generalized Genocchi
numbers.

######  Theorem 2.1.

Let $a\geq 2$ be an integer. Then for all positive integer $n$, the integer
$G_{n,a}$ is a multiple of the integer $\pi_{a}(n)$.

###### Proof.

By applying Proposition 1.1 to the IDC-series $t\mapsto
e^{(a-1)t}+e^{(a-2)t}+\dots+e^{t}+1$, we find that the function
$t\stackrel{{\scriptstyle
f}}{{\mapsto}}\frac{a}{e^{(a-1)at}+e^{(a-2)at}+\dots+e^{at}+1}$ is an IDC-
series. Now, let us explicitly express the differential coefficients of $f$ in
terms of the generalized Genocchi numbers. We have:

$\displaystyle\frac{a}{e^{(a-1)t}+e^{(a-2)t}+\dots+e^{t}+1}$ $\displaystyle=$
$\displaystyle\frac{1}{t}\cdot\frac{at}{e^{(a-1)t}+e^{(a-2)t}+\dots+e^{t}+1}$
$\displaystyle=$
$\displaystyle\frac{1}{t}\sum_{n=1}^{+\infty}G_{n,a}\frac{t^{n}}{n!}$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{+\infty}\frac{G_{n,a}}{n}\cdot\frac{t^{n-1}}{(n-1)!}.$

Then, by substituting $t$ by $at$, we get

$f(t)~{}=~{}\sum_{n=1}^{+\infty}\frac{a^{n-1}G_{n,a}}{n}\cdot\frac{t^{n-1}}{(n-1)!}.$

So, since $f$ is an IDC-series, we derive that
$\frac{a^{n-1}G_{n,a}}{n}\in{\mathbb{Z}}$ ($\forall n\in{\mathbb{N}}^{*}$);
that is

$n\mid a^{n-1}G_{n,a}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(\forall
n\in{\mathbb{N}}^{*}).$

Next, since $\pi_{a}(n)$ is a divisor of $n$ which is coprime with $a$, we
derive that $\pi_{a}(n)\mid a^{n-1}G_{n,a}$ and $\pi_{a}(n)$ is coprime with
$a^{n-1}$. Then, by Gauss’s lemma, we conclude that $\pi_{a}(n)$ divides
$G_{n,a}$, as required. ∎

By taking $a=2$ in Theorem 2.1, we immediately deduce the following corollary
which concerns the usual Genocchi numbers:

######  Corollary 2.2.

For all positive integer $n$, the Genocchi number $G_{n}$ is a multiple of the
odd part of $n$. $\square$

We now turn to study the remainders of the generalized Genocchi numbers
$G_{n,a}$ modulo $a$. The main result in this direction is the following:

######  Theorem 2.3.

For any integers $a$ and $n$, both greater than $1$, we have

$G_{n,a}~{}\equiv~{}1-\frac{n}{2}\,a\ (\mathrm{mod}\ a).$

###### Proof.

Let $a$ and $n$ be two integers, both greater than $1$. According to
Proposition 1.2, we have that:

$\displaystyle G_{n,a}$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{n-1}\binom{n}{k}B_{k}a^{k}$ (2.1) $\displaystyle=$
$\displaystyle 1-\frac{n}{2}a+\sum_{k=2}^{n-1}\binom{n}{k}B_{k}a^{k}.$

Next, by using Corollary 1.4, we have for any $k\in\\{2,3,\dots,n-1\\}$ and
any prime divisor $p$ of $a$ (so $\vartheta_{p}(a)\geq 1$):

$\displaystyle\vartheta_{p}\left(B_{k}a^{k}\right)$ $\displaystyle=$
$\displaystyle\vartheta_{p}\left(B_{k}\right)+(k-1)\vartheta_{p}(a)+\vartheta_{p}(a)$
$\displaystyle\geq$ $\displaystyle-1+(k-1)+\vartheta_{p}(a)$ $\displaystyle=$
$\displaystyle k-2+\vartheta_{p}(a)$ $\displaystyle\geq$
$\displaystyle\vartheta_{p}(a),$

showing that $B_{k}a^{k}\equiv 0\ (\mathrm{mod}\ a)$ ($\forall
k\in\\{2,3,\dots,n-1\\}$). By substituting these congruences into (2.1), we
conclude to the required congruence $G_{n,a}\equiv 1-\frac{n}{2}a\
(\mathrm{mod}\ a)$. ∎

From Theorem 2.3, we immediately derive the following corollary:

######  Corollary 2.4.

Let $a\geq 2$ be an integer.

  * •

If $a$ is odd then we have for any positive integer $n$:

$G_{n,a}~{}\equiv~{}1\ (\mathrm{mod}\ a).$

  * •

If $a$ is even then we have for any integer $n\geq 2$:

$G_{n,a}~{}\equiv~{}\begin{cases}1&\text{if $n$ is even}\\\
1+\frac{a}{2}&\text{if $n$ is odd}\end{cases}\ (\mathrm{mod}\ a).$ $\square$

######  Remark 2.5.

By applying Corollary 2.4 for $a=2$, we obtain the well-known result stating
that the usual Genocchi numbers of even positive orders are all odd.

We conclude this note by the following corollary which is an immediate
consequence of Corollary 2.4 above.

######  Corollary 2.6.

For any integers $a$ and $n$, both greater than $1$, we have:

$\gcd\left(G_{n,a},a\right)~{}\in~{}\\{1,2\\}.$

Besides, the equality $\gcd(G_{n,a},a)=2$ holds if and only if $a\equiv 2\
(\mathrm{mod}\ 4)$ and $n$ is odd. $\square$

## References

  * [1] L. Comtet. Advanced Combinatorics. The Art of Finite and Infinite Expansions, revised and enlarged ed., D. Reidel Publ. Co., Dordrecht, 1974.
  * [2] D. Dumont & D. Foata. Une propriété de symétrie des nombres de Genocchi, Bull. Soc. Math. France, 104 (1976), p. 433-451.
  * [3] B. Farhi. A new generalization of the Genocchi numbers and its consequence on the Bernoulli polynomials, preprint, 2020. Available at https://arxiv.org/abs/2012.01969.
  * [4] J. Gandhi. A conjectured representation of Genocchi numbers, Amer. Math. Monthly, 77 (1970), p. 505-506.
  * [5] A. Genocchi. Intorno all espressioni generali di numeri Bernoulliani, Annali di scienze mat. e fisiche, compilati da Barnaba Tortolini, 3 (1852), p. 395-405.
  * [6] G.H. Hardy & E.M. Wright. The theory of numbers, 5th ed., Oxford Univ. Press, London, 1979.
  * [7] N. Nielsen. Traité élémentaire des nombres de Bernoulli, Gauthier–Villars, Paris, 1923.
  * [8] J. Riordan & P. Stein. Proof of a conjecture on Genocchi numbers, Disc. Math., 5 (1973), p. 381-388.


# About optimal loss function for training
physics-informed neural networks
under respecting causality

Vasiliy A. Es’kin
Department of Radiophysics, University of Nizhny Novgorod
Nizhny Novgorod, Russia, 603950
and
Manpower IT Solutions, Nizhny Novgorod, 603140, Russia
<EMAIL_ADDRESS>
& Danil V. Davydov
Mechanical Engineering Research Institute
Russian Academy of Sciences
Nizhny Novgorod, Russia, 603155
and
Manpower IT Solutions
Nizhny Novgorod, 603140, Russia
<EMAIL_ADDRESS>
& Ekaterina D. Egorova
Huawei Nizhny Novgorod Research Center
Nizhny Novgorod, Russia
and
Institute of Applied Physics
Russian Academy of Sciences
Nizhny Novgorod, Russia, 603155
<EMAIL_ADDRESS>
& Alexey O. Malkhanov
Huawei Nizhny Novgorod Research Center
Nizhny Novgorod, Russia
<EMAIL_ADDRESS>
& Mikhail A. Akhukov
Manpower IT Solutions
Nizhny Novgorod, 603140, Russia
<EMAIL_ADDRESS>
& Mikhail E. Smorkalov
Skolkovo Institute of Science and Technology
Moscow, Russia
and
Huawei Nizhny Novgorod Research Center
Nizhny Novgorod, Russia
<EMAIL_ADDRESS>
Corresponding author: Vasiliy Alekseevich Es’kin<EMAIL_ADDRESS>
###### Abstract

A method is presented that allows to reduce a problem described by
differential equations with initial and boundary conditions to the problem
described only by differential equations. The advantage of using the modified
problem for physics-informed neural networks (PINNs) methodology is that it
becomes possible to represent the loss function in the form of a single term
associated with differential equations, thus eliminating the need to tune the
scaling coefficients for the terms related to boundary and initial conditions.
The weighted loss functions respecting causality were modified and new
weighted loss functions based on generalized functions are derived. Numerical
experiments have been carried out for a number of problems, demonstrating the
accuracy of the proposed methods.

_Keywords_ Deep Learning $\cdot$ Physics-informed Neural Networks $\cdot$
Partial differential equations $\cdot$ Predictive modeling $\cdot$
Computational physics $\cdot$ Nonlinear dynamics

##  1 Introduction

The last decade has been notable by significant achievements in the field of
machine learning such as deep learning [1]. The development of training
methods and the improvement of deep neural network architectures have led to
the creation of applications that are comparable or superior to the average
person in many fields of activity previously available only to humans [2, 3,
4]. Among these applications, machine translation systems [5, 6, 7], image
generators by description [8, 9], game systems [10, 11] and text generators
[12] can be distinguished. The current revolution in the field of deep
learning has not left aside the areas of human activities related to the
obtaining and use of scientific knowledge. Auxiliary systems for solving
scientific problems can be divided into those that are trained using
previously obtained data (from experiments or computations) and build their
predictions on this basis [13, 14, 15, 16, 17], those that leverage previously
established physical laws [18, 19, 20] and those that utilize the both
approaches [21, 22]. The first category includes the widely known Alphafold
[23] and, somewhat less famous, analysis system of interaction of laser
radiation with the medium [24, 25]. Data-driven methods are out of
consideration of this paper, and we rather focus on the second category
exclusively. Initial ideas for constraining neural networks using physical
laws have been explored in [26, 27] studies of previous century. Though more
contemporary research in the field was done few years ago in physics-informed
neural networks (PINNs) [18] with help of modern computational tools. In the
PINN approach, a neural network is trained to approximate the dependences of
physical values on spatial and temporal coordinates for a given problem
described by physical equations together with initial and boundary conditions,
and sometimes with data from computations or experiments. This method has been
used in recent years to solve a wide range of problems described by ordinary
differential equations [21, 22, 18, 28], integro-differential equations [29],
nonlinear partial differential equations [21, 22, 18, 30, 31, 32], PDE with
noisy data [33], etc. [34, 35], related to various fields of knowledge such as
thermodynamics [36, 32], hydrodynamics [37, 38, 30, 31], mechanics [39],
electrodynamics [40, 32], geophysics [41], finance [29], etc.

Despite of the PINN success, it is often necessary to adapt hyperparameters of
learning and the architecture of the neural network for each new problem. In
particular, training a neural network boils down to minimization of the loss
function, which includes terms related to differential equations, initial and
boundary conditions [21, 18, 19]. In practice, for the training procedure to
converge, it is necessary to tune the weighting coefficients for each of the
terms mentioned. This tuning is done either empirically or by finding the
optimum [42] associated with significant computational cost. In a number of
tasks, the weighing is done not only for these different terms, but also for
various collocation points chosen in the domain of interest [43]. Such
weighting is based either on the minimization of the general loss function, or
based on the features of some behavior, in particular, on the causality
principle [19], which has been efficiently used for a number of tasks and is
actively leveraged to modify other methods and accelerate their convergence
[44]. Notwithstanding the successes achieved, the PINN approach still
continues to evolve and its development requires clear and transparent methods
for weightings of coefficients of the elements of the general loss function
and the architecture of the neural network for the problem.

This work is devoted to the development and application a number of techniques
that can improve the accuracy of neural networks trained on the basis of the
PINN approach. We will briefly list the contributions made in the paper:

  1. 1.

Reformulation of the original problem described by a differential equation
accompanied with initial and boundary conditions to the form in which takes
advantage of the differential equation only.

  2. 2.

Based on the reformulated problem, the loss function with the single
significant term is proposed.

  3. 3.

Further, the modified loss function based on the causality principle is
derived to account the possibility of using it for large number of temporal
slices.

  4. 4.

Loss functions based on causality and spatial-temporal locality principles are
proposed within the framework of the formalism of generalized functions
(Heaviside step function and Dirac delta function).

  5. 5.

A number of experiments have been carried out for the above techniques for
fully connected neural networks with homogeneous and heterogeneous activation
functions.

The paper is structured as follows. In section 2, the original PINN method and
suggestions for modifications of the problems and loss function are presented.
Section 3 includes the original description of causal training and
modifications of the causal training and its improvement to spatial-temporal
training within the framework of the formalism of generalized functions. In
Section 4, numerical experiments to illustrate the accuracy and the efficiency
of presented methods are given. Finally, in the Section 5 concluding remarks
are given.

##  2 Formulation of the problem

###  2.1 Physical-informed neural networks (PINNs)

Consider nonlinear partial differential equations (PDEs), which depend on time
$t$ and spatial coordinates ${\bf x}$, and in general take the following form

${\bf u}_{t}+\mathcal{N}[{\bf u}]=0,\quad t\in[0,T],\quad{\bf x}\in\Omega,$
(1)

under the initial and boundary conditions

$\displaystyle{\bf u}(0,{\bf x})={\bf g}({\bf x}),\quad{\bf x}\in\Omega,$ (2)
$\displaystyle\mathcal{B}[{\bf u}]=0,\quad t\in[0,T],\quad{\bf
x}\in\partial\Omega,$ (3)

where ${\bf u}(t,{\bf x})$ denotes the latent solution that is governed by the
PDE system of equation (1), ${\bf u}_{t}$ is temporal derivative,
$\mathcal{N}$ is a nonlinear differential operator, $\mathcal{B}$ is a
boundary operator corresponding to boundary conditions (Dirichlet, Neumann,
Robin, periodic conditions), ${\bf{g}}({\bf x})$ is initial distribution of
${\bf u}$, $\Omega$ and $\partial\Omega$ are spatial domain and its boundary,
respectively.

According to PINN approach [18], which stands on the shoulders of the
universal approximation theorem [45], the aim of PINN is to approximate the
unknown solution ${\bf u}(t,{\bf x})$ with a deep neural network ${\bf
u}_{{\bm{\theta}}}(t,{\bf x})$, where ${\bm{\theta}}$ denote all trainable
parameters (e.g., weight and biases) of the network. Finding optimal
parameters is an optimization problem, which requires definition of a loss
function such that its minimum gives the solution of the PDE. The physical-
informed model is trained by minimizing the composite loss function which
consists of the local residuals of the differential equation over the problem
domain and its boundary as shown below:

$\mathcal{L}({\bm{\theta}})=\lambda_{ic}\mathcal{L}_{ic}({\bm{\theta}})+\lambda_{bc}\mathcal{L}_{bc}({\bm{\theta}})+\lambda_{r}\mathcal{L}_{r}({\bm{\theta}}),$
(4)

where

$\displaystyle\mathcal{L}_{ic}\left({\bm{\theta}}\right)=\frac{1}{N_{ic}}\sum_{i=1}^{N_{ic}}\left|{\bf
u}_{{\bm{\theta}}}\left(0,{\bf x}_{i}^{(ic)}\right)-{\bf g}\left({\bf
x}_{i}^{(ic)}\right)\right|^{2},$ (5)
$\displaystyle\mathcal{L}_{bc}\left(\bm{\theta}\right)=\frac{1}{N_{bc}}\sum_{i=1}^{N_{bc}}\left|\mathcal{B}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t_{i}^{(bc)},{\bf
x}_{i}^{(bc)}\right)\right|^{2},$ (6)
$\displaystyle\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\left|\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t_{i}^{(r)},{\bf
x}_{i}^{(r)}\right)\right|^{2},$ (7) $\displaystyle\mathcal{R}\left[{\bf
u}\right]:={\bf u}_{t}+\mathcal{N}\left[{\bf u}\left(t,{\bf x}\right)\right].$
(8)

Here $\left\\{{\bf x}_{i}^{(ic)}\right\\}^{N_{ic}}_{i=1}$,
$\left\\{{t}_{i}^{(bc)},{\bf x}_{i}^{(bc)}\right\\}^{N_{bc}}_{i=1}$, and
$\left\\{{t}_{i}^{(r)},{\bf x}_{i}^{(r)}\right\\}^{N_{r}}_{i=1}$ are sets of
points corresponding to initial condition domain, boundary condition domain
and PDE domain, respectively. These points can be the vertices of a fixed mesh
or can be randomly sampled at each iteration of a gradient descent algorithm.
All required gradients w.r.t. input variables ($t$ and $\bf x$) and network
parameters $\theta$ can be efficiently computed via automatic differentiation
[46] with algorithmic accuracy, which is defined by the accuracy of
computation system. Note, the hyperparameters $\lambda_{ic}$, $\lambda_{bc}$
and $\lambda_{r}$ allow for separate tuning of learning rate for each of the
loss terms in order to improve convergence of the model [47, 48].

The optimization problem can be defined as follow

${\bm{\theta}}^{*}=\underset{{\bm{\theta}}}{\arg}\min\mathcal{L}({\bm{\theta}}),$
(9)

where ${\bm{\theta}}^{*}$ are optimal parameters of the neural network which
minimize the discrepancy between the exact unknown solution ${\bf u}$ and the
approximate one ${\bf u}_{\boldsymbol{\theta}^{*}}$. Schematic diagram of the
general PINN method, which is used throughout this paper, is shown on the
Figure 1.

Let’s rewrite the loss components (5)–(7) in the continuous form for the
further analysis. These representations are

$\displaystyle\mathcal{L}_{ic}\left({\bm{\theta}}\right)=\frac{1}{V_{\Omega}}\int_{\Omega}\left|{\bf
u}_{{\bm{\theta}}}\left(0,{\bf x}\right)-{\bf g}\left({\bf
x}\right)\right|^{2}d{\bf x},$ (10)
$\displaystyle\mathcal{L}_{bc}\left(\bm{\theta}\right)=\frac{1}{T}\int_{0}^{T}\left|\mathcal{B}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t,{\bf x}\right)\right|^{2}dt,\quad{\bf
x}\in{\partial\Omega},$ (11)
$\displaystyle\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{T}\frac{1}{V_{\Omega}}\int^{T}_{0}dt\int_{\Omega}d{\bf
x}\left|\mathcal{R}\left[{\bf u}_{{\bm{\theta}}}\right]\left(t,{\bf
x}\right)\right|^{2},$ (12)

where $V_{\Omega}$ is the “volume” of domain $\Omega$.

Figure 1: Schematic diagram of the general PINN method.

###  2.2 Reformulation of the problem

####  2.2.1 Problem in time domain

Let’s modify the problem (1)–(3) to the one which is free from the initial
conditions. For this we expand the time interval from infinitely distant
moment in the past to the time $T$. On time interval $(-\infty,0]$ the value
of given process ($\bf U$) is assumed to be constant in time and has a spatial
distribution described by the initial distribution (2). The value $\bf U$ has
the following form:

${\bf U}={\bf g}({\bf x})\chi(-t)+{\bf u}(t,{\bf x})\chi(t),\quad
t\in(-\infty,T],$ (13)

where ${\bf u}(t,{\bf x})$ is unknown function that supposed to be a solution
of the system of equations (1)–(3), $\chi(t)$ is Heaviside step function [49]
defined as

$\chi(t)=\left\\{\begin{array}[]{ll}0,&t<0,\\\ 1/2,&t=0,\\\
1,&t>0.\end{array}\right.$ (14)

In order the new function $\bf U$ to be the solution at time interval
$(-\infty,0]$ the initial equation (1) should be supplied with additional
term. Such term corresponds to “force” which disappears instantly at the
moment $t=0$.

Without the loss of generality and for the sake of simplicity, we will
consider the problem in an unbounded domain, i.e. without boundary conditions.
As result the problem (1)–(3) is reduced to the following differential
equations

${\bf U}_{t}+\mathcal{N}[{\bf U}]=\mathcal{N}[{\bf g}({\bf x})]\chi(-t),\quad
t\in(-\infty,T],\quad{\bf x}\in\Omega,$ (15)

where unknown function ${\bf u}(t,{\bf x})$ is encapsulated into $\bf U$. As
previously, we approximate the unknown solution ${\bf U}(t,{\bf x})$ with a
deep neural network ${\bf U}_{{\bm{\theta}}}(t,{\bf x})={\bf g}({\bf
x})\chi(-t)+{\bf u}_{{\bm{\theta}}}(t,{\bf x})\chi(t)$.

In this problem statement, the loss function for problem (15) contains only
the term corresponding to residuals of differential equation:

$\displaystyle\mathcal{L}({\bm{\theta}})=\mathcal{L}_{r}({\bm{\theta}}),$ (16)

where

$\displaystyle\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{T}\frac{1}{V_{\Omega}}\int^{T}_{-\infty}dt\int_{\Omega}d{\bf
x}\left|\mathcal{R}^{*}\left[{\bf U}_{{\bm{\theta}}},{\bf
g}\right]\left(t,{\bf x}\right)\right|^{2},$ (17)
$\displaystyle\mathcal{R}^{*}\left[{\bf u},{\bf g}\right]:={\bf
u}_{t}+\mathcal{N}\left[{\bf u}\left(t,{\bf x}\right)\right]-\mathcal{N}[{\bf
g}({\bf x})]\chi(-t).$ (18)

Expanding the integrand in (17) we have

$\mathcal{R}^{*}\left[{\bf U}_{{\bm{\theta}}},{\bf g}\right]\left(t,{\bf
x}\right){=}\left\\{\begin{array}[]{ll}\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t,{\bf x}\right){+}\left[{\bf
u}_{{\bm{\theta}}}\left(t,{\bf x}\right){-}{\bf g}\left({\bf
x}\right)\right]\delta(t),&t\geq 0,\\\ 0&t<0.\end{array}\right.$ (19)

Here $\delta(t)$ is Dirac delta function [50]. We use the following properties
of delta function and Heaviside step function [50]

$\displaystyle\delta(x)=\frac{\partial\chi(x)}{\partial
x},\quad\chi(-x)=1-\chi(x).$ (20)

Substitution of (19) to (17) gives expression for loss in the following
representation

$\displaystyle\mathcal{L}_{r}\left(\bm{\theta}\right)$
$\displaystyle=\frac{1}{T}\frac{1}{V_{\Omega}}\lim_{t\to
0}\delta(t)\int_{\Omega}d{\bf x}\left|{\bf u}_{{\bm{\theta}}}\left(0,{\bf
x}\right){-}{\bf g}\left({\bf x}\right)\right|^{2}$
$\displaystyle+2\frac{1}{T}\frac{1}{V_{\Omega}}\int_{\Omega}d{\bf
x}\left|\mathcal{R}\left[{\bf u}_{{\bm{\theta}}}\right]\left(0,{\bf
x}\right)\left[{\bf u}_{{\bm{\theta}}}\left(0,{\bf x}\right){-}{\bf
g}\left({\bf x}\right)\right]\right|$
$\displaystyle+\frac{1}{T}\frac{1}{V_{\Omega}}\int^{T}_{0}dt\int_{\Omega}d{\bf
x}\left|\mathcal{R}\left[{\bf u}_{{\bm{\theta}}}\right]\left(t,{\bf
x}\right)\right|^{2}.$ (21)

Here we have taken into account the property of delta function which states
its dominance in the finite sum, i.e.
$\left|a\delta(x)+b\right|^{2}\simeq\left|a\right|^{2}\delta^{2}(x)+2\left|ab\right|\delta(x)+\left|b\right|^{2}$.

###  2.3 Optimization problem reformulation

Then we apply the well known representation of the Dirac delta function [50]

$\displaystyle\delta(x)=\lim_{\sigma\to 0}\frac{\sigma}{x^{2}+\sigma^{2}}.$
(22)

We obtain the following formulation of the loss function for (9)

$\displaystyle\mathcal{L}\left(\bm{\theta}\right)$
$\displaystyle=\frac{1}{T}\frac{1}{V_{\Omega}}\int^{T}_{0}dt\int_{\Omega}d{\bf
x}\left|\mathcal{R}\left[{\bf u}_{{\bm{\theta}}}\right]\left(t,{\bf
x}\right)\right|^{2}$
$\displaystyle+\frac{1}{T}\frac{1}{V_{\Omega}}\left\\{2\int_{\Omega}d{\bf
x}\left|\mathcal{R}\left[{\bf u}_{{\bm{\theta}}}\right]\left(0,{\bf
x}\right)\left[{\bf u}_{{\bm{\theta}}}\left(0,{\bf x}\right){-}{\bf
g}\left({\bf x}\right)\right]\right|\right.$
$\displaystyle\left.+\lim_{\begin{subarray}{c}\sigma\to 0\\\ t\to
0\end{subarray}}\frac{\sigma}{t^{2}+\sigma^{2}}\int_{\Omega}d{\bf x}\left|{\bf
u}_{{\bm{\theta}}}\left(0,{\bf x}\right){-}{\bf g}\left({\bf
x}\right)\right|^{2}\right\\}.$ (23)

To proceed from continuous time to the discrete time representation, we divide
the time interval $[0,T]$ into $N_{t}$ intervals and use the rectangular rule
of integration for this mesh with step size $T/N_{t}$. We take
$\sigma=T/N_{t}$ for the approximation of the delta function in the time
domain. This representation of delta function is valid as it preserves its
fundamental property [50] [$\int_{-\infty}^{\infty}f(x)\delta(x-a)dx=f(a)$]
for the numerical integration. Finally, we get the following loss function

$\displaystyle\mathcal{L}\left(\bm{\theta}\right)$
$\displaystyle=\frac{1}{N_{t}}\frac{1}{V_{\Omega}}\sum\limits_{i=0}^{N_{t}}\int_{\Omega}d{\bf
x}\left|\mathcal{R}\left[{\bf u}_{{\bm{\theta}}}\right]\left(t_{i},{\bf
x}\right)\right|^{2}$
$\displaystyle+\frac{1}{T}\frac{1}{V_{\Omega}}\left\\{2\int_{\Omega}d{\bf
x}\left|\mathcal{R}\left[{\bf u}_{{\bm{\theta}}}\right]\left(t_{0},{\bf
x}\right)\left[{\bf u}_{{\bm{\theta}}}\left(t_{0},{\bf x}\right){-}{\bf
g}\left({\bf x}\right)\right]\right|\right.$
$\displaystyle\left.+\frac{N_{t}}{T}\int_{\Omega}d{\bf x}\left|{\bf
u}_{{\bm{\theta}}}\left(t_{0},{\bf x}\right){-}{\bf g}\left({\bf
x}\right)\right|^{2}\right\\}.$ (24)

Here we take $t_{0}=0$. Note, the multiplier ${N_{t}}$ of third term of (24)
coincides with coefficient $\lambda_{ic}$ empirically found in [19, 20, 51] in
the case of $N_{t}=100$.

With this representation, the approximation of delta function is associated
with a numerical integration scheme. For a number of tasks though, an
approximation of the delta function is possible that is not tied to
integration rules. To achieve that the value $\sigma$ in (22) can be taken as
an estimated relaxation time of the initial conditions to zero, or of doubling
or halving values of initial conditions. Note that obtained approximation of
the delta function can be used as a weight $\lambda_{ic}$ in the original
training scheme described by loss terms (5)–(7). Example of deriving such
approximation for specific problem is provided in section 4.3.1.

####  2.3.1 MAE loss

Interestingly, if instead of MSE loss we use MAE loss for optimization problem
(9), then no additional representation of the Dirac delta function is required
and instead of loss (23) we have the following function

$\displaystyle\mathcal{L}\left(\bm{\theta}\right)$
$\displaystyle=\frac{1}{T}\frac{1}{V_{\Omega}}\left(\int^{T}_{0}dt\int_{\Omega}d{\bf
x}\left|\mathcal{R}\left[{\bf u}_{{\bm{\theta}}}\right]\left(t,{\bf
x}\right)\right|+\int_{\Omega}d{\bf x}\left|{\bf
u}_{{\bm{\theta}}}\left(0,{\bf x}\right){-}{\bf g}\left({\bf
x}\right)\right|\right).$ (25)

Based on our experiments, training a neural network with such a loss function
converges extremely weak to the required solution, despite the simplicity of
the equation (25). For this reason, in all of the experiments presented in
section 4, MSE loss has been used.

###  2.4 Problem in spatial-temporal domain

We can reformulate the problem in spatial domain in the same manner. For
greater clarity and brevity, and without violating generality, we will
consider a problem with dependence on only one spatial coordinate $x$.

Let spatial domain be determined by a coordinate $x$ ($x\in[X_{1},X_{2}]$),
and problem (1)–(3) is subject to the following boundary conditions

$\displaystyle{\bf u}(t,X_{1})={\bf f_{1}}(t),\quad{\bf u}(t,X_{2})={\bf
f_{2}}(t).$ (26)

Then, artificially extended solution of (1)–(3) takes the following form

$\displaystyle{\bf U}$ $\displaystyle={\bf g}({x})\chi(-t)+{\bf
u}(t,{x})\chi(t)\left[\chi(x-X_{1})-\chi(x-X_{2})\right]$
$\displaystyle+\left[{\bf f_{1}}(t)\chi(X_{1}-x)+{\bf
f_{2}}(t)\chi(x-X_{2})\right]\chi(t),$ (27) $\displaystyle{\text{for the
domain }}t\in(-\infty,T],\quad x\in(-\infty,\infty).$

In this case, the problem (1)–(3) is reduced to the following differential
equations

$\displaystyle{\bf U}_{t}+\mathcal{N}[{\bf U}]$
$\displaystyle=\mathcal{N}[{\bf g}({x})]\chi(-t)$
$\displaystyle+\left[\frac{\partial{\bf f_{1}}(t)}{\partial
t}\chi(X_{1}-x)+\frac{\partial{\bf f_{2}}(t)}{\partial
t}\chi(x-X_{2})\right]\chi(t)$ $\displaystyle+\left\\{\mathcal{N}\left[{\bf
f_{1}}(t)\right]\chi(x-X_{1})+\mathcal{N}\left[{\bf
f_{2}}(t)\right]\chi(X_{2}-x)\right\\}\chi(t),\quad t\in(-\infty,T],\quad
x\in(-\infty,\infty).$ (28)

The value $\mathcal{R}^{*}\left[{\bf U}_{{\bm{\theta}}},{\bf
g}\right]\left(t,{x}\right)$ has form similar to (19)

$\mathcal{R}^{*}\left[{\bf U}_{{\bm{\theta}}},{\bf
g}\right]\left(t,{x}\right){=}\left\\{\begin{array}[]{lll}\mathcal{R}\left[{\bf
U}^{*}_{{\bm{\theta}}}\right]\left(t,{x}\right)&{+}\left[{\bf
u}_{{\bm{\theta}}}\left(t,{x}\right){-}{\bf
g}\left({x}\right)\right]\delta(t)&\\\ &-\left[\frac{\partial{\bf
f_{1}}(t)}{\partial t}\chi(X_{1}-x)+\frac{\partial{\bf f_{2}}(t)}{\partial
t}\chi(x-X_{2})\right]&\\\ &-\left\\{\mathcal{N}\left[{\bf
f_{1}}(t)\right]\chi(x-X_{1})+\mathcal{N}\left[{\bf
f_{2}}(t)\right]\chi(X_{2}-x)\right\\},&t\geq 0,\quad x\in[X_{1},X_{2}],\\\
0,&&t<0,\quad\text{or }x\notin[X_{1},X_{2}],\end{array}\right.$ (29)

where

$\displaystyle{\bf U}^{*}_{{\bm{\theta}}}$ $\displaystyle={\bf
u}_{{\bm{\theta}}}(t,{x})\left[\chi(x-X_{1})-\chi(x-X_{2})\right]+{\bf
f_{1}}(t)\chi(X_{1}-x)+{\bf f_{2}}(t)\chi(x-X_{2}).$ (30)

The representation of $\mathcal{R}\left[{\bf U}^{*}_{{\bm{\theta}}}\right]$ in
terms of $\mathcal{R}\left[{\bf u}_{{\bm{\theta}}}\right]$ similarly to
equation (19) depends on the specific type of differential equation. Example
of such representation can be found in section 4.2.1.

##  3 Causal training for physics-informed neural network

It was suggested in [19] to reformulate the loss function so that it would
respect the physical causality when solving PDEs with the help of PINNs:

$\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{N_{t}}\sum_{i=0}^{N_{t}}{w}_{i}\mathcal{L}_{r}\left(t_{i},{\bm{\theta}}\right),$
(31)

where weights $w_{i}$ are

$w_{i}=\exp\left[-\varepsilon\sum_{k=0}^{i-1}\mathcal{L}_{r}\left(t_{k},{\bm{\theta}}\right)\right],\quad\text{for
}i=1,2,3,\ldots,N_{t},\text{and }w_{0}=1,$ (32)

and where $\varepsilon$ is causality parameter that controls the steepness of
the weights $w_{i}$. The temporal residual loss
[$\mathcal{L}_{r}\left(t_{i},{\bm{\theta}}\right)$] for the total loss (23) is

$\displaystyle\mathcal{L}_{r}\left(t_{i},{\bm{\theta}}\right)=$
$\displaystyle\frac{1}{V_{\Omega}}\int_{\Omega}d{\bf
x}\left\\{\left|\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t_{i},{\bf x}\right)\right|^{2}\right.$
$\displaystyle+2\frac{N_{t}}{T}\delta_{0i}\left|\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t_{0},{\bf x}\right)\left[{\bf
u}_{{\bm{\theta}}}\left(t_{0},{\bf x}\right){-}{\bf g}\left({\bf
x}\right)\right]\right|$
$\displaystyle\left.+\frac{N_{t}^{2}}{T^{2}}\delta_{0i}\left|{\bf
u}_{{\bm{\theta}}}\left(t_{0},{\bf x}\right){-}{\bf g}\left({\bf
x}\right)\right|^{2}\right\\},$ (33)

$\delta_{0i}$ is the Kronecker symbol.

In the case when domain $\Omega$ consists of single spatial dimension $x$
($x\in[X_{1},X_{2}]$, $V_{\Omega}=X_{2}-X_{1}$) the temporal residual loss has
discrete form under regular mesh with $N_{x}$th nodes

$\displaystyle\mathcal{L}_{r}\left(t_{i},{\bm{\theta}}\right)=$
$\displaystyle\frac{1}{N_{x}}\sum_{n=0}^{N_{x}}\left\\{\left|\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t_{i},{x_{n}}\right)\right|^{2}\right.$
$\displaystyle+2\frac{N_{t}}{T}\delta_{0i}\left|\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t_{0},{x_{n}}\right)\left[{\bf
u}_{{\bm{\theta}}}\left(t_{0},{x_{n}}\right){-}{\bf
g}\left({x_{n}}\right)\right]\right|$
$\displaystyle\left.+\frac{N_{t}^{2}}{T^{2}}\delta_{0i}\left|{\bf
u}_{{\bm{\theta}}}\left(t_{0},{x_{n}}\right){-}{\bf
g}\left({x_{n}}\right)\right|^{2}\right\\}.$ (34)

The resulting weighted residual loss is

$\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{N_{t}}\sum_{i=0}^{N_{t}}\exp\left[-\varepsilon\sum_{k=0}^{i-1}\mathcal{L}_{r}\left(t_{k},{\bm{\theta}}\right)\right]\mathcal{L}_{r}\left(t_{i},{\bm{\theta}}\right).$
(35)

Authors of [19] recognized that the weights $w_{i}$ should be large — and
therefore allow the minimization of
$\mathcal{L}_{r}\left(t_{i},{\bm{\theta}}\right)$ — only if all residuals
$\left\\{\mathcal{L}_{r}\left(t_{k},{\bm{\theta}}\right)\right\\}_{k=0}^{i}$
before $t_{i}$ are minimized properly, and vice versa. As a consequence,
$\mathcal{L}_{r}\left(t_{i},{\bm{\theta}}\right)$ is not minimized unless all
previous residuals
$\left\\{\mathcal{L}_{r}\left(t_{k},{\bm{\theta}}\right)\right\\}_{k=0}^{i-1}$
decrease to some small value such that $w_{i}$ is large enough.

###  3.1 Modification of the causal training for physics-informed neural
network

Despite the causal training method significantly improves the convergence of
neural network learning this approach does not scale for time grids with the
growth of number of nodes. This can be explained by the fact that the exponent
of the weighting function (32) contains the cumulative sum. This sum is
increasing with the increase of the $N_{t}$ and such behavior hinders learning
of neural network for large $t$. Note that this feature in row of cases can be
lead to the best convergence of training in compared to other methods of
weighting.

This drawback becomes obvious when trying to write the residual loss (35) in
continuous form

$\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{T}\int_{0}^{T}\lim_{N_{t}\to\infty}\left[\exp\left(-\varepsilon\frac{N_{t}}{t}\int_{0}^{t}\mathcal{L}_{r}\left(t^{\prime},{\bm{\theta}}\right)dt^{\prime}\right)\right]\mathcal{L}_{r}\left(t,{\bm{\theta}}\right)dt.$
(36)

This integral tends to zero for any $\varepsilon$ and time due to unlimited
increase of $N_{t}$. Such disadvantage prevents mathematical analysis of the
loss function using standard methods based on integral and differential
calculus.

In order to mitigate the above drawback, the residual loss (35) can be
rewritten in the following form

$\displaystyle\mathcal{L}_{r}\left(\bm{\theta}\right)$
$\displaystyle=\frac{1}{N_{t}}\sum_{i=0}^{N_{t}}w_{i}\mathcal{L}_{r}\left(t_{i},{\bm{\theta}}\right),$
(37) $\displaystyle w_{i}$ $\displaystyle=\exp\left[-\varepsilon
s_{i}\right],$ (38) $\displaystyle s_{i}$
$\displaystyle=\frac{1}{i}\sum_{k=0}^{i-1}\mathcal{L}_{r}\left(t_{k},{\bm{\theta}}\right).$
(39)

The continuous form of residual loss (37) is non zero and has the following
representation

$\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{T}\int_{0}^{T}dt\exp\left[-\frac{\varepsilon}{t}\int_{0}^{t}dt^{\prime}\mathcal{L}_{r}\left(t^{\prime},{\bm{\theta}}\right)\right]\mathcal{L}_{r}\left(t,{\bm{\theta}}\right).$
(40)

###  3.2 Alternatives of the causal training for physics-informed neural
network

###  3.3 Heaviside–weighted residual loss

Let’s consider residual loss weighted using Heaviside step function

$\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{T}\int_{0}^{T}\left\\{1-\chi\left[S(t)\right]\right\\}\mathcal{L}_{r}\left(t,{\bm{\theta}}\right)dt,$
(41)

where

$S(t)=\frac{1}{t}\int_{0}^{t}\mathcal{L}_{r}\left(t^{\prime},{\bm{\theta}}\right)dt^{\prime}.$
(42)

In this case term $1-\chi(x)$ in the integrand of equation (41) does not tend
to zero under condition $x\leq 0$. It means that causality condition is
satisfied. To calculate the loss we use the following well known
representations of the Heaviside step function [50]

$\chi(x)=\lim_{\varepsilon\to\infty}\frac{1}{1+\exp\left(-2\varepsilon
x\right)}.$ (43)

The equation (43) leads to the following representation of residual loss

$\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{T}\int_{0}^{T}dt\frac{1}{1+\exp\left[2\varepsilon
S(t)\right]}\mathcal{L}_{r}\left(t,{\bm{\theta}}\right),$ (44)

which has the discrete form (31) with weights

$\displaystyle
w_{i}=\left[1+\exp\left({2\varepsilon}s_{i}\right)\right]^{-1}.$ (45)

####  3.3.1 Dirac–weighted residual loss

Let’s consider the extreme form of the weighted residual loss

$\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{T}\int_{0}^{T}\delta\left[S(t)\right]\mathcal{L}_{r}\left(t,{\bm{\theta}}\right)dt.$
(46)

Such loss function has well-pronounced property that it prohibits training for
large time stamps without have been previously trained for small ones.

To proceed with the discrete form of this loss function we will use the
following well known representation of the Dirac delta function [50]

$\displaystyle\delta(x)=\lim_{\sigma^{2}\to
0}\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{x^{2}}{2\sigma^{2}}\right).$
(47)

The integral form of residual loss in this case is

$\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{T}\int_{0}^{T}dt\exp\left\\{-\varepsilon\left[\frac{1}{t}\int_{0}^{t}dt^{\prime}t^{\prime}\mathcal{L}_{r}\left(t^{\prime},{\bm{\theta}}\right)\right]^{2}\right\\}\mathcal{L}_{r}\left(t,{\bm{\theta}}\right).$
(48)

Note, here the causality properties of weights are further enhanced with
multiplying by $t$ inside the inner integral. The discrete form of (48) in
case of equidistant time slices is

$\displaystyle\mathcal{L}_{r}\left(\bm{\theta}\right)=\frac{1}{N_{t}}\sum_{i=0}^{N_{t}}w_{i}\mathcal{L}_{r}\left(t_{i},{\bm{\theta}}\right),$
(49) $\displaystyle w_{i}=\exp\left(-\varepsilon\tilde{s}_{i}^{2}\right),$
(50)
$\displaystyle\tilde{s}_{i}=\frac{1}{i}\sum_{k=0}^{i-1}k\mathcal{L}_{r}\left(t_{k},{\bm{\theta}}\right).$
(51)

When deriving the equation (51) we have taken into account that
$t_{k}=kT/N_{t}$ and multiplier $T/N_{t}$ of $\tilde{s}_{i}$ was “hidden” in
causality parameter $\varepsilon$ in (50).

###  3.4 Expanding of the causal training approach to the spatial-temporal
domain

For problem (2.4) the loss function has the following representation

$\displaystyle\mathcal{L}\left(\bm{\theta}\right)$
$\displaystyle=\frac{1}{T}\frac{1}{V_{\Omega}}\int^{T}_{0}dt\int_{\Omega}d{x}\left|\mathcal{R}^{*}\left[{\bf
U}^{*}_{{\bm{\theta}}}\right]\left(t,{x}\right)\right|^{2}$
$\displaystyle+\frac{1}{T}\frac{1}{V_{\Omega}}\left\\{2\int_{\Omega}d{x}\left|\mathcal{R}^{*}\left[{\bf
U}^{*}_{{\bm{\theta}}}\right]\left(0,{x}\right)\left[{\bf
u}_{{\bm{\theta}}}\left(0,{x}\right){-}{\bf
g}\left({x}\right)\right]\right|\right.$
$\displaystyle\left.+\lim_{\begin{subarray}{c}\sigma\to 0\\\ t\to
0\end{subarray}}\frac{\sigma}{t^{2}+\sigma^{2}}\int_{\Omega}d{x}\left|{\bf
u}_{{\bm{\theta}}}\left(0,{x}\right){-}{\bf
g}\left({x}\right)\right|^{2}\right\\}.$ (52)

The weighted residual loss in the spatial-temporal domain can be defined as
follows

$\displaystyle\mathcal{L}^{*}_{r}\left(\bm{\theta}\right)=\frac{1}{N_{t}}\sum_{i=0}^{N_{t}}w_{i}\mathcal{L}^{*}_{r}\left(t_{i},{\bm{\theta}}\right),$
(53) $\displaystyle w_{i}=\exp\left[-\varepsilon s_{i}\right],\quad
s_{i}=\frac{1}{i}\sum_{k=0}^{i-1}\mathcal{L}^{*}_{r}\left(t_{k},{\bm{\theta}}\right),$
(54)

where

$\displaystyle\mathcal{L}^{*}_{r}\left(t_{i},{\bm{\theta}}\right)$
$\displaystyle=\frac{1}{N_{x}}\sum_{n=0}^{N_{x}}\left(w_{n}^{(\to)}+w_{n}^{(\leftarrow)}\right)\left\\{\left|\mathcal{R}\left[{\bf
U}^{*}_{{\bm{\theta}}}\right]\left(t_{i},{x_{n}}\right)\right|^{2}\right.$
$\displaystyle+2\frac{N_{t}}{T}\delta_{0i}\left|\mathcal{R}\left[{\bf
U}^{*}_{{\bm{\theta}}}\right]\left(t_{0},{x_{n}}\right)\left[{\bf
u}_{{\bm{\theta}}}\left(t_{0},{x_{n}}\right){-}{\bf
g}\left({x_{n}}\right)\right]\right|$
$\displaystyle\left.+\frac{N_{t}^{2}}{T^{2}}\delta_{0i}\left|{\bf
u}_{{\bm{\theta}}}\left(t_{0},{x_{n}}\right){-}{\bf
g}\left({x_{n}}\right)\right|^{2}\right\\},$ (55)

and

$\displaystyle w_{n}^{(\to)}$ $\displaystyle=\exp\left[-\varepsilon
s^{(\to)}_{n}\right],\quad w_{n}^{(\leftarrow)}=\exp\left[-\varepsilon
s^{(\leftarrow)}_{n}\right],$ $\displaystyle s^{(\to)}_{n}$
$\displaystyle=\frac{1}{n}\sum_{k=0}^{n-1}\tilde{\mathcal{L}}_{r}\left(x_{k},{\bm{\theta}}\right),\quad
s^{(\leftarrow)}_{n}=\frac{1}{N_{x}-n}\sum_{k=0}^{N_{x}-n-1}\tilde{\mathcal{L}}_{r}\left(x_{N_{x}-k},{\bm{\theta}}\right).$
(56)

The spatial residual loss is

$\displaystyle\tilde{\mathcal{L}}_{r}\left(x_{k},{\bm{\theta}}\right)$
$\displaystyle=\frac{1}{N_{t}}\sum_{i=0}^{N_{t}}\left|\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t_{i},{x_{k}}\right)\right|^{2}.$ (57)

Here superscripts $(\to)$ and $(\leftarrow)$ mean weighting in the forward and
backward directions along the $x$ axis, respectively (i.e. from $X_{1}$ to
$X_{2}$ and from $X_{2}$ to $X_{1}$).

###  3.5 Training procedure

Based on the above modifications of the problem, loss function and causality
weights, algorithms 1 and 2 present a causal and a spatial-temporal local
training algorithms for PINNs, respectively.

Data: Sets of coordinates $(t,{\bf x})$ and values ${\bf u}$ of points
corresponding to initial and boundary condition domains and set of coordinates
of collocation points

Result: Trained NN ${\bf u}_{\boldsymbol{\theta}}$ which approximates the
solution of PDE

Initialize temporal weights $w_{i}$ with 1, causal parameter $\varepsilon$ and
threshold parameter $\delta_{w}$ with chosen values.

Use $N$ epochs of a gradient descent algorithm to update the parameters
$\boldsymbol{\theta}$ as follows

for  _$n=1,\dots,N$_ do

1\. Calculate and update the temporal weight $w_{i}$ given by equations (38),
(45) or (50) using the temporal residual loss (34)

if  _$\underset{i}{\min}(\exp\left[-\varepsilon s_{i}\right]) >\delta_{w}$_
then

$\varepsilon\leftarrow m_{\varepsilon}\times\varepsilon$

2\. Calculate the loss using equations (16) and (31)

3\. Update the parameters $\theta$ via gradient descent

${\bm{\theta}}_{n}\leftarrow{\bm{\theta}}_{n-1}-\eta\nabla_{{\bm{\theta}}}\mathcal{L}\left(\bm{\theta}_{n-1}\right).$
(58)

From numerical experiments we found that the most promising hyperparameters
are $\delta_{w}=0.99$, initial value of $\varepsilon$=$10^{-3}$,
multiplication factor $m_{\varepsilon}=2$. Here $\eta$ is the learning rate of
gradient descent.

Algorithm 1 Training of physics-informed neural networks at time domain

Data: Sets of coordinates $(t,{\bf x})$ and values ${\bf u}$ of points
corresponding to initial and boundary condition domains and set of coordinates
of collocation points

Result: Trained NN ${\bf u}_{\boldsymbol{\theta}}$ which approximates the
solution of PDE

Initialize temporal and spatial weights $w_{i}$, $w_{i}^{(\to)}$, and
$w_{i}^{(\leftarrow)}$ with 1, causal parameter $\varepsilon$ and threshold
parameter $\delta_{w}$ with chosen values.

Use $N$ epochs of a gradient descent algorithm to update the parameters
$\boldsymbol{\theta}$ as follows

for  _$n=1,\dots,N$_ do

1\. Calculate and update the spatial weights $w_{i}^{(\to)}$ and
$w_{i}^{(\leftarrow)}$ given by equation (56) using the spatial residual loss
(57)

2\. Calculate and update the temporal weight $w_{i}$ given by equations (38),
(45) or (50) using the temporal residual loss (55)

if  _$\underset{i}{\min}(\exp\left[-\varepsilon s_{i}\right]) >\delta_{w}$ &
$\underset{i}{\min}\left(\exp\left[-\varepsilon
s_{i}^{(\to)}\right]\right)>\delta_{w}$ &
$\underset{i}{\min}\left(\exp\left[-\varepsilon
s_{i}^{(\leftarrow)}\right]\right)>\delta_{w}$_ then

$\varepsilon\leftarrow m_{\varepsilon}\times\varepsilon$

3\. Calculate the loss using equations (16) and (54)

4\. Update the parameters $\theta$ via gradient descent

${\bm{\theta}}_{n}\leftarrow{\bm{\theta}}_{n-1}-\eta\nabla_{{\bm{\theta}}}\mathcal{L}\left(\bm{\theta}_{n-1}\right).$
(59)

Algorithm 2 Training of the physics-informed neural networks at spatio-
temporal domain

Some details of the training procedures are explained below.

To save computational resources and to accelerate convergence for problems
with periodic boundary conditions along $x$ axis let’s construct a Fourier
feature embedding of input coordinates in the following form

$v(x)=\left[1,\cos(k_{x}x),\sin(k_{x}x),\cos(2k_{x}x),\sin(2k_{x}x),...,\cos(mk_{x}x),\sin(mk_{x}x)\right],$
(60)

where $k_{x}=2\pi/L$, $L$ is the minimal spatial period along $x$ direction,
and $m$ is some non-negative integer. With such embedding any trained
approximation of problem solution
${\bm{u}}_{\boldsymbol{\theta}}\left[v(x)\right]$ exactly satisfies periodic
conditions (see details in [19, 52]).

The networks had been trained via stochastic gradient descent using the Adam
optimizer [53] and different schedulers (see details in the table 5).

##  4 Numerical experiments

In this section several numerical experiments are presented to demonstrate the
accuracy and performance aforementioned methods.

To evaluate the accuracy of the approximate solution obtained with the help of
PINN method, the values of the solution of (1) predicted by the neural network
at given points are compared with the values calculated on the basis of
classical high-precision numerical methods. As a measure of accuracy, the
relative total ${\mathbb{L}}_{2}$ error of prediction is taken, which can be
expressed with the following relation

$\epsilon_{\rm error}=\left\\{\frac{1}{N_{e}}\sum^{N_{e}}_{i=1}\left[{\bf
u}_{\boldsymbol{\theta}}(t_{i},{\bf x}_{i})-{\bf u}(t_{i},{\bf
x}_{i})\right]^{2}\right\\}^{1/2}\times\left\\{\frac{1}{N_{e}}\sum^{N_{e}}_{i=1}\left[{\bf
u}(t_{i},{\bf x}_{i})\right]^{2}\right\\}^{-1/2},$ (61)

where $\left\\{{t}_{i},{\bf x}_{i}\right\\}^{N_{e}}_{i=1}$ is the set of
evaluation points taken from the domain $[0,T]\times\Omega$, ${\bf
u}_{\boldsymbol{\theta}}$ and ${\bf u}$ are the predicted and reference
solutions respectively.

For our experiments we used Pytorch [54] version 1.12.1 and the training was
carried out on a node with Nvidia Tesla V100 GPU.

###  4.1 Allen–Cahn equation

As in [18, 19] let’s consider a classical phase field modeled by one-
dimensional Allen–Cahn equation with periodic boundary condition

$\displaystyle u_{t}-0.0001u_{xx}+5u^{3}-5u=0,\quad t\in[0,1],\quad
x\in[-1,1],$ $\displaystyle u(x,0)=x^{2}\cos(\pi x),$ $\displaystyle
u(t,-1)=u(t,1),\quad u_{x}(t,-1)=u_{x}(t,1).$ (62)

We represent the latent variable $u$ by a fully-connected neural network
$u_{\boldsymbol{\theta}}$ with hyperbolic tangent activation function, 4
hidden layers with 128 neurons per each layer. For simplicity, a uniform mesh
of size $100\times 256$ is constructed in computational domain
$[0,1]\times[-1,1]$, yielding $N_{ic}=256$ initial points and $N_{r}=25600$
collocation points for enforcing the PDE residual. Fourier expansion of the
input coordinates in the embedding is determined by the first 20 spatial
harmonics ($m=20$). For this problem, the logarithm of (16) was used as the
loss function for training $u_{\boldsymbol{\theta}}$. We have chosen such an
approach because based on our experiments it demonstrates faster and more
stable convergence. The optimization problem in this case of loss function is

$\displaystyle{\bm{\theta}}^{*}=\underset{{\bm{\theta}}}{\arg}\min\log\left[{\mathcal{L}}(\boldsymbol{\theta})\right].$
(63)

The reference solution can be generated using the Chebfun package [55], or
dataset from supplementary materials for [19] can be used.

The experiments were carried out for various weights described by equations
(38), (45) and (50). Hyperparameters used for training are presented in table
5 and table 6. In table 1 we also report the accuracy for Allen–Cahn problem
achieved using other methods from the literature [18, 43, 56, 57, 19, 20, 51].
As can be seen from this table the best result is achieved with Dirac delta
function and MLP with a relative ${\mathbb{L}}_{2}$ error $6.29\times
10^{-5}$. Results of this experiment are shown on the Figure 2. The predicted
solution coincides with the ground truth almost perfectly.

Figure 2: Allen–Cahn equation. (a) is reference solution, (b) is prediction of a trained physics-informed neural network based on algorithm 1, (c) is absolute difference between reference solution and predicted solution. The relative error $\epsilon_{\rm error}$ is $6.29\times 10^{-5}$. (d), (e) and (f) are comparison of the predicted (red dash lines) and reference solutions (blue solid lines) corresponding to three temporal snapshots at $t=0.0$, $t=0.5$ and $t=1.0$, respectively. Method  | Relative ${\mathbb{L}}_{2}$ error
---|---
Original formulation by Raissi et al. [18] | $4.98\times 10^{-1}$
Adaptive time sampling [43] | $2.33\times 10^{-2}$
Self-attention [56] | $2.10\times 10^{-2}$
Time marching [57] | $1.68\times 10^{-2}$
Causal training (MLP) [19] | $1.43\times 10^{-3}$
Modified MLP with conservation laws [20] | $5.42\times 10^{-4}$
Causal training (weights (38))  | $\bf{5.08\times 10^{-4}}$
Heaviside function causal training  |
(weights (45)) | $\bf{3.34\times 10^{-4}}$
Causal training (modified MLP) [19] | $1.39\times 10^{-4}$
DASA-PINN [51] | $8.57\times 10^{-5}$
Dirac delta function causal training  |
(weights (50)) | $\bf{6.29\times 10^{-5}}$
Table 1: Allen–Cahn equation: Relative ${\mathbb{L}}_{2}$ errors obtained by
different approaches

###  4.2 Korteweg–De Vries equation

Consider the one-dimensional Korteweg–De Vries equation

$\displaystyle u_{t}+\eta uu_{x}+\mu^{2}u_{xxx}=0,\quad t\in[0,1],\quad
x\in[-1,1],$ (64) $\displaystyle u(x,0)=\cos(\pi x),$ (65) $\displaystyle
u(t,-1)=u(t,1).$ (66)

We used the classical parameters of this problem $\eta=1$ and $\mu=0.022$
[58]. We represent the latent variable $u$ by a fully-connected neural network
$u_{\boldsymbol{\theta}}$ with $tanh$ activation function, 3 hidden layers and
128 neurons per hidden layer. A uniform mesh of size $100\times 512$ was
constructed in computational domain $[0,1]\times[-1,1]$, yielding $N_{ic}=512$
initial points and $N_{r}=51200$ collocation points for enforcing the PDE
residual at the algorithm 1 with weights given by (38) and weights base on
[19]. For the algorithm 2 with weights (50) the mesh is of size $100\times
256$ with $N_{r}=25600$ collocation points. $\lambda_{ic}=200$,
$\lambda_{r}=5$ are used in vanilla causal training [19]. For all scenarios
the number of training epochs is equal to $300000$. The reference solution was
generated using the pseudo-spectral method [58].

To apply algorithm 2 with weights (50) we need to derive the functions
$\mathcal{R}\left[{\bf U}^{*}_{{\bm{\theta}}}\right]$ for current problem.
This derivation is presented in the next subsection.

####  4.2.1 Spatial-temporal local training

To use the algorithm 2 we need to represent the value $\mathcal{R}\left[{\bf
U}^{*}_{{\bm{\theta}}}\right]$ in terms of $\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]$ for $t\geq 0$. For this purpose let use expression
(30).

$\displaystyle\mathcal{R}\left[{\bf
U}^{*}_{{\bm{\theta}}}\right]=\frac{\partial{\bf
U}^{*}_{{\bm{\theta}}}}{\partial t}+\eta{\bf
U}^{*}_{{\bm{\theta}}}\frac{\partial{\bf U}^{*}_{{\bm{\theta}}}}{\partial
x}+\mu^{2}\frac{\partial^{3}{\bf U}^{*}_{{\bm{\theta}}}}{\partial x^{3}}.$
(67)

$\displaystyle\mathcal{R}\left[{\bf U}^{*}_{{\bm{\theta}}}\right]$
$\displaystyle=\left(\frac{\partial{\bf u}_{{\bm{\theta}}}}{\partial
t}+\eta{\bf u}_{{\bm{\theta}}}\frac{\partial{\bf u}_{{\bm{\theta}}}}{\partial
x}+\mu^{2}\frac{\partial^{3}{\bf u}_{{\bm{\theta}}}}{\partial
x^{3}}\right)\left[\chi(x-X_{1})-\chi(x-X_{2})\right]$
$\displaystyle+\left[\frac{\partial{\bf f_{1}}(t)}{\partial
t}\chi(X_{1}-x)+\frac{\partial{\bf f_{2}}(t)}{\partial t}\chi(x-X_{2})\right]$
$\displaystyle+\eta{\bf f_{1}}(t)\left[{\bf u}_{{\bm{\theta}}}(t,{x})-{\bf
f_{1}}(t)\right]\delta(x-X_{1})$ $\displaystyle+\eta{\bf f_{2}}(t)\left[{\bf
u}_{{\bm{\theta}}}(t,{x})-{\bf f_{2}}(t)\right]\delta(x-X_{2})$
$\displaystyle+3\mu^{2}\left[\frac{\partial{\bf u}_{{\bm{\theta}}}}{\partial
x}\frac{\partial\delta(x-X_{1})}{\partial x}+\frac{\partial^{2}{\bf
u}_{{\bm{\theta}}}}{\partial x^{2}}\delta(x-X_{1})\right]$
$\displaystyle-3\mu^{2}\left[\frac{\partial{\bf u}_{{\bm{\theta}}}}{\partial
x}\frac{\partial\delta(x-X_{2})}{\partial x}+\frac{\partial^{2}{\bf
u}_{{\bm{\theta}}}}{\partial x^{2}}\delta(x-X_{2})\right]$
$\displaystyle+\mu^{2}\left[{\bf u}_{{\bm{\theta}}}-{\bf
f_{1}}(t)\right]\frac{\partial^{2}\delta(x-X_{1})}{\partial
x^{2}}+\mu^{2}\left[{\bf f_{2}}(t)-{\bf
u}_{{\bm{\theta}}}\right]\frac{\partial^{2}\delta(x-X_{2})}{\partial x^{2}}.$
(68)

By applying the fundamental equation that defines the derivatives of the delta
function [50]

$\int_{-\infty}^{\infty}f(x)\frac{\partial^{n}\delta(x)}{\partial
x^{n}}dx=-\int_{-\infty}^{\infty}\frac{\partial f(x)}{\partial
x}\frac{\partial^{n-1}\delta(x)}{\partial x^{n-1}}dx,$ (69)

we can simplify equation (68) to form

$\displaystyle\mathcal{R}\left[{\bf U}^{*}_{{\bm{\theta}}}\right]$
$\displaystyle=\left(\frac{\partial{\bf u}_{{\bm{\theta}}}}{\partial
t}+\eta{\bf u}_{{\bm{\theta}}}\frac{\partial{\bf u}_{{\bm{\theta}}}}{\partial
x}+\mu^{2}\frac{\partial^{3}{\bf u}_{{\bm{\theta}}}}{\partial
x^{3}}\right)\left[\chi(x-X_{1})-\chi(x-X_{2})\right]$
$\displaystyle+\left[\frac{\partial{\bf f_{1}}(t)}{\partial
t}\chi(X_{1}-x)+\frac{\partial{\bf f_{2}}(t)}{\partial t}\chi(x-X_{2})\right]$
$\displaystyle+\eta{\bf f_{1}}(t)\left[{\bf u}_{{\bm{\theta}}}(t,{x})-{\bf
f_{1}}(t)\right]\delta(x-X_{1})$ $\displaystyle+\eta{\bf f_{2}}(t)\left[{\bf
u}_{{\bm{\theta}}}(t,{x})-{\bf f_{2}}(t)\right]\delta(x-X_{2})$
$\displaystyle+\mu^{2}\frac{\partial^{2}{\bf u}_{{\bm{\theta}}}}{\partial
x^{2}}\left[\delta(x-X_{1})-\delta(x-X_{2})\right].$ (70)

Expression (29) in the domain $t\geq 0$ has the following representation

$\displaystyle\mathcal{R}^{*}\left[{\bf U}_{{\bm{\theta}}},{\bf
g}\right]\left(t,{x}\right)$ $\displaystyle{=}\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t,{x}\right){+}\left[{\bf
u}_{{\bm{\theta}}}\left(t,{x}\right){-}{\bf
g}\left({x}\right)\right]\delta(t)$ $\displaystyle+\eta{\bf
f_{1}}(t)\left[{\bf u}_{{\bm{\theta}}}(t,{x})-{\bf
f_{1}}(t)\right]\delta(x-X_{1})$ $\displaystyle+\eta{\bf f_{2}}(t)\left[{\bf
u}_{{\bm{\theta}}}(t,{x})-{\bf f_{2}}(t)\right]\delta(x-X_{2})$
$\displaystyle+\mu^{2}\frac{\partial^{2}{\bf u}_{{\bm{\theta}}}}{\partial
x^{2}}\left[\delta(x-X_{1})-\delta(x-X_{2})\right].$ (71)

For the periodic boundary conditions with minimum spatial period $X_{2}-X_{1}$
equation (71) is simplified to

$\displaystyle\mathcal{R}^{*}\left[{\bf U}_{{\bm{\theta}}},{\bf
g}\right]\left(t,{x}\right)$ $\displaystyle{=}\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t,{x}\right){+}\left[{\bf
u}_{{\bm{\theta}}}\left(t,{x}\right){-}{\bf
g}\left({x}\right)\right]\delta(t)$ $\displaystyle+2\eta{\bf
u}_{{\bm{\theta}}}(t,X_{2})\left[{\bf u}_{{\bm{\theta}}}(t,{x})-{\bf
u}_{{\bm{\theta}}}(t,X_{2})\right]\delta(x-X_{1}).$ (72)

In this case the residual loss (17) is given by

$\displaystyle\mathcal{L}_{r}\left(\bm{\theta}\right)$
$\displaystyle=\frac{1}{T}\frac{1}{V_{\Omega}}\int^{T}_{0}dt\int_{X_{1}}^{X_{2}}d{x}\left|\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t,{x}\right)\right|^{2}$
$\displaystyle+\frac{1}{T}\frac{1}{V_{\Omega}}\left\\{2\int_{X_{1}}^{X_{2}}d{x}\left|\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(0,{x}\right)\left[{\bf
u}_{{\bm{\theta}}}\left(0,{x}\right){-}{\bf
g}\left({x}\right)\right]\right|\right.$
$\displaystyle\left.+\lim_{\begin{subarray}{c}\sigma\to 0\\\ t\to
0\end{subarray}}\frac{\sigma}{t^{2}+\sigma^{2}}\int_{X_{1}}^{X_{2}}d{x}\left|{\bf
u}_{{\bm{\theta}}}\left(0,{x}\right){-}{\bf
g}\left({x}\right)\right|^{2}\right\\}$
$\displaystyle+\frac{1}{T}\frac{1}{V_{\Omega}}\left\\{4\eta\int_{0}^{T}d{t}\left|\mathcal{R}\left[{\bf
u}_{{\bm{\theta}}}\right]\left(t,X_{1}\right){\bf
u}_{{\bm{\theta}}}(t,X_{2})\left[{\bf u}_{{\bm{\theta}}}(t,X_{1})-{\bf
u}_{{\bm{\theta}}}(t,X_{2})\right]\right|\right.$
$\displaystyle+\lim_{\begin{subarray}{c}\sigma\to 0\\\ x\to
X_{1}\end{subarray}}\frac{\sigma}{x^{2}+\sigma^{2}}4\eta^{2}\int_{0}^{T}d{t}\left|{\bf
u}_{{\bm{\theta}}}(t,X_{2})\left[{\bf u}_{{\bm{\theta}}}(t,x)-{\bf
u}_{{\bm{\theta}}}(t,X_{2})\right]\right|^{2}$
$\displaystyle\left.+4\eta^{2}\lim_{\begin{subarray}{c}\sigma\to 0\\\ x\to
X_{1}\end{subarray}}\frac{\sigma}{x^{2}+\sigma^{2}}\lim_{\begin{subarray}{c}\tilde{\sigma}\to
0\\\ t\to
0\end{subarray}}\frac{\tilde{\sigma}}{t^{2}+\tilde{\sigma}^{2}}\left|{\bf
u}_{{\bm{\theta}}}(t,X_{2})\left[{\bf u}_{{\bm{\theta}}}(t,x)-{\bf
u}_{{\bm{\theta}}}(t,X_{2})\right]\right|\left|{\bf
u}_{{\bm{\theta}}}\left(t,{x}\right){-}{\bf
g}\left({x}\right)\right|\right\\}.$ (73)

Figure 3: Korteweg–De Vries equation. (a) is reference solution, (b) is prediction of a trained physics-informed neural network based on algorithm 2, (c) is the absolute difference between reference solution and predicted one. The relative error $\epsilon_{\rm error}$ is $6.84\times 10^{-3}$. (d), (e) and (f) are comparison of the predicted (red dash lines) and reference solutions (blue solid lines) corresponding to the three temporal snapshots at $t=0.0$, $t=0.5$ and $t=1.0$, respectively. Method  | Relative ${\mathbb{L}}_{2}$ error
---|---
Causal training (algorithm 1 with weights (38)) | $3.11\times 10^{-2}$
Vanilla causal training (algorithm 1 with weights from [19]) | $2.03\times 10^{-2}$
Dirac delta function causal training |
(algorithm 2 with weights (50)) | ${6.84\times 10^{-3}}$
Dirac delta function causal training |
(algorithm 2 with weights (50) for the PAF MLP) | ${2.45\times 10^{-3}}$
Table 2: Korteweg–DeVries equation: Relative ${\mathbb{L}}_{2}$ errors
obtained by different approaches

Training hyperparameters are presented in the table 5 and table 6. Table 2
demonstrates the accuracy of Korteweg–De Vries problem solution for different
approaches. As can be seen the best result is achieved with Dirac delta
function and algorithm 2 for MLP with a resulting relative ${\mathbb{L}}_{2}$
error $6.84\times 10^{-3}$. Results of this experiment are shown on the Figure
3. The predicted solution coincides with the ground truth almost perfectly.

####  4.2.2 MLP with heterogeneous activation functions

Despite the relativity good results obtained for problem (64)–(66) with the
help of described fully connected neural network, the architecture discussed
above is not the best for this problem. As noted in [59], although the
universal approximation theorem [45] guarantees the existence of a neural
approximator for the given problem, it does not provide any mechanism for
choosing or constructing the most efficient neural network architecture for
PINN training.

Let use the approach described in [59], which allows to design the
architecture of a neural network which demonstrates better results for this
problem. This approach is named physical activation function PINN (PAF PINN)
which is based on the known solutions or behavior of unknown solution $u$ for
the system of differential equations under consideration.

The fundamental solution of a differential equation (64) in an unbounded
medium can be represented as a soliton [60, 61]

$\displaystyle u(x,t)=c\ {\rm sech}^{2}\left(ax-bt-d\right),$ (74)

where $a$, $b$, $c$ and $d$ are constants, which depend on the initial
condition and parameters $\eta$ and $\mu$, ${\rm
sech}(\cdot)=\cosh^{-1}(\cdot)$ is hyperbolic secant. In accordance with the
concept of PAF [59], the activation function for one of the last layers of the
neural network should be the function ${\rm sech}^{2}$.

In addition, we take into account that in the system described by equation
(64)–(66) there can be several solitons whose interactions satisfy the
nonlinear superposition principle [62]. The interactions of solitons are
defined by Bäcklund transform which include $tanh$ function for solitons of
the limited values. Therefore, as the activation function of the last layer,
it is necessary to take the hyperbolic tangent.

Taking into account the noted features of the known solutions of equation
(64), the following neural network architecture was proposed:

  1. 1.

the input layer consists of 2 neurons (for $x$ and $t$ data);

  2. 2.

the second layer consists of 120 neurons with a linear activation function;

  3. 3.

the third layer consists of 128 neurons of which 10 have a linear activation
function, and 118 neurons have a hyperbolic tangent as activation function;

  4. 4.

the fourth layer is the same as third;

  5. 5.

the fifth layer consists of 128 neurons of which 10 have a linear activation
function, and 118 neurons have ${\rm sech}^{2}$ as activation function;

  6. 6.

the sixth layer is the same as third;

  7. 7.

the seventh layer consists of one linear neuron.

A schematic representation of this neural network is shown in the panel (g) of
the figure 4.

Hyperparameters of Dirac delta function causal training for this neural
network are presented in the table 5 and table 6. The training consists of $3$
iterations with $3\times 10^{4}$ epochs for each of first two iteration, and
$9\times 10^{4}$ epochs for the last one. The relative error $\epsilon_{\rm
error}$ after training at the first, second and third iterations equals to
$8.45\times 10^{-3}$, $4.03\times 10^{-3}$ and $2.45\times 10^{-3}$,
respectively. Results of this experiment are shown on the Figure 4. As can be
seen the predicted solution coincides with the ground truth with great
accuracy.

Figure 4: Korteweg–De Vries equation. (a) is reference solution, (b) is
prediction of a trained physics-informed neural network based on the algorithm
2, (c) is the absolute difference between reference solution and predicted
solution. The relative error $\epsilon_{\rm error}$ is $2.45\times 10^{-3}$.
(d), (e) and (f) are comparison of the predicted (red dash lines) and
reference solutions (blue solid lines) corresponding to the three temporal
snapshots at $t=0.0$, $t=0.5$ and $t=1.0$, respectively. (g) is schematic
representation of the architecture of the neural network.

Based on the experiments we have come to the conclusion that the neural
network with heterogeneous activation functions not only has better accuracy
and speed of convergence compared with a neural network with homogeneous
activation functions, but also superiorly extrapolate the solution in space-
time domain where the training was not performed. Figure 5 shows the solutions
for the MLP (b) described in the previous subsection and for PAF MLP (e) for
time interval $t\in[0,2]$. It can be seen, PAF MLP demonstrates significantly
better qualitative behavior of the solution in the domain $t\in[1,2]$, in
comparison with the behavior of MLP in the same area. The relative error
quantitatively represents these results: PAF MLP has relative error
$\epsilon_{\rm error}=3.28\times 10^{-1}$, whereas MLP has relative error
$\epsilon_{\rm error}=6.61\times 10^{-1}$.

Figure 5: Korteweg–De Vries equation. (a) is reference solution, (b) is
prediction of a trained physics-informed neural network based on algorithm 1
(best result for MLP, $\epsilon_{\rm error}=6.84\times 10^{-3}$ for
$t\in[0,1]$ and $\epsilon_{\rm error}=6.61\times 10^{-1}$ for $t\in[0,2]$),
(c) is absolute value of difference of reference solution and predicted
solution. (d)–(f) are same for the PAF MLP ($\epsilon_{\rm error}=2.45\times
10^{-3}$ for $t\in[0,1]$ and $\epsilon_{\rm error}=3.28\times 10^{-1}$ for
$t\in[0,2]$).

###  4.3 Petrov–Kudrin equations

Consider the Petrov–Kudrin system of equations [63], which describes the
propagation of intense electromagnetic waves in a nonlinear nondispersive
medium,

$\displaystyle E_{\rho}-{\varepsilon_{1}^{-1/2}}H_{\tau}=0,$ $\displaystyle
H_{\rho}+{\rho^{-1}}H-\varepsilon_{1}^{1/2}\exp\left(\alpha
E\right)E_{\tau}=0,$ (75)

where $\rho\in[0,5]$, $\tau\in[0,4.75]$, $\alpha=1$, $\varepsilon_{1}=2$. We
will consider the initial problem with initial values

$\displaystyle E(0,\rho)=\left[1+\rho^{2}\exp\left(\alpha
E\right)\right]^{-3/2},\quad H(0,\rho)=0.$ (76)

This problem has the exact solution [63] in the following form

$\displaystyle E(\tau,\rho)={\rm
Re}\left\\{\left[\left(1-i\theta\right)^{2}+\rho^{2}\exp\left(\alpha
E\right)\right]^{-3/2}\left(1-i\theta\right)\right\\},$ $\displaystyle
H(\tau,\rho)=\varepsilon_{1}^{1/2}\rho\exp\left(\alpha E\right){\rm
Re}\left\\{i\left[\left(1-i\theta\right)^{2}+\rho^{2}\exp\left(\alpha
E\right)\right]^{-3/2}\right\\}.$ (77)

Here $\theta=\tau+\alpha\rho H/(2\varepsilon_{1}^{1/2})$.

We represent the latent variables $E$ and $H$ with the help of a fully-
connected neural network $u_{\boldsymbol{\theta}}$ with $tanh$ activation
function, 6 hidden layers and 64 neurons per hidden layer. For simplicity, a
uniform mesh of size $100\times 256$ was constructed in the computational
domain $[0,4.75]\times[0,5]$, yielding $N_{ic}=512$ initial points and
$N_{r}=25600$ collocation points for enforcing the PDE residual in algorithm 1
with weights (38), (50) or with weights from [19]. For all scenarios the
number of epochs is supposed to be $300000$. As in the case of the Allen–Cahn
problem, for this problem, the logarithm of (16) and the optimization problem
63 were used for training $u_{\boldsymbol{\theta}}$.

Method  | Relative ${\mathbb{L}}_{2}$ error
---|---
Vanilla causal training (algorithm 1 with weights from [19]) | $8.01\times 10^{-3}$
Causal training (algorithm 1 with weights (38)) | $7.82\times 10^{-3}$
Dirac delta function causal training |
(algorithm 1 with weights (50)) | ${7.69\times 10^{-3}}$
Table 3: Petrov–Kudrin equations: Relative ${\mathbb{L}}_{2}$ errors obtained
by different approaches

To use the algorithm 1 it is necessary to get away from the limit
representation of delta function towards the discrete form. The required
transformation is described in the following subsection.

####  4.3.1 Approximation of Dirac delta function

According to the equation (22) let approximate the delta function in the (23)
with function $1/\sigma$. It is assumed that the $\sigma$ value is the doubled
time during which the initial distribution of the field decreases or increases
twice. For the problem (75),(76) we need to find dependence of field $E$ on
$\tau$ at $\rho=0$. This value is further denoted with $\tilde{E}$. In
addition, we suppose that the field $E$ depends on the coordinate $\rho$
similar to $E$ in (76). Taking into account the latter circumstance, we obtain
the following equation for the $\tilde{E}$ from the equations (75)

$\frac{\partial}{\partial\tau}\left(e^{\alpha\tilde{E}}\frac{\partial\tilde{E}}{\partial\tau}\right)=-6e^{\alpha\tilde{E}}.$
(78)

The solution of the equation (78) is

$\tilde{E}(\tau)=\frac{1}{\alpha}\log\left[e^{\alpha}\cos\left(\tau\sqrt{6\alpha}\right)\right].$
(79)

Here we used initial value $\tilde{E}=1$. It follows from the equation (79)
that the decay time of the field $\tilde{E}$ is two times

$\tilde{\tau}=\arccos\left(e^{-\alpha/2}\right)/\sqrt{6\alpha}.$ (80)

It follows from this equation that

$\delta(\tau)\simeq 1/(2\tilde{\tau}).$ (81)

Figure 6: Petrov–Kudrin equations: the first row is $E$, the second row is
$H$. (a) and (d) are reference solutions, (b) and (e) are predictions of a
trained physics-informed neural network based on the algorithm 1, (c) and (f)
are absolute values of difference between reference and predicted solutions.
The total relative error $\epsilon_{\rm error}$ is $7.69\times 10^{-3}$.
Figure 7: Petrov–Kudrin equations: (a), (b) and (c) are comparison of the
predicted (red dash lines) and reference solutions (blue solid lines) of
electric field $E$ corresponding to the three temporal snapshots at
$\tau=0.0$, $\tau=2.375$ and $\tau=4.75$, respectively. (d), (e) and (f) are
the same dependences for the magnetic field $H$

Training hyperparameters for algorithm 1 with weights from [19], with weights
(38) and with weights (50) are presented in the table 5 and table 6. Table 3
presents the accuracy of the solution for Petrov–Kudrin problem obtained with
the help of these approaches. The best result is achieved with Dirac delta
function and algorithm 1 for MLP with relative ${\mathbb{L}}_{2}$ error
$7.69\times 10^{-3}$. Results of this experiment are shown on the Figures 6
and 7. The predicted solution coincides with the ground truth with great
accuracy.

##  5 Conclusion

In this article, we presented several approaches to neural network training
within the framework of the PINN concept. These approaches simplify the
construction and analysis of the loss function and improve the training
convergence. We have also proposed the neural network architectures which are
more relevant for the problems under consideration.

We have reformulated the original problem described by the differential
equation and the initial and boundary conditions into the problem described
only by the differential equation. The major advantage of such a trick is that
it becomes possible to reduce the loss function to the single term associated
with the differential equations, thus eliminating the need to tune the scaling
coefficients for the terms associated with the boundary and initial
conditions. Based on this approach we have derived closed-form expressions for
the MAE and MSE losses. Approximations of the Dirac delta function, which is
part of the losses in the MSE formulation, have been proposed. These
approximations are determined by the discretization method of the domain or by
the behavior of the expected solution at the initial time and at the
boundaries of the domain.

Inspired by [19] we have proposed loss functions based on causality and
spatial-temporal locality principles within the framework of the formalism of
generalized functions such as Heaviside step and Dirac delta functions. It has
also been shown that the causality training can be extended to the spatial
domains and can be transformed into the spatial-temporal training.

Finally, a modification of the MLP model for the Korteweg–De Vries problem has
been proposed based on [59]. This modification is based on the fundamental
solution of the Korteweg–De Vries equation and the nonlinear superposition
principle for the solitons. It has been shown that such a neural network with
heterogeneous activation functions can be trained much faster and the obtained
results have significantly better accuracy compared to MLP with homogeneous
activation functions. In addition, such neural network is much better able to
extrapolate to new domains.

Numerical experiments have been performed for a number of problems and the
accuracies of the proposed methods are given. The causal trainings with
weights based on the Dirac delta function were achieved the best results in
all carried out the numerical experiments.

##  Appendix A Nomenclature

The summarizes the main notations, abbreviations and symbols are given in
table 4.

Notation | Description
---|---
MLP | Multilayer perceptron
PDE | Partial differential equation
PINN | Physic-informed neural network
MAE | Mean absolute error
MSE | Mean squared error
DASA-PINN | Differentiable adversarial self-adaptive pointwise weighing scheme for PINN
PAF | physical activation function
${\bf u}$ | solution of PDE
$\bf U$ | artificially extended solution of PDE
${\mathcal{N}}$ | nonlinear differential operator
${\mathcal{B}}$ | boundary operator
$\mathcal{R}$ | PDE residual
${\bf u}_{\boldsymbol{\theta}}$ | neural network representation of the PDE solution
${\boldsymbol{\theta}}$ | vector of the trainable parameters of the neural network
$N_{t}$ | number of intervals mesh along the time axis
$N_{x}$ | number of intervals mesh along the $x$ axis
$w_{i}$ | residual weight at the time $t_{i}$
$w_{n}^{(\to)}$ | residual weight at the coordinate $x_{n}$ in case of the beginning calculation of weight at low boundary
$w_{n}^{(\leftarrow)}$ | residual weight at the coordinate $x_{n}$ in case of the beginning calculation of weight at upper boundary
$\varepsilon$ | causality parameter
$\delta_{w}$ | threshold for increasing a causality parameter
${\mathcal{L}\left(t,{\boldsymbol{\theta}}\right)}$ | temporal residual loss
${\mathcal{L}^{*}\left(t,{\boldsymbol{\theta}}\right)}$ | temporal residual loss in the spatial-temporal domain
${\mathcal{L}(\boldsymbol{\theta})}$ | aggregate training loss
Table 4: Nomenclature

##  Appendix B Hyperparameters of approaches

Tables 5 and 6 summarize the training and network hyperparameters,
respectively, for all numerical experiments.

The optimizer used in all experiments is Adam, and parameters $\theta$ are
initialized with the Xavier scheme [64].

PDE Problem  | Weights | $N$ | $\eta_{s}$ | $\eta_{\min}$ | $\varepsilon_{\rm int}$ | $\delta_{w}$ | Scheduler
---|---|---|---|---|---|---|---
Allen–Cahn | (38) | $3\times 10^{5}$ | $3\times 10^{-3}$ | $10^{-12}$ | $10^{-3}$ | $0.99$ | ExponentialLR
(51) | $10^{-5}$ | $0.95$
(50) | $6\times 10^{5}$ | $10^{-3}$ | — | $10^{4}$ | $0.99$ | StepLR${}^{0.95}_{5000}$
Korteweg–De Vries | [19] | $3\times 10^{5}$ | $3\times 10^{-3}$ | — | $10^{-2}$ | $0.99$ | StepLR${}^{0.9}_{5000}$
(38)
(50)
Korteweg–De Vries (PAF) | (50) | $3\times 10^{4}$ | $3\times 10^{-3}$ | $10^{-8}$ | $10^{-5}$ | $0.99$ | CosineAnnealingLR
$10^{-4}$ | $0.16$
$9\times 10^{4}$ | $3\times 10^{-5}$ | $0.33$
Petrov–Kudrin | [19] | $3\times 10^{5}$ | $10^{-3}$ | — | $10^{-2}$ | $0.99$ | StepLR${}^{0.9}_{5000}$
(38)
(50)
Table 5: Training hyperparameters

Here $N$ is number of epochs, $\eta_{s}$ is initial value of learning rate,
$\eta_{\min}$ is minimum value of learning rate, $\varepsilon_{\rm int}$ is
the initial value of the causality parameter. For the ExponentialLR scheduler
the a decay–rate is calculated with following formula

$\gamma=\left(\eta_{\min}/\eta_{s}\right)^{1/N}.$

StepLR${}^{\gamma}_{E}$ means decay–rate is $\gamma$ every $E$ training
epochs.

PDE Problem  | Architecture |  Hidden layers |  Neurons per layer |  Activation function | $N_{t}$ | $N_{x}$
---|---|---|---|---|---|---
Allen–Cahn | MLP | 4 | 128 | $\tanh$ | 100 | 256
Korteweg–De Vries | MLP | 3 | 128 | $\tanh$ | 100 | 512 for [19] and (38)
|  |  |  |  |  | 256 for (50)
Korteweg–De Vries | PAF MLP | 5 | 120 or 128 | linear, $\tanh$ or ${\rm sech}^{2}$ | 100 | 512 for (50)
Petrov–Kudrin | MLP | 6 | 64 | $\tanh$ | 100 | 256
Table 6: The models and mesh parameters

Here $N_{t}$ and $N_{x}$ are number of values along time and $x$ coordinate
axes, respectively.

## References

  * [1] L. Alzubaidi, J. Zhang, A. J. Humaidi, A. A. Dujaili, Y. Duan, O. A. Shamma, J. Santamaría, M. A. Fadhel, M. A. Amidie, and L. Farhan, _Review of deep learning: concepts, CNN architectures, challenges, applications, future directions_. Springer International Publishing, 2021. [Online]. Available: https://doi.org/10.1186/s40537-021-00444-8
  * [2] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _Nature_ , vol. 521, no. 7553, pp. 436–444, may 2015. [Online]. Available: http://www.nature.com/articles/nature14539
  * [3] P. Linardatos, V. Papastefanopoulos, and S. Kotsiantis, “Explainable AI: A Review of Machine Learning Interpretability Methods,” _Entropy_ , vol. 23, no. 1, 2021. [Online]. Available: https://www.mdpi.com/1099-4300/23/1/18
  * [4] S. Khan, M. Naseer, M. Hayat, S. W. Zamir, F. S. Khan, and M. Shah, “Transformers in Vision: A Survey,” _ACM Comput. Surv._ , vol. 54, no. 10s, sep 2022. [Online]. Available: https://doi.org/10.1145/3505244
  * [5] A. Vaswani, S. Bengio, E. Brevdo, F. Chollet, A. N. Gomez, S. Gouws, L. Jones, L. Kaiser, N. Kalchbrenner, N. Parmar, R. Sepassi, N. Shazeer, and J. Uszkoreit, “Tensor2Tensor for Neural Machine Translation,” 2018. [Online]. Available: https://arxiv.org/abs/1803.07416
  * [6] Q. Wang, B. Li, T. Xiao, J. Zhu, C. Li, D. F. Wong, and L. S. Chao, “Learning Deep Transformer Models for Machine Translation,” 2019. [Online]. Available: https://arxiv.org/abs/1906.01787
  * [7] S. Yao and X. Wan, “Multimodal transformer for multimodal machine translation,” in _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_. Online: Association for Computational Linguistics, Jul. 2020, pp. 4346–4350. [Online]. Available: https://aclanthology.org/2020.acl-main.400
  * [8] A. Ramesh, M. Pavlov, G. Goh, S. Gray, C. Voss, A. Radford, M. Chen, and I. Sutskever, “Zero-Shot Text-to-Image Generation,” 2021. [Online]. Available: https://arxiv.org/abs/2102.12092
  * [9] S. Frolov, T. Hinz, F. Raue, J. Hees, and A. Dengel, “Adversarial text-to-image synthesis: A review,” _Neural Networks_ , vol. 144, pp. 187–209, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0893608021002823
  * [10] D. Silver, T. Hubert, J. Schrittwieser, I. Antonoglou, M. Lai, A. Guez, M. Lanctot, L. Sifre, D. Kumaran, T. Graepel, T. Lillicrap, K. Simonyan, and D. Hassabis, “A general reinforcement learning algorithm that masters chess, shogi, and Go through self-play,” _Science_ , vol. 362, no. 6419, pp. 1140–1144, 2018. [Online]. Available: https://www.science.org/doi/abs/10.1126/science.aar6404
  * [11] J. Schrittwieser, I. Antonoglou, T. Hubert, K. Simonyan, L. Sifre, S. Schmitt, A. Guez, E. Lockhart, D. Hassabis, T. Graepel, T. Lillicrap, and D. Silver, “Mastering Atari, Go, chess and shogi by planning with a learned model,” _Nature_ , vol. 588, no. 7839, pp. 604–609, dec 2020. [Online]. Available: http://www.nature.com/articles/s41586-020-03051-4
  * [12] L. Ouyang, J. Wu, X. Jiang, D. Almeida, C. L. Wainwright, P. Mishkin, C. Zhang, S. Agarwal, K. Slama, A. Ray, J. Schulman, J. Hilton, F. Kelton, L. Miller, M. Simens, A. Askell, P. Welinder, P. Christiano, J. Leike, and R. Lowe, “Training language models to follow instructions with human feedback,” 2022\. [Online]. Available: https://arxiv.org/abs/2203.02155
  * [13] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, “Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators,” _Nature Machine Intelligence_ , vol. 3, no. 3, pp. 218–229, mar 2021. [Online]. Available: https://doi.org/10.1038%2Fs42256-021-00302-5
  * [14] Z. Li, D. Z. Huang, B. Liu, and A. Anandkumar, “Fourier Neural Operator with Learned Deformations for PDEs on General Geometries,” 2022. [Online]. Available: https://arxiv.org/abs/2207.05209
  * [15] M. A. Krinitskiy, V. M. Stepanenko, A. O. Malkhanov, and M. E. Smorkalov, “A General Neural-Networks-Based Method for Identification of Partial Differential Equations, Implemented on a Novel AI Accelerator,” _Supercomputing Frontiers and Innovations_ , vol. 9, no. 3, sep 2022. [Online]. Available: https://superfri.org/index.php/superfri/article/view/439
  * [16] V. Fanaskov and I. Oseledets, “Spectral Neural Operators,” 2022. [Online]. Available: https://arxiv.org/abs/2205.10573
  * [17] O. Ovadia, A. Kahana, P. Stinis, E. Turkel, and G. E. Karniadakis, “ViTO: Vision Transformer-Operator,” mar 2023. [Online]. Available: http://arxiv.org/abs/2303.08891
  * [18] M. Raissi, P. Perdikaris, and 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_ , vol. 378, pp. 686–707, 2019. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0021999118307125
  * [19] S. Wang, S. Sankaran, and P. Perdikaris, “Respecting causality is all you need for training physics-informed neural networks,” 2022. [Online]. Available: http://arxiv.org/abs/2203.07404
  * [20] H. Wang, X. Qian, Y. Sun, and S. Song, “A Modified Physics Informed Neural Networks for Solving the Partial Differential Equation with Conservation Laws,” _—_. [Online]. Available: https://ssrn.com/abstract=4274376
  * [21] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations,” no. Part I, pp. 1–22. [Online]. Available: https://arxiv.org/abs/1711.10561
  * [22] ——, “Physics Informed Deep Learning (Part II): Data-driven Discovery of Nonlinear Partial Differential Equations,” no. Part II, pp. 1–19. [Online]. Available: https://arxiv.org/abs/1711.10566
  * [23] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko, A. Bridgland, C. Meyer, S. A. A. Kohl, A. J. Ballard, A. Cowie, B. Romera-Paredes, S. Nikolov, R. Jain, J. Adler, T. Back, S. Petersen, D. Reiman, E. Clancy, M. Zielinski, M. Steinegger, M. Pacholska, T. Berghammer, S. Bodenstein, D. Silver, O. Vinyals, A. W. Senior, K. Kavukcuoglu, P. Kohli, and D. Hassabis, “Highly accurate protein structure prediction with AlphaFold,” _Nature_ , vol. 596, no. 7873, pp. 583–589, aug 2021. [Online]. Available: https://www.nature.com/articles/s41586-021-03819-2
  * [24] J. Lin, F. Haberstroh, S. Karsch, and D. Andreas, “Applications of object detection networks at high-power laser systems and experiments,” pp. 1–10, 2022\. [Online]. Available: http://arxiv.org/abs/2210.02539
  * [25] M. J. V. Streeter, C. Colgan, C. C. Cobo, C. Arran, E. E. Los, R. Watt, N. Bourgeois, L. Calvin, J. Carderelli, N. Cavanagh, and et al., “Laser wakefield accelerator modelling with variational neural networks,” _High Power Laser Science and Engineering_ , vol. 11, p. e9, 2023. [Online]. Available: https://www.cambridge.org/core/journals/high-power-laser-science-and-engineering/article/laser-wakefield-accelerator-modelling-with-variational-neural-networks/272206304D99ABB614B6F5A57B389B29
  * [26] D. C. Psichogios and L. H. Ungar, “A hybrid neural network-first principles approach to process modeling,” _Aiche Journal_ , vol. 38, pp. 1499–1511, 1992.
  * [27] I. Lagaris, A. Likas, and D. Fotiadis, “Artificial neural networks for solving ordinary and partial differential equations,” _IEEE Transactions on Neural Networks_ , vol. 9, no. 5, pp. 987–1000, 1998. [Online]. Available: https://doi.org/10.1109%2F72.712178
  * [28] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and A. Edelman, “Universal Differential Equations for Scientific Machine Learning,” pp. 1–55, jan 2020. [Online]. Available: http://arxiv.org/abs/2001.04385
  * [29] L. Yuan, Y.-Q. Ni, X.-Y. Deng, and S. Hao, “A-PINN: Auxiliary physics informed neural networks for forward and inverse problems of nonlinear integro-differential equations,” _Journal of Computational Physics_ , vol. 462, p. 111260, aug 2022. [Online]. Available: https://doi.org/10.1016/j.jcp.2022.111260https://linkinghub.elsevier.com/retrieve/pii/S0021999122003229
  * [30] X. Jin, S. Cai, H. Li, and G. E. Karniadakis, “NSFnets (Navier-Stokes flow nets): Physics-informed neural networks for the incompressible Navier-Stokes equations,” _Journal of Computational Physics_ , vol. 426, p. 109951, 2021\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021999120307257
  * [31] L. Zhao, Z. Li, Z. Wang, B. Caswell, J. Ouyang, and G. E. Karniadakis, “Active- and transfer-learning applied to microscale-macroscale coupling to simulate viscoelastic flows,” _Journal of Computational Physics_ , vol. 427, p. 110069, 2021. [Online]. Available: https://doi.org/10.1016/j.jcp.2020.110069
  * [32] E. Kharazmi, Z. Zhang, and G. E. Karniadakis, “hp-VPINNs: Variational physics-informed neural networks with domain decomposition,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 374, p. 113547, 2021. [Online]. Available: https://doi.org/10.1016/j.cma.2020.113547
  * [33] L. Yang, X. Meng, and G. E. Karniadakis, “B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data,” _Journal of Computational Physics_ , vol. 425, p. 109913, 2021. [Online]. Available: https://doi.org/10.1016/j.jcp.2020.109913
  * [34] S. Cuomo, V. S. di Cola, F. Giampaolo, G. Rozza, M. Raissi, and F. Piccialli, “Scientific Machine Learning through Physics-Informed Neural Networks: Where we are and What’s next,” jan 2022. [Online]. Available: http://arxiv.org/abs/2201.05624
  * [35] G. Pang, M. D’Elia, M. Parks, and G. E. Karniadakis, “nPINNs: Nonlocal physics-informed neural networks for a parametrized nonlocal universal Laplacian operator. Algorithms and applications,” _Journal of Computational Physics_ , vol. 422, p. 109760, 2020. [Online]. Available: https://doi.org/10.1016/j.jcp.2020.109760
  * [36] R. G. Patel, I. Manickam, N. A. Trask, M. A. Wood, M. Lee, I. Tomas, and E. C. Cyr, “Thermodynamically consistent physics-informed neural networks for hyperbolic systems,” _Journal of Computational Physics_ , vol. 449, p. 110754, jan 2022. [Online]. Available: https://doi.org/10.1016/j.jcp.2021.110754https://linkinghub.elsevier.com/retrieve/pii/S0021999121006495
  * [37] S. Cai, Z. Mao, Z. Wang, M. Yin, and G. E. Karniadakis, “Physics-informed neural networks (PINNs) for fluid mechanics: a review,” _Acta Mechanica Sinica_ , vol. 37, no. 12, pp. 1727–1738, dec 2021. [Online]. Available: https://link.springer.com/10.1007/s10409-021-01148-1
  * [38] S. Thakur, M. Raissi, H. Mitra, and A. Ardekani, “Temporal Consistency Loss for Physics-Informed Neural Networks,” pp. 1–14, jan 2023. [Online]. Available: http://arxiv.org/abs/2301.13262
  * [39] B. Moseley, A. Markham, and T. Nissen-Meyer, “Solving the wave equation with physics-informed deep learning,” jun 2020. [Online]. Available: http://arxiv.org/abs/2006.11894
  * [40] G. Lin, P. Hu, F. Chen, X. Chen, J. Chen, J. Wang, and Z. Shi, “BINet: Learning to Solve Partial Differential Equations with Boundary Integral Networks,” pp. 1–27, oct 2021. [Online]. Available: http://arxiv.org/abs/2110.00352
  * [41] Q. He, D. Barajas-Solano, G. Tartakovsky, and A. M. Tartakovsky, “Physics-informed neural networks for multiphysics data assimilation with application to subsurface transport,” _Advances in Water Resources_ , vol. 141, p. 103610, jul 2020. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0309170819311649
  * [42] T. Qin, A. Beatson, D. Oktay, N. McGreivy, and R. P. Adams, “Meta-PDE: Learning to Solve PDEs Quickly Without a Mesh,” nov 2022. [Online]. Available: http://arxiv.org/abs/2211.01604
  * [43] C. L. Wight and J. Zhao, “Solving Allen-Cahn and Cahn-Hilliard Equations using the Adaptive Physics Informed Neural Networks,” jul 2020. [Online]. Available: http://arxiv.org/abs/2007.04542
  * [44] A. Daw, J. Bu, S. Wang, P. Perdikaris, and A. Karpatne, “Mitigating Propagation Failures in PINNs using Evolutionary Sampling,” jul 2022. [Online]. Available: http://arxiv.org/abs/2207.02338
  * [45] K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” _Neural Networks_ , vol. 2, no. 5, pp. 359–366, 1989. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0893608089900208
  * [46] A. Griewank and A. Walther, _Evaluating Derivatives_ , 2nd ed. Society for Industrial and Applied Mathematics, 2008\. [Online]. Available: https://epubs.siam.org/doi/abs/10.1137/1.9780898717761
  * [47] S. Wang, Y. Teng, and P. Perdikaris, “Understanding and Mitigating Gradient Flow Pathologies in Physics-Informed Neural Networks,” _SIAM Journal on Scientific Computing_ , vol. 43, no. 5, pp. A3055—-A3081, 2021. [Online]. Available: https://doi.org/10.1137/20M1318043
  * [48] S. Wang, X. Yu, and P. Perdikaris, “When and why PINNs fail to train: A neural tangent kernel perspective,” _Journal of Computational Physics_ , vol. 449, p. 110768, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S002199912100663X
  * [49] B. Davies, _Integral Transforms and Their Applications_ , 3rd ed. Springer New York, NY, 2008. [Online]. Available: https://link.springer.com/book/10.1007/978-1-4684-9283-5#bibliographic-information
  * [50] A. I. Saichev and W. A. Woyczynski, _Distributions in the Physical and Engineering Sciences, Volume 1: Distributional and Fractal Calculus, Integral Transforms and Wavelets_ , 1st ed. Birkhäuser Cham, 1996. [Online]. Available: https://link.springer.com/book/10.1007/978-3-319-97958-8
  * [51] G. Zhang, H. Yang, F. Zhu, Y. Chen, and X. Zheng, “Dasa-Pinns: Differentiable Adversarial Self-Adaptive Pointwise Weighting Scheme for Physics-Informed Neural Networks,” _SSRN Electronic Journal_ , 2023. [Online]. Available: https://www.ssrn.com/abstract=4376049
  * [52] S. Dong and N. Ni, “A method for representing periodic functions and enforcing exactly periodic boundary conditions with deep neural networks,” _Journal of Computational Physics_ , vol. 435, p. 110242, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021999121001376
  * [53] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” 2014. [Online]. Available: https://arxiv.org/abs/1412.6980
  * [54] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, _PyTorch: An Imperative Style, High-Performance Deep Learning Library_. Red Hook, NY, USA: Curran Associates Inc., 2019.
  * [55] T. A. Driscoll, N. Hale, and L. N. Trefethen, _Chebfun Guide_. Pafnuty Publications, Oxford, 2014.
  * [56] L. McClenny and U. Braga-Neto, “Self-Adaptive Physics-Informed Neural Networks using a Soft Attention Mechanism,” sep 2020. [Online]. Available: http://arxiv.org/abs/2009.04544
  * [57] R. Mattey and S. Ghosh, “A novel sequential method to train physics informed neural networks for Allen–Cahn and Cahn–Hilliard equations,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 390, p. 114474, feb 2022. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0045782521006939
  * [58] N. J. Zabusky and M. D. Kruskal, “Interaction of "Solitons" in a Collisionless Plasma and the Recurrence of Initial States,” _Phys. Rev. Lett._ , vol. 15, pp. 240–243, Aug 1965. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.15.240
  * [59] J. Abbasi and P. Ø. Andersen, “Physical Activation Functions (PAFs): An Approach for More Efficient Induction of Physics into Physics-Informed Neural Networks (PINNs).” [Online]. Available: https://arxiv.org/abs/2205.14630
  * [60] D. J. Korteweg and G. de Vries, “XLI. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves,” _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , vol. 39, no. 240, pp. 422–443, may 1895. [Online]. Available: https://www.tandfonline.com/doi/full/10.1080/14786449508620739
  * [61] A. F. Vakakis, Ed., _Normal Modes and Localization in Nonlinear Systems_. Dordrecht: Springer Netherlands, 2001. [Online]. Available: http://link.springer.com/10.1007/978-94-017-2452-4
  * [62] K. Brauer, “The Korteweg-de Vries Equation: History, exact Solutions, and graphical Representation,” 2000. [Online]. Available: https://www.math.arizona.edu/~gabitov/teaching/141/math_485/KDV.pdf
  * [63] E. Y. Petrov and A. V. Kudrin, “Exact Axisymmetric Solutions of the Maxwell Equations in a Nonlinear Nondispersive Medium,” _Phys. Rev. Lett._ , vol. 104, p. 190404, May 2010. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.104.190404
  * [64] X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” in _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , ser. Proceedings of Machine Learning Research, Y. W. Teh and M. Titterington, Eds., vol. 9. Chia Laguna Resort, Sardinia, Italy: PMLR, 13–15 May 2010, pp. 249–256. [Online]. Available: https://proceedings.mlr.press/v9/glorot10a.html