Researchers have discussed physics informed neural networks for a l...
Partial differential equations (PDEs) arise in many natural phenome...
In the small data regime, it is very easy for neural networks or ot...
The key point made here is that neural networks that can leverage p...
The universal function approximation theorem of neural networks - t...
The elegant thing about this - is that automatic differentiation is...
"This simple yet powerful construction allows us to tackle a wide r...
Burger's equation was first introduced in 1915 by Harry Bateman, an...
One very attractive feature of this method, is that the training da...
One of the biggest advantage of PINNs is that they are simple to co...
Besides specifying the models u(x, t) and f(x, t), learning the par...
Clearly in this simple example of Burger's equation, PINNs converge...
Physics Informed Deep Learning (Part I): Data-driven
Solutions of Nonlinear Partial Differential Equations
Maziar Raissi
1
, Paris Perdikaris
2
, and George Em Karniadakis
1
1
Division of Applied Mathematics, Brown University,
Providence, RI, 02912, USA
2
Department of Mechanical Engineering and Applied Mechanics,
University of Pennsylvania,
Philadelphia, PA, 19104, USA
Abstract
We introduce physics informed neural networks – neural networks that
are trained to solve supervised learning tasks while respecting any given law
of physics described by general nonlinear partial differential equations. In
this two part treatise, we present our developments in the context of solving
two main classes of problems: data-driven solution and data-driven discovery
of partial differential equations. Depending on the nature and arrangement
of the available data, we devise two distinct classes of algorithms, namely
continuous time and discrete time models. The resulting neural networks
form a new class of data-efficient universal function approximators that nat-
urally encode any underlying physical laws as prior information. In this first
part, we demonstrate how these networks can be used to infer solutions to
partial differential equations, and obtain physics-informed surrogate models
that are fully differentiable with respect to all input coordinates and free
parameters.
Keywords: Data-driven scientific computing, Machine learning, Predictive
modeling, Runge-Kutta methods, Nonlinear dynamics
1. Introduction
With the explosive growth of available data and computing resources, re-
cent advances in machine learning and data analytics have yielded transfor-
mative results across diverse scientific disciplines, including image recognition
[1], natural language processing [2], cognitive science [3], and genomics [4].
Preprint submitted to Journal Name November 30, 2017
arXiv:1711.10561v1 [cs.AI] 28 Nov 2017
However, more often than not, in the course of analyzing complex physical,
biological or engineering systems, the cost of data acquisition is prohibitive,
and we are inevitably faced with the challenge of drawing conclusions and
making decisions under partial information. In this small data regime, the
vast majority of state-of-the art machine learning techniques (e.g., deep/-
convolutional/recurrent neural networks) are lacking robustness and fail to
provide any guarantees of convergence.
At first sight, the task of training a deep learning algorithm to accurately
identify a nonlinear map from a few – potentially very high-dimensional –
input and output data pairs seems at best naive. Coming to our rescue, for
many cases pertaining to the modeling of physical and biological systems,
there a exist a vast amount of prior knowledge that is currently not being
utilized in modern machine learning practice. Let it be the principled physical
laws that govern the time-dependent dynamics of a system, or some empirical
validated rules or other domain expertise, this prior information can act as
a regularization agent that constrains the space of admissible solutions to a
manageable size (for e.g., in incompressible fluid dynamics problems by dis-
carding any non realistic flow solutions that violate the conservation of mass
principle). In return, encoding such structured information into a learning
algorithm results in amplifying the information content of the data that the
algorithm sees, enabling it to quickly steer itself towards the right solution
and generalize well even when only a few training examples are available.
The first glimpses of promise for exploiting structured prior information
to construct data-efficient and physics-informed learning machines have al-
ready been showcased in the recent studies of [5, 6, 7]. There, the authors
employed Gaussian process regression [8] to devise functional representations
that are tailored to a given linear operator, and were able to accurately infer
solutions and provide uncertainty estimates for several prototype problems
in mathematical physics. Extensions to nonlinear problems were proposed
in subsequent studies by Raissi et. al. [9, 10] in the context of both in-
ference and systems identification. Despite the flexibility and mathematical
elegance of Gaussian processes in encoding prior information, the treatment
of nonlinear problems introduces two important limitations. First, in [9, 10]
the authors had to locally linearize any nonlinear terms in time, thus limit-
ing the applicability of the proposed methods to discrete-time domains and
compromising the accuracy of their predictions in strongly nonlinear regimes.
2
Secondly, the Bayesian nature of Gaussian process regression requires certain
prior assumptions that may limit the representation capacity of the model
and give rise to robustness/brittleness issues, especially for nonlinear prob-
lems [11].
1.1. Problem setup and summary of contributions
In this work we take a different approach by employing deep neural net-
works and leverage their well known capability as universal function ap-
proximators [12]. In this setting, we can directly tackle nonlinear problems
without the need for committing to any prior assumptions, linearization, or
local time-stepping. We exploit recent developments in automatic differen-
tiation [13] – one of the most useful but perhaps underused techniques in
scientific computing – to differentiate neural networks with respect to their
input coordinates and model parameters to obtain physics informed neural
networks. Such neural networks are constrained to respect any symmetry,
invariance, or conservation principles originating from the physical laws that
govern the observed data, as modeled by general time-dependent and non-
linear partial differential equations. This simple yet powerful construction
allows us to tackle a wide range of problems in computational science and
introduces a potentially disruptive technology leading to the development
of new data-efficient and physics-informed learning machines, new classes of
numerical solvers for partial differential equations, as well as new data-driven
approaches for model inversion and systems identification.
The general aim of this work is to set the foundations for a new paradigm
in modeling and computation that enriches deep learning with the long-
standing developments in mathematical physics. These developments are
presented in the context of two main problem classes: data-driven solution
and data-driven discovery of partial differential equations. To this end, let
us consider parametrized and nonlinear partial differential equations of the
general form
u
t
+ N[u; λ] = 0,
where u(t, x) denotes the latent (hidden) solution and N[·; λ] is a nonlinear
operator parametrized by λ. This setup encapsulates a wide range of prob-
lems in mathematical physics including conservation laws, diffusion processes,
advection-diffusion-reaction systems, and kinetic equations. As a motivating
example, the one dimensional Burgers’ equation [14] corresponds to the case
3
where N[u; λ] = λ
1
uu
x
−λ
2
u
xx
and λ = (λ
1
, λ
2
). Here, the subscripts denote
partial differentiation in either time or space. Given noisy measurements of
the system, we are interested in the solution of two distinct problems. The
first problem is that of predictive inference, filtering and smoothing, or data
driven solutions of partial differential equations [9, 5] which states: given
fixed model parameters λ what can be said about the unknown hidden state
u(t, x) of the system? The second problem is that of learning, system identi-
fication, or data-driven discovery of partial differential equations [10, 6, 15]
stating: what are the parameters λ that best describe the observed data?
In this first part of our two-part treatise, we focus on computing data-
driven solutions to partial differential equations of the general form
u
t
+ N[u] = 0, x ∈ Ω, t ∈ [0, T ], (1)
where u(t, x) denotes the latent (hidden) solution, N[·] is a nonlinear differ-
ential operator, and Ω is a subset of R
D
. In what follows, we put forth two
distinct classes of algorithms, namely continuous and discrete time models,
and highlight their properties and performance through the lens of different
benchmark problems. All code and data-sets accompanying this manuscript
are available at https://github.com/maziarraissi/PINNs.
2. Continuous Time Models
We define f(t, x) to be given by the left-hand-side of equation (1); i.e.,
f := u
t
+ N[u], (2)
and proceed by approximating u(t, x) by a deep neural network. This as-
sumption along with equation (2) result in a physics informed neural net-
work f (t, x). This network can be derived by applying the chain rule for
differentiating compositions of functions using automatic differentiation [13].
2.1. Example (Burgers’ Equation)
As an example, let us consider the Burgers’ equation. This equation arises
in various areas of applied mathematics, including fluid mechanics, nonlinear
acoustics, gas dynamics, and traffic flow [14]. It is a fundamental partial
differential equation and can be derived from the Navier-Stokes equations for
the velocity field by dropping the pressure gradient term. For small values of
4
the viscosity parameters, Burgers’ equation can lead to shock formation that
is notoriously hard to resolve by classical numerical methods. In one space
dimension, the Burger’s equation along with Dirichlet boundary conditions
reads as
u
t
+ uu
x
− (0.01/π)u
xx
= 0, x ∈ [−1, 1], t ∈ [0, 1], (3)
u(0, x) = −sin(πx),
u(t, −1) = u(t, 1) = 0.
Let us define f(t, x) to be given by
f := u
t
+ uu
x
− (0.01/π)u
xx
,
and proceed by approximating u(t, x) by a deep neural network. To highlight
the simplicity in implementing this idea we have included a Python code
snippet using Tensorflow [16]; currently one of the most popular and well
documented open source libraries for machine learning computations. To
this end, u(t, x) can be simply defined as
def u(t, x):
u = neural_net(tf.concat([t,x],1), weights, biases)
return u
Correspondingly, the physics informed neural network f(t, x) takes the form
def f(t, x):
u = u(t, x)
u_t = tf.gradients(u, t)[0]
u_x = tf.gradients(u, x)[0]
u_xx = tf.gradients(u_x, x)[0]
f = u_t + u*u_x - (0.01/tf.pi)*u_xx
return f
The shared parameters between the neural networks u(t, x) and f (t, x) can
be learned by minimizing the mean squared error loss
MSE = MSE
u
+ MSE
f
, (4)
5
where
MSE
u
=
1
N
u
N
u
X
i=1
|u(t
i
u
, x
i
u
) − u
i
|
2
,
and
MSE
f
=
1
N
f
N
f
X
i=1
|f(t
i
f
, x
i
f
)|
2
.
Here, {t
i
u
, x
i
u
, u
i
}
N
u
i=1
denote the initial and boundary training data on u(t, x)
and {t
i
f
, x
i
f
}
N
f
i=1
specify the collocations points for f(t, x). The loss MSE
u
corresponds to the initial and boundary data while MSE
f
enforces the struc-
ture imposed by equation (3) at a finite set of collocation points.
In all benchmarks considered in this work, the total number of training
data N
u
is relatively small (a few hundred up to a few thousand points), and
we chose to optimize all loss functions using L-BFGS; a quasi-Newton, full-
batch gradient-based optimization algorithm [17]. For larger data-sets a more
computationally efficient mini-batch setting can be readily employed using
stochastic gradient descent and its modern variants [18, 19]. Despite the
fact that there is no theoretical guarantee that this procedure converges to
a global minimum, our empirical evidence indicates that, if the given partial
differential equation is well-posed and its solution is unique, our method is
capable of achieving good prediction accuracy given a sufficiently expressive
neural network architecture and a sufficient number of collocation points N
f
.
This general observation deeply relates to the resulting optimization land-
scape induced by the mean square error loss of equation 4, and defines an
open question for research that is in sync with recent theoretical develop-
ments in deep learning [20, 21]. Here, we will test the robustness of the
proposed methodology using a series of systematic sensitivity studies that
accompany the numerical results presented in the following.
Figure 1 summarizes our results for the data-driven solution of the Burg-
ers equation. Specifically, given a set of N
u
= 100 randomly distributed
initial and boundary data, we learn the latent solution u(t, x) by training all
3021 parameters of a 9-layer deep neural network using the mean squared
error loss of (4). Each hidden layer contained 20 neurons and a hyperbolic
tangent activation function. In general, the neural network should be given
sufficient approximation capacity in order to accommodate the anticipated
6
complexity of u(t, x). However, in this example, our choice aims to highlight
the robustness of the proposed method with respect to the well known issue
of over-fitting. Specifically, the term in MSE
f
in equation (4) acts as a reg-
ularization mechanism that penalizes solutions that do not satisfy equation
(3). Therefore, a key property of physics informed neural networks is that
they can be effectively trained using small data sets; a setting often encoun-
tered in the study of physical systems for which the cost of data acquisition
may be prohibitive.
The top panel of Figure 1 shows the predicted spatio-temporal solution
u(t, x), along with the locations of the initial and boundary training data.
We must underline that, unlike any classical numerical method for solv-
ing partial differential equations, this prediction is obtained without any
sort of discretization of the spatio-temporal domain. The exact solution for
this problem is analytically available [14], and the resulting prediction error
is measured at 6.7 · 10
−4
in the relative L
2
-norm. Note that this error is
about two orders of magnitude lower than the one reported in our previous
work on data-driven solution of partial differential equation using Gaussian
processes [9]. A more detailed assessment of the predicted solution is pre-
sented in the bottom panel of figure 1. In particular, we present a compari-
son between the exact and the predicted solutions at different time instants
t = 0.25, 0.50, 0.75. Using only a handful of initial and boundary data, the
physics informed neural network can accurately capture the intricate non-
linear behavior of the Burgers’ equation that leads to the development of a
sharp internal layer around t = 0.4. The latter is notoriously hard to ac-
curately resolve with classical numerical methods and requires a laborious
spatio-temporal discretization of equation (3).
To further analyze the performance of our method, we have performed
the following systematic studies to quantify its predictive accuracy for differ-
ent number of training and collocation points, as well as for different neural
network architectures. In table 1 we report the resulting relative L
2
error
for different number of initial and boundary training data N
u
and different
number of collocation points N
f
, while keeping the 9-layer network archi-
tecture fixed. The general trend shows increased prediction accuracy as the
total number of training data N
u
is increased, given a sufficient number of
collocation points N
f
. This observation highlights a key strength of physics
informed neural networks: by encoding the structure of the underlying phys-
7
0.0 0.2 0.4 0.6 0.8
t
−1.0
−0.5
0.0
0.5
1.0
x
u(t, x)
Data (100 points)
−0.75
−0.50
−0.25
0.00
0.25
0.50
0.75
−1 0 1
x
−1
0
1
u(t, x)
t = 0.25
−1 0 1
x
−1
0
1
u(t, x)
t = 0.50
Exact Prediction
−1 0 1
x
−1
0
1
u(t, x)
t = 0.75
Figure 1: Burgers’ equation: Top: Predicted solution u(t, x) along with the initial and
boundary training data. In addition we are using 10,000 collocation points generated using
a Latin Hypercube Sampling strategy. Bottom: Comparison of the predicted and exact
solutions corresponding to the three temporal snapshots depicted by the white vertical
lines in the top panel. The relative L
2
error for this case is 6.7 ·10
−4
. Model training took
approximately 60 seconds on a single NVIDIA Titan X GPU card.
ical law through the collocation points N
f
, one can obtain a more accurate
and data-efficient learning algorithm.
1
Finally, table 2 shows the resulting
relative L
2
for different number of hidden layers, and different number of
neurons per layer, while the total number of training and collocation points
is kept fixed to N
u
= 100 and N
f
= 10, 000, respectively. As expected, we
observe that as the number of layers and neurons is increased (hence the
capacity of the neural network to approximate more complex functions), the
1
Note that the case N
f
= 0 corresponds to a standard neural network model, i.e., a
neural network that does not take into account the underlying governing equation.
8
N
u
N
f
2000 4000 6000 7000 8000 10000
20 2.9e-01 4.4e-01 8.9e-01 1.2e+00 9.9e-02 4.2e-02
40 6.5e-02 1.1e-02 5.0e-01 9.6e-03 4.6e-01 7.5e-02
60 3.6e-01 1.2e-02 1.7e-01 5.9e-03 1.9e-03 8.2e-03
80 5.5e-03 1.0e-03 3.2e-03 7.8e-03 4.9e-02 4.5e-03
100 6.6e-02 2.7e-01 7.2e-03 6.8e-04 2.2e-03 6.7e-04
200 1.5e-01 2.3e-03 8.2e-04 8.9e-04 6.1e-04 4.9e-04
Table 1: Burgers’ equation: Relative L
2
error between the predicted and the exact solution
u(t, x) for different number of initial and boundary training data N
u
, and different number
of collocation points N
f
. Here, the network architecture is fixed to 9 layers with 20 neurons
per hidden layer.
Layers
Neurons
10 20 40
2 7.4e-02 5.3e-02 1.0e-01
4 3.0e-03 9.4e-04 6.4e-04
6 9.6e-03 1.3e-03 6.1e-04
8 2.5e-03 9.6e-04 5.6e-04
Table 2: Burgers’ equation: Relative L
2
error between the predicted and the exact solution
u(t, x) for different number of hidden layers and different number of neurons per layer.
Here, the total number of training and collocation points is fixed to N
u
= 100 and N
f
=
10, 000, respectively.
predictive accuracy is increased.
2.2. Example (Shr¨odinger Equation)
This example aims to highlight the ability of our method to handle pe-
riodic boundary conditions, complex-valued solutions, as well as different
types of nonlinearities in the governing partial differential equations. The
one-dimensional nonlinear Schr¨odinger equation is a classical field equation
that is used to study quantum mechanical systems, including nonlinear wave
propagation in optical fibers and/or waveguides, Bose-Einstein condensates,
and plasma waves. In optics, the nonlinear term arises from the intensity
dependent index of refraction of a given material. Similarly, the nonlinear
term for Bose-Einstein condensates is a result of the mean-field interactions
of an interacting, N-body system. The nonlinear Schr¨odinger equation along
9
with periodic boundary conditions is given by
ih
t
+ 0.5h
xx
+ |h|
2
h = 0, x ∈ [−5, 5], t ∈ [0, π/2], (5)
h(0, x) = 2 sech(x),
h(t, −5) = h(t, 5),
h
x
(t, −5) = h
x
(t, 5),
where h(t, x) is the complex-valued solution. Let us define f(t, x) to be given
by
f := ih
t
+ 0.5h
xx
+ |h|
2
h,
and proceed by placing a complex-valued neural network prior on h(t, x).
In fact, if u denotes the real part of h and v is the imaginary part, we
are placing a multi-out neural network prior on h(t, x) =
u(t, x) v(t, x)
.
This will result in the complex-valued (multi-output) physic informed neural
network f(t, x). The shared parameters of the neural networks h(t, x) and
f(t, x) can be learned by minimizing the mean squared error loss
MSE = MSE
0
+ MSE
b
+ MSE
f
, (6)
where
MSE
0
=
1
N
0
N
0
X
i=1
|h(0, x
i
0
) − h
i
0
|
2
,
MSE
b
=
1
N
b
N
b
X
i=1
|h
i
(t
i
b
, −5) − h
i
(t
i
b
, 5)|
2
+ |h
i
x
(t
i
b
, −5) − h
i
x
(t
i
b
, 5)|
2
,
and
MSE
f
=
1
N
f
N
f
X
i=1
|f(t
i
f
, x
i
f
)|
2
.
Here, {x
i
0
, h
i
0
}
N
0
i=1
denotes the initial data, {t
i
b
}
N
b
i=1
corresponds to the colloca-
tion points on the boundary, and {t
i
f
, x
i
f
}
N
f
i=1
represents the collocation points
on f(t, x). Consequently, MSE
0
corresponds to the loss on the initial data,
MSE
b
enforces the periodic boundary conditions, and MSE
f
penalizes the
Schr¨odinger equation not being satisfied on the collocation points.
In order to assess the accuracy of our method, we have simulated equation
10
(5) using conventional spectral methods to create a high-resolution data set.
Specifically, starting from an initial state h(0, x) = 2 sech(x) and assuming
periodic boundary conditions h(t, −5) = h(t, 5) and h
x
(t, −5) = h
x
(t, 5), we
have integrated equation (5) up to a final time t = π/2 using the Chebfun
package [22] with a spectral Fourier discretization with 256 modes and a
fourth-order explicit Runge-Kutta temporal integrator with time-step ∆t =
π/2 · 10
−6
. Under our data-driven setting, all we observe are measurements
{x
i
0
, h
i
0
}
N
0
i=1
of the latent function h(t, x) at time t = 0. In particular, the train-
ing set consists of a total of N
0
= 50 data points on h(0, x) randomly parsed
from the full high-resolution data-set, as well as N
b
= 50 randomly sampled
collocation points {t
i
b
}
N
b
i=1
for enforcing the periodic boundaries. Moreover,
we have assumed N
f
= 20, 000 randomly sampled collocation points used
to enforce equation (5) inside the solution domain. All randomly sampled
point locations were generated using a space filling Latin Hypercube Sam-
pling strategy [23].
Here our goal is to infer the entire spatio-temporal solution h(t, x) of the
Schr¨odinger equation (5). We chose to jointly represent the latent func-
tion h(t, x) = [u(t, x) v(t, x)] using a 5-layer deep neural network with
100 neurons per layer and a hyperbolic tangent activation function. Fig-
ure 2 summarizes the results of our experiment. Specifically, the top panel
of figure 2 shows the magnitude of the predicted spatio-temporal solution
|h(t, x)| =
p
u
2
(t, x) + v
2
(t, x), along with the locations of the initial and
boundary training data. The resulting prediction error is validated against
the test data for this problem, and is measured at 1.97 · 10
−3
in the rela-
tive L
2
-norm. A more detailed assessment of the predicted solution is pre-
sented in the bottom panel of Figure 2. In particular, we present a compar-
ison between the exact and the predicted solutions at different time instants
t = 0.59, 0.79, 0.98. Using only a handful of initial data, the physics informed
neural network can accurately capture the intricate nonlinear behavior of the
Schr¨odinger equation.
One potential limitation of the continuous time neural network models
considered so far, stems from the need to use a large number of collocation
points N
f
in order to enforce physics informed constraints in the entire spatio-
temporal domain. Although this poses no significant issues for problems
in one or two spatial dimensions, it may introduce a severe bottleneck in
higher dimensional problems, as the total number of collocation points needed
11
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
t
−5
0
5
x
|h(t, x)|
Data (150 points)
0.5
1.0
1.5
2.0
2.5
3.0
3.5
−5 0 5
x
0
5
|h(t, x)|
t = 0.59
−5 0 5
x
0
5
|h(t, x)|
t = 0.79
Exact Prediction
−5 0 5
x
0
5
|h(t, x)|
t = 0.98
Figure 2: Shr¨odinger equation: Top: Predicted solution |h(t, x)| along with the initial and
boundary training data. In addition we are using 20,000 collocation points generated using
a Latin Hypercube Sampling strategy. Bottom: Comparison of the predicted and exact
solutions corresponding to the three temporal snapshots depicted by the dashed vertical
lines in the top panel. The relative L
2
error for this case is 1.97 · 10
−3
.
to globally enforce a physics informed constrain (i.e., in our case a partial
differential equation) will increase exponentially. In the next section, we put
forth a different approach that circumvents the need for collocation points
by introducing a more structured neural network representation leveraging
the classical Runge-Kutta time-stepping schemes [24].
3. Discrete Time Models
Let us apply the general form of Runge-Kutta methods with q stages [24]
to equation (1) and obtain
u
n+c
i
= u
n
− ∆t
P
q
j=1
a
ij
N[u
n+c
j
], i = 1, . . . , q,
u
n+1
= u
n
− ∆t
P
q
j=1
b
j
N[u
n+c
j
].
(7)
12
Here, u
n+c
j
(x) = u(t
n
+ c
j
∆t, x) for j = 1, . . . , q. This general form en-
capsulates both implicit and explicit time-stepping schemes, depending on
the choice of the parameters {a
ij
, b
j
, c
j
}. Equations (7) can be equivalently
expressed as
u
n
= u
n
i
, i = 1, . . . , q,
u
n
= u
n
q+1
,
(8)
where
u
n
i
:= u
n+c
i
+ ∆t
P
q
j=1
a
ij
N[u
n+c
j
], i = 1, . . . , q,
u
n
q+1
:= u
n+1
+ ∆t
P
q
j=1
b
j
N[u
n+c
j
].
(9)
We proceed by placing a multi-output neural network prior on
u
n+c
1
(x), . . . , u
n+c
q
(x), u
n+1
(x)
. (10)
This prior assumption along with equations (9) result in a physics informed
neural network that takes x as an input and outputs
u
n
1
(x), . . . , u
n
q
(x), u
n
q+1
(x)
. (11)
3.1. Example (Burgers’ Equation)
To highlight the key features of the discrete time representation we revisit
the problem of data-driven solution of the Burgers’ equation. For this case,
the nonlinear operator in equation (9) is given by
N[u
n+c
j
] = u
n+c
j
u
n+c
j
x
− (0.01/π)u
n+c
j
xx
,
and the shared parameters of the neural networks (10) and (11) can be
learned by minimizing the sum of squared errors
SSE = SSE
n
+ SSE
b
, (12)
where
SSE
n
=
q+1
X
j=1
N
n
X
i=1
|u
n
j
(x
n,i
) − u
n,i
|
2
,
13
and
SSE
b
=
q
X
i=1
|u
n+c
i
(−1)|
2
+ |u
n+c
i
(1)|
2
+ |u
n+1
(−1)|
2
+ |u
n+1
(1)|
2
.
Here, {x
n,i
, u
n,i
}
N
n
i=1
corresponds to the data at time t
n
. The Runge-Kutta
scheme now allows us to infer the latent solution u(t, x) in a sequential fash-
ion. Starting from initial data {x
n,i
, u
n,i
}
N
n
i=1
at time t
n
and data at the
domain boundaries x = −1 and x = 1, we can use the aforementioned loss
function (12) to train the networks of (10), (11), and predict the solution at
time t
n+1
. A Runge-Kutta time-stepping scheme would then use this predic-
tion as initial data for the next step and proceed to train again and predict
u(t
n+2
, x), u(t
n+3
, x), etc., one step at a time.
In classical numerical analysis, these steps are usually confined to be small
due to stability constraints for explicit schemes or computational complexity
constrains for implicit formulations [24]. These constraints become more se-
vere as the total number of Runge-Kutta stages q is increased, and, for most
problems of practical interest, one needs to take thousands to millions of such
steps until the solution is resolved up to a desired final time. In sharp con-
trast to classical methods, here we can employ implicit Runge-Kutta schemes
with an arbitrarily large number of stages at effectively no extra cost.
2
This
enables us to take very large time steps while retaining stability and high pre-
dictive accuracy, therefore allowing us to resolve the entire spatio-temporal
solution in a single step.
The result of applying this process to the Burgers’ equation is presented
in figure 3. For illustration purposes, we start with a set of N
n
= 250 initial
data at t = 0.1, and employ a physics informed neural network induced by an
implicit Runge-Kutta scheme with 500 stages to predict the solution at time
t = 0.9 in a single step. The theoretical error estimates for this scheme predict
a temporal error accumulation of O(∆t
2q
) [24], which in our case translates
into an error way below machine precision, i.e., ∆t
2q
= 0.8
1000
≈ 10
−97
. To
our knowledge, this is the first time that an implicit Runge-Kutta scheme
2
To be precise, it is only the number of parameters in the last layer of the neural
network that increases linearly with the total number of stages.
14