# Cell-average based neural network method for hyperbolic and parabolic
partial differential equations

Changxin Qiu<EMAIL_ADDRESS>Jue Yan111Research work of the author is
supported by National Science Foundation grant DMS-1620335 and Simons
Foundation grant 637716<EMAIL_ADDRESS>Department of Mathematics, Iowa
State University, Ames, IA 50011, USA

###### Abstract

Motivated by finite volume scheme, a cell-average based neural network method
is proposed. The method is based on the integral or weak formulation of
partial differential equations. A simple feed forward network is forced to
learn the solution average evolution between two neighboring time steps.
Offline supervised training is carried out to obtain the optimal network
parameter set, which uniquely identifies one finite volume like neural network
method. Once well trained, the network method is implemented as a finite
volume scheme, thus is mesh dependent. Different to traditional numerical
methods, our method can be relieved from the explicit scheme CFL restriction
and can adapt to any time step size for solution evolution. For Heat equation,
first order of convergence is observed and the errors are related to the
spatial mesh size but are observed independent of the mesh size in time. The
cell-average based neural network method can sharply evolve contact
discontinuity with almost zero numerical diffusion introduced. Shock and
rarefaction waves are well captured for nonlinear hyperbolic conservation
laws.

###### keywords:

Machine learning neural network; Finite volume methods; Nonlinear hyperbolic
conservation laws.

##  1 Introduction

In this paper, we develop cell-average based neural network (CANN) method
solving time dependent hyperbolic and parabolic partial differential equations
(PDEs)

$u_{t}+f(u)_{x}=\mu u_{xx}.$ (1.1)

Our general idea is to follow available numerical schemes to build up neural
network methods. We consider to combine the powerful machine learning
mechanism of neural networks and the principles of classical numerical methods
to develop suitable neural network methods for partial differential equations
(1.1).

Machine learning with neural networks have achieved tremendous success in
image classification, text, videos and speech recognition [1, 2, 3, 4] for the
last three decades. In the last few years, connections between differential
equations and machine learning have been established, i.e. [5, 6, 7, 8, 9, 10,
11]. Very recently, machine learning neural networks have also been explored
to directly solve partial differential equations, for which a better solver
may be obtained or it can assist to improve the performance of current
numerical methods.

Nonlinear convection diffusion equation (1.1) may not be complex and hard to
solve in general. It can be used as a prototype or model equation for the more
complicated Euler and Navier-Stokes equations for fluid dynamics. There exist
quite a few numerical methods successfully developed for (1.1). Several
numerical challenges are still present. For example, a lot of ongoing efforts
are toward investigating implicit or semi-implicit methods to obtain efficient
solvers. It turns out, once well trained, our cell-average or finite volume
based neural network method can be relieved from the explicit scheme CFL
restriction and can adapt large time step size (i.e. $\Delta t=4\Delta x$) for
solution evolution, even being implemented as an explicit method.

###  1.1 Related works

Depending on the ways of neural networks applied and the goals, neural network
methods can be roughly classified into two groups. One group is to design and
apply neural network methods to directly solve partial differential equations.
The second group is to have neural networks applied to assist and improve
available numerical methods.

One popular class is to find best solution representation in terms of neural
networks. With network input vector as $x$ and $t$, such methods have the
advantage of automatic differentiation, mesh free and can be applied to solve
many types of PDEs. Approximation power of neural networks [12, 13, 14, 15,
16] are explored with these methods. For solving PDEs, we have the early works
of [17, 18], the works of [19, 20] and the popular physics-informed (PINN)
methods of [21, 22, 23, 24]. To improve the PINN efficiency on larger domain,
an extreme and distributed network is considered in [25]. Application to
incompressible Navier-Stokes equations is studied in [26]. We also have the
works of [27] and [28, 29], for which weak formulations are applied in the
loss function in stead of PDEs explicitly enforced on collocation points.
Comparison to reduced basis method is studied in [30] and performance
comparison to finite volume and discontinuous Galerkin methods are considered
in [31]. We refer to [32, 33] for PINN convergence studies.

Neural network methods have been found rather successful for solving high
dimensional PDEs. We refer to the early work of [34] and other discussions in
[35, 36, 37, 38]. Other studies include designing specific networks for
elliptic type PDEs or solving inverse problems, see [39, 40, 41]. In [42],
convolutional networks are used for estimating the mean and variance for
uncertainty quantification. We also have the works of [43] and [44], for which
method of lines approach is explored with Fourier basis considered and
Residual networks applied to evolve the dynamical system. Different to PINN or
related methods in which a global solution with neural networks is sought, our
method is similar to finite volume scheme thus is mesh dependent and is a
local solver.

The other group is to combine networks with classical numerical methods for
performance improvement. We have early result of [45] applying neural networks
as trouble-cell indicator. Deep reinforcement network in [46] is explored to
estimate the weights and enhance WENO schemes performance. In [47], network is
applied for identifying suitable amount of artificial viscosity added. Neural
network of [48] is applied to speed up iterative solver for elliptic type
PDEs. Furthermore, we have WENO schemes augmented with convolutional networks
for shock detection in [49] and DG methods with imaging edge detection
technique with convolutional network explored for strong shock detection in
[50]. We further have the work of [51] in which network is applied estimating
total variation bounded constants for DG methods.

###  1.2 Motivation and our approach

The cell-average based neural network method is closely related to finite
volume scheme. Let’s review first order upwind finite volume method for linear
advection equation

$u_{t}+u_{x}=0,$ (1.2)

which motivates the design of our cell-average based method. Integrating the
advection equation (1.2) over one computational cell
$I_{j}=[x_{j-1/2},x_{j+1/2}]$, with cell average notation
$\bar{u}_{j}=\frac{1}{\Delta x}\int_{I_{j}}u(x,t)~{}dx$ and numerical flux
$\hat{u}$ introduced at $x_{j\pm 1/2}$, we obtain
$\frac{\partial\bar{u}_{j}}{\partial
t}+\frac{\hat{u}_{j+1/2}-\hat{u}_{j-1/2}}{\Delta x}=0$. This is the starting
point of finite volume schemes. Adapting upwind for the numerical flux and
forward Euler for time discretization, we have the following first order
upwind finite volume method

$\frac{\bar{u}^{n+1}_{j}-\bar{u}^{n}_{j}}{\Delta
t}+\frac{\bar{u}^{n}_{j}-\bar{u}^{n}_{j-1}}{\Delta x}=0.$ (1.3)

For cell-average based neural network method, we consider to integrate the
advection equation (1.2) over rectangle box of $I_{j}\times(t_{n},t_{n+1})$,
which covers both the computational cell in space and the sub interval in
time. The integral or weak formulation of the advection equation is obtained
as

$\bar{u}_{j}(t_{n+1})=\bar{u}_{j}(t_{n})-\frac{1}{\Delta
x}\int^{t_{n+1}}_{t_{n}}\int_{I_{j}}u_{x}~{}dxdt.$ (1.4)

Exact solution of (1.2) satisfies the above integral formulation. With extra
regularity on the solution, the above integral format is equivalent to the
advection equation (1.2) given in differentiation format. A simple feed
forward neural network is adapted to approximate spatial variable related term
$\mathcal{N}(\overrightarrow{V}_{j}^{in};\Theta)\approx\int^{t_{n+1}}_{t_{n}}\int_{I_{j}}u_{x}~{}dxdt$.
Here $\overrightarrow{V}_{j}^{in}$ is the network input vector and $\Theta$ is
the network parameter set. And we obtain the following cell-average based
neural network method for (1.2)

$\bar{v}_{j}^{n+1}=\bar{v}_{j}^{n}+\mathcal{N}(\overrightarrow{V}_{j}^{in};\Theta^{\star}),\quad\forall
j\quad\mbox{and}\quad\forall n.$ (1.5)

Setting $\bar{v}_{j}^{n}=\bar{u}_{j}(t_{n})$ and by minimizing the difference
between network output $\bar{v}_{j}^{n+1}$ and the target
$\bar{u}_{j}(t_{n+1})$, optimal parameter set $\Theta^{\star}$ can be obtained
through offline supervised learning.

In a word, through intensive training we force the neural network to learn the
solution average evolution between two neighboring time steps. The
functionality of the parameter set $\Theta^{\star}$ is similar to the scheme
definition of upwind finite volume method (1.3), which are the coefficients of
cell-average based neural network scheme (1.5). For nonlinear convection
diffusion equation (1.1), same framework is applied that is summarized below

$\bar{u}_{j}(t_{n+1})=\bar{u}_{j}(t_{n})-\frac{1}{\Delta
x}\int^{t_{n+1}}_{t_{n}}\int_{I_{j}}\left\\{f(u)_{x}-\mu
u_{xx}\right\\}~{}dxdt\quad\approx\quad\bar{v}_{j}^{n+1}=\bar{v}_{j}^{n}+\mathcal{N}(\overrightarrow{V}_{j}^{in};\Theta^{\star}).$

Notice that we do not discretize or approximate differential terms involving
spatial variable. In stead, having the neural network handle all spatial
variable related differentiation and integration approximation is the major
idea of our method. In Figure 1 we list two simulations with our CANN method.
Shock and contact discontinuity are both sharply captured.

Figure 1: (a) & (b) Inviscid Burgers’ equation
$u_{t}+(\frac{u^{2}}{2})_{x}=0$; (c) & (d) Linear advection equation
$u_{t}+u_{x}=0$

Below we summarize the results and features of cell-average based neural
network method. Some are quite outstanding and are not common to classical
numerical methods.

  * 1.

Adapt any time step size, i.e. $\Delta t=8\Delta x$, even being an explicit
scheme

  * 2.

For Heat equation, errors are independent of time step size $\Delta t$.

  * 3.

Introduce almost zero artificial numerical diffusion for contact discontinuity
evolution

Due to some mysterious reason, the neural network method is able to catch up
solution information around the next time level $t_{n+1}$ thus allows large
time step size evolution. Different to classical numerical methods for which
we design a scheme first, here the network itself finds the _best scheme_ for
us.

The organization of the article is the following. In section §2.1, we
introduce the definition of cell-average based neural network method and
highlight its connection to finite volume schemes. In section §2.2, we
emphasize the training process and summarize the learning data difference
between linear and nonlinear PDEs. In Section §2.3, we list the major results
of cell-average based neural network method. Sequence of numerical examples
are presented in section §3. Final conclusion remarks are given in section §4.

##  2 Neural network solver

###  2.1 Problem setup, motivation and cell-averaged neural network method

We consider to develop finite volume or cell-average based neural network
(CANN) method solving partial differential equations (PDEs)

$u_{t}=\mathcal{L}(u),~{}~{}~{}(x,t)\in(a,b)\times\Re^{+}.$ (2.1)

Here $t$ and $x$ denote the time and spatial variables and $(a,b)$ is the
spatial domain. Differentiation operator $\mathcal{L}$ is introduced to
represent a generic first order hyperbolic or second order parabolic
differentiation operator. For example, we have
$\mathcal{L}(u)=(\frac{u^{2}}{2})_{x}$ for the inviscid Burgers equation and
$\mathcal{L}(u)={u}_{xx}$ for the Heat equation.

Our neural network method is mesh dependent and motivated by finite volume
method. Once well trained. the cell-average based neural network method can be
applied solving PDEs (2.1) as a regular finte volume scheme. Having a uniform
partition of $(a,b)$ into $J$ cells and $\Delta x=\frac{b-a}{J}$ is adapted as
the cell size. Denoting $x_{1/2}=a$, $x_{J+1/2}=b$, we have
$[a,b]=\bigcup^{J}_{j=1}I_{j}$ with $I_{j}=[x_{j-1/2},x_{j+1/2}]$ as one
computational cell. Furthermore, we have partition in time and adapt $\Delta
t$ for the time step size and we have $t_{n}=n\times\Delta t$ with $t_{0}=0$.

Now integrate the partial differential equation (2.1) over the computational
cell $I_{j}$ and time interval $[t_{n},t_{n+1}]$, we have

$\int^{t_{n+1}}_{t_{n}}\int_{I_{j}}~{}u_{t}~{}dxdt=\int^{t_{n+1}}_{t_{n}}\int_{I_{j}}~{}\mathcal{L}(u)~{}dxdt.$
(2.2)

With the definition of cell average $\bar{u}_{j}(t)=\frac{1}{\Delta
x}\int_{I_{j}}u(x,t)~{}dx$, equation (2.2) can be integrated out as

$\bar{u}_{j}(t_{n+1})-\bar{u}_{j}(t_{n})=\frac{1}{\Delta
x}\int^{t_{n+1}}_{t_{n}}\int_{I_{j}}~{}\mathcal{L}(u)~{}dxdt.$ (2.3)

The equation above is nothing but the integral format or a weak formulation of
(2.1). It is equivalent to the original partial differential equation (2.1)
with extra differentiation condition applied. The integral format (2.3) is the
starting point at where we design our neural network method. The idea of cell-
average based neural network method is to apply a simple fully connected
network $\mathcal{N}(\cdot;\Theta)$ to approximate the right hand side of
(2.3)

$\mathcal{N}(\cdot;\Theta)\approx\frac{1}{\Delta
x}\int^{t_{n+1}}_{t_{n}}\int_{I_{j}}~{}\mathcal{L}(u)~{}dxdt,$ (2.4)

where $\Theta$ denotes the network parameter set of all weight matrices and
biases.

Figure 2: Illustration of cell-average based neural network method

Given the solution averages $\left\\{\bar{u}^{n}_{j}\right\\}$ at time level
$t_{n}$, we apply following neural network to approximate the solution average
$\bar{u}^{n+1}_{j}$ at next time level $t_{n+1}$

$\bar{v}_{j}^{out}=\bar{v}_{j}^{in}+\mathcal{N}(\overrightarrow{V}_{j}^{in};\Theta).$
(2.5)

With $\bar{v}_{j}^{in}=\bar{u}^{n}_{j}$ and comparing (2.5) and the integral
format (2.3) of the PDEs, we have

$\bar{v}_{j}^{out}\approx\bar{u}^{n+1}_{j}.$

Vector $\overrightarrow{V}_{j}^{in}$ as the input vector of the network
$\mathcal{N}(\overrightarrow{V}_{j}^{in};\Theta)$ is an important component
that should be carefully chosen, see Figure 2. Its general format is given as

$\overrightarrow{V}_{j}^{in}=\Big{[}\bar{u}^{n}_{j-p},\cdots,\bar{u}^{n}_{j-1},\bar{u}^{n}_{j},\bar{u}^{n}_{j+1},\cdots,\bar{u}^{n}_{j+q}\Big{]}^{T},$
(2.6)

where we include the left $p$ cell averages and right $q$ cell averages of
$\bar{u}^{n}_{j}$ in the input vector. The suitable stencil or the $p$ and $q$
values in (2.6) determine the effectiveness of the neural network method
approximating the solution average $\bar{u}_{j}^{n+1}$ at the next time level.

Let’s use advection equation $u_{t}+u_{x}=0$ as an example to demonstrate the
relation between finite volume method and the cell-averaged neural network
method. Forwar Euler upwind finite volume scheme can be rewritten in terms of
vector multiplication format as

$\bar{u}^{n+1}_{j}-\bar{u}^{n}_{j}=\frac{\Delta t}{\Delta
x}\left(\bar{u}^{n}_{j}-\bar{u}^{n}_{j-1}\right)=\begin{bmatrix}\frac{\Delta
t}{\Delta x}&-\frac{\Delta t}{\Delta
x}\end{bmatrix}\begin{bmatrix}\bar{u}^{n}_{j}\\\
\bar{u}^{n}_{j-1}\end{bmatrix}.$

Now consider the cell-average based neural network method for $u_{t}+u_{x}=0$.
We adapt upwind mechanism and choose input vector as
$\overrightarrow{V}_{j}^{in}=[\bar{u}^{n}_{j},\bar{u}^{n}_{j-1}]^{T}$, see
illustration in Figure 2. The goal of the neural network is to obtain optimal
parameter set $\Theta$ through minimizing the error between network output
$\bar{v}_{j}^{out}$ and the given target $\bar{u}^{n+1}_{j}$. The parameter
set $\Theta$ is similar to the coefficients of a finite volume scheme, for
example the $\frac{\Delta t}{\Delta x}$ before $\bar{u}^{n}_{j}$ and
$-\frac{\Delta t}{\Delta x}$ before $\bar{u}^{n}_{j-1}$ of the upwind scheme.
We highlight the network parameter set $\Theta$ interact with input vector
$\overrightarrow{V}_{j}^{in}$ in a nonlinear fashion, even for a linear
differential equation. Thus neural network method behaves differently to
finite volume schemes. Once well trained, the optimal parameter set
$\Theta^{*}$ uniquely determines the _scheme definition_ of a neural network
method. We then apply this neural network solver at any cell location (index
$j$) and at any time level $t_{n}$ to evolve the solution average to next time
level $t_{n+1}$ as an explicit scheme.

Motivated by the well established finite volume schemes for hyperbolic PDEs,
we adapt upwind or method of characteristic mechanism to choose suitable
network input vector $\overrightarrow{V}_{j}^{in}$ for hyperbolic PDEs. And we
adapt a symmetric mechanism similar to central scheme for Heat equation to
choose suitable network input vector for parabolic PDEs. The guideline of
picking a suitable stencil ($p$ and $q$ in (2.6)) for the network input vector
$\overrightarrow{V}_{j}^{in}$ is summarized below.

Guideline on the choice of network input vector $\overrightarrow{V}_{j}^{in}$:

  * 1.

Hyperbolic PDEs: include the characteristic or domain of dependence

  * 2.

Parabolic PDEs: a symmetric stencil mechanism

In this paper, we consider a standard fully connected neural network with $M$
($M\geq 3$) layers. The input and output vectors of the network are the first
and last layer. Among the total $M$ layers, the interior $(M-2)$ are the
hidden layers. Thus the minimum structure of the neural network involves one
hidden layer with $M=3$. We have $n_{i}$ ($i=1,\cdots,M$) denote the number of
neurons in each layer. The first layer is the input vector with $n_{1}=p+q+1$
as its dimension and the last layer is the output vector with $n_{M}=1$ as its
dimension. The abstract goal of machine learning is to find a function
$\mathcal{N}:R^{p+q+1}\to R^{1}$ such that $\mathcal{N}(\cdot;\Theta)$
accurately approximates $\frac{1}{\Delta
x}\int^{t_{n+1}}_{t_{n}}\int_{I_{j}}~{}\mathcal{L}(u)~{}dxdt$, the right hand
side of (2.3). The optimal parameter set $\Theta$ of the network
$\mathcal{N}(\cdot;\Theta)$ will be obtained by training the network
intensively over the given data set.

Learning data are collected in the form of pairs. Each pair refers to the
solution averages at two neighboring time level. The training data set is
denoted as

$S=\left\\{\left(\bar{u}_{j}^{n},\bar{u}_{j}^{n+1}\right),\,\,j=1,\cdots,J\right\\}_{n=0}^{m},$
(2.7)

which are solution averages obtained from highly accurate numerical method for
the PDE. We highlight the training data pairs are solution averages collected
over the spatial domain and from time levels $t_{0}$ to $t_{m+1}$.

The fully neural network $\mathcal{N}(\cdot;\Theta)$ consists of $M$ layers.
Each two consecutive layers is connected with an affine linear transformation
and a point-wise nonlinear activation function. The function or the mapping
$\mathcal{N}(\cdot;\Theta)$ is a composition of following operators,

$\mathcal{N}(\cdot;\Theta)=(\sigma_{M}\circ\textbf{W}_{M-1})\circ\cdots\circ(\sigma_{2}\circ\textbf{W}_{1}).$
(2.8)

where $\circ$ stands for operator composition. We have $\bf{W}_{i}$ denoting
the linear transformation operator or the weight matrix connecting the neurons
from $i$-th layer to $(i+1)$-th layer. The parameter set is further augmented
with the biases vectors. We have $\sigma_{i}:R\to R$ as the activation
function ($i\geq 2$), which is applied to each neuron of the $i$-th layer in a
component-wise fashion. In this paper we apply $tanh(x)$ function as the
activation function. Specifically $\sigma_{i}=tanh(x)$ is applied between all
layers, except to the output layer for which we have $\sigma_{M}(x)=x$.

###### Definition 1

A cell-average based neural network method is uniquely determined by the
following four components: (1) the choice of spatial mesh size $\Delta x$; (2)
the choice of time step size $\Delta t$; (3) the choice of network input
vector $\overrightarrow{V}_{j}^{in}$ of (2.6); and (4) the number of hidden
layers and neurons per layer the corresponding structure of the neural
network.

We further highlight that our cell-average based neural network method, once
well trained, will be implemented as a regular explicit finite volume scheme.
Notice the cell-average based neural network method is designed to approximate
the weak or integral format of the partial differential equations (2.3), not
the original partial differential equations (2.1) given in differentiation
format.

###  2.2 Training Process

In this section, we discuss how to train the network to obtain optimal
parameter set $\Theta^{*}$ such that the neural network (2.5) can accurately
approximate the solution average evolution
$\bar{u}^{n}_{j}\to\bar{u}^{n+1}_{j}$. This is achieved by applying
$\bar{v}^{in}_{j}=\bar{u}_{j}^{n}$ in (2.5) to obtain network output
$\bar{v}_{j}^{out}$, comparing with the target $\bar{u}_{j}^{n+1}$, and then
looping among the data set $S$ of (2.7) to minimize the error or the squared
loss function

$L_{j,t_{n}}(\Theta)=(\bar{v}_{j}^{out}-\bar{u}_{j}^{n+1})^{2},$ (2.9)

for all $j=1,\cdots,J$ and for all $n=0,\cdots,m$. This choice of loss
function defined over one single data pair is referred as the stochastic or
approximate gradient descent method. Notice for $j=1$ or $j=J$ or those close
to boundary cells, the stencil or network input vector
$\overrightarrow{V}_{j}^{in}$ of (2.6) requires averages values of $p$ cells
to the left and $q$ cells to the right of the current cell. In this paper we
only consider Dirichlet or periodic boundary conditions. Thus we simply copy
solution averages from inside the domain implementing periodic boundary
condition. And we assign exact solution values to those out of domain ghost
cells solving Burgers equation Riemann problems.

Recall we have $S$ of (2.7) denoting the training data set, which include
solution averages spread over the spatial domain with index $j$ and up to time
levels $t_{m+1}$. Numerical tests show that for linear partial differential
equations, training data of $(t_{0},t_{1})$ with $m=0$ in $S$ is sufficient
for obtaining optimal parameter set $\Theta$. For nonlinear partial
differential equations, i.e. the Burgers equation, multiple time levels ($m>0$
in $S$) data are necessary to guarantee the neural network effectively
learning the solution evolution mechanism, for example capturing the right
shock speed.

Figure 3: Network parameter one epoch training for cell-average based neural
network (CANN) method

As discussed in section 2.1, the network parameter set $\Theta$ behave more or
less as the coefficients of a finite volume scheme. Parameter set $\Theta$ is
independent of at what cell location and at what time level the neural network
method is applied. Thus we should minimize the loss function over all spatial
index $j$ and at all time levels in $S$ to obtain the optimal parameter set
$\Theta^{*}$ or the _coefficients_ of the neural network method.

One epoch of iteration is defined as the following. The weights and biases are
first randomly generated from normal distribution around zero. The parameter
set $\Theta$ are updated sequentially through space and time with all data
pairs applied _once_ in the training set $S$. In a word, we start with first
time level pair $(t_{0},t_{1})$ and go through all spatial cells and
iteratively update the parameter set $\Theta^{j,t_{0}}$ through the gradient
descent direction of (2.9). Then we move on to the next time level pairs
$(t_{n},t_{n+1})$ ($n\geq 1$) and go through all spatial cells again updating
the parameter set $\Theta^{j,t_{n}}$. Below is an illustration of one epoch
definition, for which we use solution averages pairs in the training set $S$
once iterating the network parameter, see Figure 3.

$\Theta^{j,t_{n}}_{old}\longrightarrow\Theta^{j,t_{n}}_{new},\quad
j=1,\cdots,J,\quad n=0,\cdots,m.$

The superscripts $j$ and $t_{n}$ notations are added to specify the iteration
or update is processed at what cell and at what time level. They all refer to
the update of the same optimized parameter set $\Theta^{*}$. Learning rate
$\alpha=0.01$ or $\alpha=0.001$ is applied, when updating through stochastic
gradient descent method.

In this paper, we have integer $K$ introduced denoting the total number of
epochs iterations involved in training. Instead of iterating through one epoch
defined above, at each time level we run total $K$ iterations to update the
parameter set $\Theta$ with one iteration defined as the one round of updating
over spatial domain (index $j$). We then continue the iteration of $\Theta$ to
the next time level, which can be summarized as

$\Theta^{j,t_{n},i}\longrightarrow\Theta^{j,t_{n},i+1},\quad
j=1,\cdots,J,\quad i=1,\cdots,K,\quad n=0,\cdots,m.$ (2.10)

In the end, the parameter set $\Theta$ have been iteratively updated with
total $J\times K\times(m+1)$ times. Specifically for the last time level pair
$(t_{m},t_{m+1})$, we record the squared $L_{2}$ error defined below
corresponding to iteration

$L^{2}_{2}(t_{m+1})=\sum^{J}_{j=1}\Delta
xL_{j,t_{m}}(\Theta)=\sum^{J}_{j=1}(\bar{v}_{j}^{out}-\bar{u}_{j}^{m+1})^{2}\Delta
x.$ (2.11)

We output the squared $L_{2}$ error of (2.11) corresponding to iteration index
$i=1,\cdots,K$ to demonstrate the effectiveness of cell-average based neural
network method. Now we conclude the section with comment on the minimum size
of training data set $S$ of (2.7), which are purely lab results observed from
numerical tests.

###### Remark 1

For linear partial differential equations, one time level solution averages in
the training set $S$ corresponding to $(t_{0},t_{1})$ or $m=0$ in (2.7), is
sufficient for obtaining an effective neural network. For nonlinear partial
differential equations, it is necessary to include multiple time levels
($m>0$) of solution averages in the training set $S$ to have the neural
network learn the evolution mechanism successfully.

For example in the numerical section when we solve inviscid Burgers’ equation
with smooth $\sin(x)$ initial in Example 3.3.1, we have solution averages data
pairs up to $t=2$ included in the training set $S$. This is because a shock
starts to develop at $t=1$. We need to have more time levels solution averages
included in the training set to make sure the neural network solver be able to
learn the shock capturing mechanism.

###  2.3 Implementation and summary of cell-average based neural network
method

With the optimal weights and biases $\Theta^{*}$ obtained and the neural
network $\mathcal{N}(\overrightarrow{V}_{j}^{in};\Theta^{*})$ well defined and
available, the cell-average based neural network method can be implemented as
a regular _explicit finite volume scheme_ as below

$\bar{v}_{j}^{n+1}=\bar{v}_{j}^{n}+\mathcal{N}(\overrightarrow{V}_{j}^{n};\Theta^{*}),\quad\forall
j=1,\cdots,J,\quad\forall n=0,1,2,\cdots.$ (2.12)

Again, with the spatial and time step sizes of $\Delta x$ and $\Delta t$ and
the previously chosen network input vector of
$\overrightarrow{V}_{j}^{n}=\Big{[}\bar{v}^{n}_{j-p},\cdots,\bar{v}^{n}_{j-1},\bar{v}^{n}_{j},\bar{v}^{n}_{j+1},\cdots,\bar{v}^{n}_{j+q}\Big{]}^{T}$,
together with the network optimal parameter set $\Theta^{*}$, we have a
complete definition of a neural network method.

For cell-average based neural network method (2.12), we assign ghost cell
values and apply boundary conditions as a regular finite volume method. Again,
we only consider Dirichlet or periodic boundary conditions in this paper. A
generic stencil include $p$ cells to the left and $q$ cells to the right to
evolve current cell to the next time level. Similar to the training process
discussed in section 2.2, we simply copy solution averages from inside the
domain for periodic boundary conditions. And we assign exact solution values
to those out of domain ghost cells for Burgers’ equation Riemann problems,
i.e. for Dirichlet boundary conditions.

One amazing result is that cell-average based neural network method can be
_relieved_ from the CFL restriction on time step size, especially for
parabolic PDEs. With explicit discretization in time, classical numerical
method requires time step size to be as small as $\Delta t\approx(\Delta
x)^{2}$, which is very expensive for multi-dimensional problems. It turns out
neural network method can adapt to any time step size $\Delta t$ that is
independent of spatial size $\Delta x$, and still a stable method can be
obtained. Once well trained, the neural network method can efficiently and
accurately evolve the solution forward in time as an explicit method.

Recall our neural network method is based on the following integral format of
the partial differential equation $u_{t}+\mathcal{L}(u)=0$

$\bar{u}_{j}(t_{n+1})-\bar{u}_{j}(t_{n})=\frac{1}{\Delta
x}\int^{t_{n+1}}_{t_{n}}\int_{I_{j}}~{}\mathcal{L}(u)~{}dxdt,$

with $\mathcal{L}$ as the differentiation operator on the spatial variable. It
is equivalent to the PDEs given in differentiation format. Due to some
mysterious reason, the neural network method is able to catch up solution
information around time level $t_{n+1}$ thus allows large time step size
evolution as an implicit method. Below we summarize the major result of our
CANN method.

###### Lemma 1

Even being an explicit scheme, cell-average based neural network method (2.12)
can adapt to any time step size $\Delta t$ for hyperbolic and parabolic PDEs
(2.1). With $\Delta t=\Delta x$, first order of accuracy is obtained with
neural network method (2.12) for linear advection and linear convection
diffusion equations. For Heat equation, neural network method errors depend
only on spatial mesh size $\Delta x$ and are independent of time step size
$\Delta t$.

Besides accuracy and error behavior studies, we also carry our a series of
numerical tests with cell-average based neural network method. Now we list the
advantages of our CANN method that are not common for classical numerical
methods.

  * 1.

Allow large time step size evolution, i.e. $\Delta t=8\Delta x$ for linear
hyperbolic PDEs.

  * 2.

For Heat equation, errors are independent of time step size. Similar errors
are obtained for $\Delta t=\Delta x$, $\Delta t=2\Delta x$ and $\Delta
t=4\Delta x$ when same spatial mesh $\Delta x$ is adapted.

  * 3.

Introduce almost zero artificial numerical diffusion for contact discontinuity
propagation.

  * 4.

Introduce little numerical dissipation and dispersive errors after long time
run.

###### Remark 2

Once one cell-average based neural network is well trained and available, it
can be applied to solve same PDE associated with different initials and over
different domains.

###### Remark 3

It remains unknown how to choose a suitable neural network architecture in
terms of number of hidden layers and neurons per layer. Numerical tests show
one or two hidden layers with a few neurons work well. We further mention that
some network structure may lead to extremely small errors and can be even
regarded as a perfect solver, for which regular order of convergence over
refined mesh error analysis does not hold anymore.

##  3 Numerical Example

In this section, we carry out a series of numerical tests to check out the
accuracy and capability of cell-average based neural network methods. We start
with linear advection equation, Heat equation and convection diffusion
equations to investigate if or not order of convergence can be observed. Then
we move on to nonlinear hyperbolic conservation law to test the capability of
neural network method capturing shock and rarefaction waves propagation.

Over the section we adapt $T$ as a generic final time at where we compute the
errors and orders. As mentioned previously, time step size $\Delta t$ is
chosen before and is a part of the definition of a neural network method
(2.12). Thus we always have $T$ as an integer multiple times of $\Delta t$.
Below we list the $L_{2}$ and $L_{\infty}$ errors formula as used in finite
volume methods.

$Error_{L_{2}}(T)=\sqrt{\sum^{J}_{j=1}(\bar{v}_{j}(T)-\bar{u}_{j}(T))^{2}\Delta
x}$ (3.1)
$Error_{L_{\infty}}(T)=\max^{J}_{j=1}|\bar{v}_{j}(T)-\bar{u}_{j}(T)|$ (3.2)

Again, we have $\bar{v}_{j}(T)$ denote our CANN method (2.12) solution on cell
$j$ and at final time $T$. We have $\bar{u}_{j}(T)$ denote either the exact
solution average or the reference solution average on cell $j$ and at time $T$
that is obtained from a highly accurate numerical method.

As marked down in Definition 1, four components of $\Delta x$, $\Delta t$,
network input vector $\overrightarrow{V}_{j}^{in}$ and the structure of neural
network (number of hidden layers and neurons per layer) together identify one
neural network solver (2.12). Most of the time we follow the principle listed
in section 2.1 to choose suitable stencil width or the network input vector of
$\overrightarrow{V}_{j}^{in}$. We highlight that for linear PDEs one time
level $(t_{0},t_{1})$ solution average data pair is applied training the
network and for nonlinear PDEs multiple time levels solution average data pair
are needed to obtain an effective neural network solver, see Remark 1. We also
mention that squared $L_{2}$ error of (2.11) around $10^{-8}$ or smaller is
used as stop condition in training.

###  3.1 Linear advection equation

In this subsection, we focus on linear hyperbolic equation of

$u_{t}+u_{x}=0.$ (3.3)

Even the above linear equation is a simple model, quite a few problems can be
tested to evaluate the capability of a numerical method. We consider three
problems for advection (3.3). One is about smooth function evolution and we
check if order of convergence can be observed. Then we study the contact
discontinuity propagation problem. It is not a trivial test, since numerical
methods tend to either generate smeared out approximations or numerical
oscillations around the discontinuity. For the third test, we have the neural
network method simulate wave propagation after long time run. Neural network
method is quite stable and produces little dispersive and dissipation errors
after long time simulation.

###### Example 3.1.1

_smooth wave propagation_

In this example we solve (3.3) with initial condition $u(x,0)=\sin(x)$.
Spatial domain is taken as $D=[0,2\pi]$ and the exact solution
$u(x,t)=\sin(x-t)$ is a smooth wave propagating from left to right. We
consider a simple neural network with input vector of

$\overrightarrow{V}_{j}^{in}=\Big{[}\bar{u}^{n}_{j-1},\bar{u}^{n}_{j}\Big{]}^{T},$
(3.4)

that is similar to the upwind finite volume scheme. The network picked
consists of 2 hidden layers with 5 neurons per layer. Total iterations of
$K=5\times 10^{5}$ is applied to train the networks. Periodic boundary
conditions are applied. We carry out two accuracy tests, for which we compute
the $L_{2}$ and $L_{\infty}$ errors of (3.1) and (3.2) at final time $T=\pi$.

Case I:

For the first case study we mimic the accuracy check of a standard numerical
method. We consider four spatial meshes of $\Delta x=\pi/10$, $\pi/20,\pi/40,$
and $\pi/80$. For each cell size $\Delta x$, time step size $\Delta t=\Delta
x$ is taken correspondingly. The four well trained neural network solvers are
very similar to each other except the mesh size. All four network solvers are
used to solve the smooth wave propagation to final time $T=\pi$ and $L_{2}$
and $L_{\infty}$ errors are computed and listed in Table 1.

$\Delta x$ | $L_{2}$ | order | $L_{\infty}$ | order
---|---|---|---|---
$\pi/10$ | $1.8756e^{-2}$ |  | $1.0237e^{-2}$ |
$\pi/20$ | $8.0830e^{-3}$ | 1.21 | $4.7403e^{-3}$ | 1.11
$\pi/40$ | $1.5547e^{-3}$ | 2.39 | $9.7037e^{-4}$ | 2.29
$\pi/80$ | $6.3500e^{-4}$ | 1.28 | $3.9838e^{-4}$ | 1.28
Table 1: Errors and orders of neural network methods for Example 3.1.1 (Case
I), $\Delta t=\Delta x$

Case II:

For this case we keep the spatial mesh size $\Delta x=\pi/40$ fixed, but
consider three different time step sizes of $\Delta t=2\Delta x$, $\Delta
t=5\Delta x$ and $\Delta t=8\Delta x$. The three neural network solvers apply
same upwind like network input vector of (3.5). The only difference is the
time step size. We list the $L_{2}$ and $L_{\infty}$ errors computed at final
time $T=\pi$ in Table 2. Three network solvers give similar errors and it is
not clear if the network errors relate to $\Delta t$. For all three time step
sizes, i.e. $\Delta t=8\Delta x$, the principle of having domain of dependence
included is not followed. We simply apply (3.5) as the input vector that is
like the upwind scheme. These settings conflict with method of characteristic.
But all neural network solvers work well and give errors similar to those in
case I.

$\Delta t$ | $L_{2}$ | $L_{\infty}$
---|---|---
2$\Delta x$ | $7.0431e^{-3}$ | $4.1544e^{-3}$
5$\Delta x$ | $9.2344e^{-3}$ | $5.3628e^{-3}$
8$\Delta x$ | $6.9895e^{-3}$ | $3.1796e^{-3}$
Table 2: Errors and orders of neural network methods for Example 3.1.1 (Case
II), $\Delta x=\pi/40$ fixed with varying $\Delta t$

###### Example 3.1.2

_contact discontinuity_

In this example, we solve advection equation (3.3) with initial condition

$u(x,0)=\left\\{\begin{aligned} 1~{}~{}x\leq 0,\\\
2~{}~{}x>0,\end{aligned}\right.$

over domain $D=[-1,4]$. Periodic boundary condition is applied. The contact
discontinuity is initially located at $x=0$, which moves back into the domain
after $t>4$. Same network input vector of (3.5) that is similar to upwind
scheme is considered. One hidden layer of 10 neurons is the chosen network
structure and a total of $K=5\times 10^{5}$ iterations are applied training
the network. We have $\Delta x=\Delta t=\frac{1}{20}$. Snapshots of our CANN
method simulation are presented in Figure 4. The contact discontinuity is
sharply resolved, even after one period of evolution at $t=5$. No oscillation
is observed and there is almost no artificial diffusion introduced with the
neural network method. Numerical simulation with CANN method behaves better
than many numerical methods that tend to generate smeared out simulations.

Figure 4: contact discontinuity evolution (Example 3.1.2) with neural network
method

###### Example 3.1.3

_long time evolution_

In this example we solve advection equation (3.3) with initial condition

$u(x,0)=\sin(x),$

over domain $D=[0,2\pi]$. We choose $\Delta x=\frac{2\pi}{100}$ and $\Delta
t=4\Delta x$. Long time simulation with CANN method is considered. To include
the domain of dependence, following neural network input vector is taken

$\overrightarrow{V}_{j}^{in}=\Big{[}\bar{u}^{n}_{j-6},\bar{u}^{n}_{j-5},\bar{u}^{n}_{j-4},\bar{u}^{n}_{j-3},\bar{u}^{n}_{j-2},\bar{u}^{n}_{j-1},\bar{u}^{n}_{j}\Big{]}^{T}.$
(3.5)

The chosen network is with 1 hidden layer and 10 neurons. A total of
$K=5\times 10^{5}$ iterations are applied. We compute $L_{2}$ and $L_{\infty}$
errors at different time locations, which are listed in Table 3. After several
periods, neural network method is still able to accurately capture the
solution evolution and gives small dispersive and dissipation errors.

| $t=4\pi/5$ | $t=2\pi$ | $t=4\pi$ | $t=8\pi$
---|---|---|---|---
$L_{2}$ | $9.1869e^{-14}$ | $1.0045e^{-12}$ | $1.8625e^{-10}$ | $7.1709e^{-6}$
$L_{\infty}$ | $1.2620e^{-13}$ | $1.0600e^{-12}$ | $1.7158e^{-10}$ | $5.3325e^{-6}$
Table 3: $L_{2}$ and $L_{\infty}$ errors (Example 3.1.3) after long time
simulation with neural network method

###  3.2 Linear convection diffusion equation

In this section, we solve linear convection diffusion equation with cell-
average based neural network method. It is further confirmed that we can be
relieved from explicit scheme small time step size restriction. For Heat
equation, numerical tests show the errors only relate to spatial mesh size
$\Delta x$ and is independent of time step size $\Delta t$.

###### Example 3.2.1

_Heat equation_

In this example, we consider solving the Heat equation

$u_{t}=u_{xx},$

with initial $u(x,0)=\sin(\pi x)$ and over domain $D=[0,1]$. Exact solution is
$u(x,t)=e^{-\pi^{2}t}\sin(\pi x)$. Different to hyperbolic PDEs, we have
infinite speed of propagation for Heat equation. Here we follow a similar to
central scheme symmetric mechanism to pick up the network input vector as

$\overrightarrow{V}_{j}^{in}=\Big{[}\bar{u}^{n}_{j-3},\bar{u}^{n}_{j-2},\bar{u}^{n}_{j-1},\bar{u}^{n}_{j},\bar{u}^{n}_{j+1},\bar{u}^{n}_{j+2},\bar{u}^{n}_{j+3}\Big{]}^{T}.$
(3.6)

The chosen network is with 2 hidden layers and 15 neurons per layer.
Iterations are run up to $K=10^{5}$ times to optimize the network parameter
set. Periodic boundary conditions are applied. Final time $T=0.1$ is adapted
to compute $L_{2}$ and $L_{\infty}$ errors of (3.1) and (3.2). We carry out
three accuracy tests. For the first and third tests, we always choose $\Delta
t=\Delta x$ and have $\Delta x$ refined to check out order of convergence. In
the second test we have $\Delta x$ fixed and choose $\Delta t$ as a multiple
times of $\Delta x$. Same network structure is applied, but different choices
of $\Delta t$, $\Delta x$ or network input vector
$\overrightarrow{V}_{j}^{in}$ are considered in each test.

Case I:

For this case, we investigate the convergence order of neural network method.
We consider four settings of $\Delta
x=\frac{1}{40},\frac{1}{80},\frac{1}{160},\frac{1}{320}$. For each cell size
$\Delta x$, we choose $\Delta t=\Delta x$ correspondingly. The four network
solvers are similar to each other except the mesh size. We solve the Heat
equation with each neural network solver to $T=0.1$ and compute the $L_{2}$
and $L_{\infty}$ errors. In Table 4 we list all errors and orders. Clean first
order of convergence is observed with CANN method.

$\Delta x$ | $L_{2}$ | order | $L_{\infty}$ | order
---|---|---|---|---
$1/40$ | $8.6949e^{-3}$ |  | $2.0873e^{-2}$ |
$1/80$ | $4.5270e^{-3}$ | 0.94 | $1.4104e^{-2}$ | 0.57
$1/160$ | $2.4736e^{-3}$ | 0.87 | $7.2650e^{-3}$ | 0.96
$1/320$ | $1.2894e^{-3}$ | 0.94 | $3.7860e^{-3}$ | 0.94
Table 4: Errors and orders of neural network methods for Heat equation (Case
I), $\Delta t=\Delta x$

Case II:

For this case we fix the spatial mesh size $\Delta x=1/160$ but apply three
time step sizes of $\Delta t=\Delta x$, $\Delta t=2\Delta x$ and $\Delta
t=4\Delta x$. We compute the $L_{2}$ and $L_{\infty}$ errors with each neural
network solver at final time $T=0.1$. Errors are listed in Table 5. The three
network solvers give similar errors. The accuracy of CANN method seems to be
independent of time step size $\Delta t$.

$\Delta t$ | $L_{2}$ | $L_{\infty}$
---|---|---
4$\Delta x$ | $2.1981e^{-3}$ | $6.8272e^{-3}$
2$\Delta x$ | $2.4969e^{-3}$ | $7.2399e^{-3}$
$\Delta x$ | $2.4736e^{-3}$ | $7.2650e^{-3}$
Table 5: Errors of neural network methods for Heat equation (Case II), $\Delta
x=1/160$ fixed, different $\Delta t$

Case III:

Motivated by infinite speed of propagation, we wonder if the error and
accuracy may be improved with increased stencil width of the network input
vector, the $p$ and $q$ values of (2.6). Three spatial mesh sizes of $\Delta
x=1/40,1/80,1/160$ with $\Delta t=\Delta x$ are studied. Different to case I,
we gradually increase the input vector stencil width with refined mesh. We
have $p=q=2$ for $\Delta x=1/40$, $p=q=4$ for $\Delta x=1/80$ and $p=q=8$ for
$\Delta x=1/160$. The three networks use same spatial cells to evolve the
solution average. Errors are listed in Table 6. There is no sign of
improvement with wider stencil included.

$\Delta x$ | $L_{2}$ | $L_{\infty}$
---|---|---
1/40 | $7.1179e^{-3}$ | $1.8046e^{-2}$
1/80 | $6.3502e^{-3}$ | $1.8319e^{-2}$
1/160 | $6.4138e^{-3}$ | $2.1510e^{-2}$
Table 6: Errors of neural network methods for Heat equation (Case III),
$\Delta t=\Delta x$, wider stencil on refined mesh

###### Example 3.2.2

_Linear convection diffusion equation_

In this subsection, we consider linear convection-diffusion equation

$u_{t}=u_{xx}+u_{x},~{}~{}~{}x\in D,~{}~{}t\geq 0,$

with initial $u(x,0)=\sin(x)$. Spatial domain is $D=[0,2\pi]$. We further
check whether the first order of accuracy can be obtained. We choose same
network input vector of (3.6) as for Heat equation. Periodic boundary
conditions are considered. Exact solution is available with
$u(x,t)=e^{-t}\sin(x+t)$. The picked network structure involves 1 hidden layer
and 15 neurons. Number of iterations is taken as $K=5\times 10^{6}$. We adapt
final time $T=\pi/4$ to compute all errors and orders.

For the order of convergence test, four spatial mesh sizes of $\Delta
x=\frac{\pi}{40},\frac{\pi}{80},\frac{\pi}{160},\frac{\pi}{320}$ are studied.
For each $\Delta x$, we choose $\Delta t=\Delta x$ correspondingly. Again, the
four neural network solvers are similar to each other except the mesh size. In
Table 7, we list the $L_{2}$ and $L_{\infty}$ errors. Roughly first order of
convergence is observed.

$\Delta x$ | $L_{2}$ | order | $L_{\infty}$ | order
---|---|---|---|---
$\pi/40$ | $5.3013e^{-3}$ |  | $3.2397e^{-3}$ |
$\pi/80$ | $1.3801e^{-3}$ | 1.94 | $7.5132e^{-4}$ | 2.11
$\pi/160$ | $5.0091e^{-4}$ | 1.46 | $2.6084e^{-4}$ | 1.52
$\pi/320$ | $2.6771e^{-4}$ | 0.91 | $1.5045e^{-4}$ | 0.89
Table 7: Errors and orders of neural network methods for linear convection
diffusion equation, $\Delta t=\Delta x$.

The second test is similar to the Case II of Example 3.2.1. We fix the spatial
size $\Delta x=\pi/160$ and vary the time step size from $\Delta t=\Delta x$,
$\Delta t=2\Delta x$ to $\Delta t=4\Delta x$. The computed $L_{2}$ and
$L_{\infty}$ errors of the three network solvers are listed in Table 8. The
errors seem to be related to time step size $\Delta t$, even the relationship
is not clear. Again, the cell-average based neural network allows large time
step. Notice the choice of $\Delta t=4\Delta x$ involves a time step size
roughly $200$ times bigger than the regular CFL restriction of $\Delta
t\approx\Delta x^{2}$.

$\Delta t$ | $L_{2}$ | $L_{\infty}$
---|---|---
4$\Delta x$ | $1.0249e^{-3}$ | $5.8832e^{-4}$
2$\Delta x$ | $8.2225e^{-4}$ | $4.2640e^{-4}$
$\Delta x$ | $5.0091e^{-4}$ | $2.6084e^{-4}$
$\Delta x/2$ | $3.7288e^{-4}$ | $2.0151e^{-4}$
Table 8: Errors of neural network methods for linear convection diffusion
equation, $\Delta x=\pi/160$ fixed, different $\Delta t$

###  3.3 Nonlinear hyperbolic equation

In this section, we investigate the effectiveness of neural network method
solving nonlinear conservation laws. We consider the inviscid Burgers’
equation

$u_{t}+\left(\frac{u^{2}}{2}\right)_{x}=0,~{}~{}(x,t)\in D\times R^{+}.$ (3.7)

Four benchmark problems are studied to illustrate the computational challenges
of nonlinear hyperbolic PDEs. We first consider the smooth initial $\sin(x)$
wave which develops into a shock after finite time evolution. Then we study
two Riemann problems with two piece wise constants initials that will develop
into a shock and a rarefaction wave. Last example involves three piece wise
constants initial that will develop into the interaction between rarefaction
wave and shock. In all four examples, Dirichlet boundary conditions are
applied.

For nonlinear problems, multiple time levels of data pairs are necessary and
applied training the neural network. For all four examples, time step size
$\Delta t=0.1$ is taken. We roughly have $\Delta t\approx 2\Delta x$. Solution
cell average values up to $t=2.0$, as listed below

$S=\left\\{\left(\bar{u}_{j}^{n},\bar{u}_{j}^{n+1}\right),\,\,j=1,\cdots,J\right\\}_{n=0}^{m=19},$
(3.8)

are included in the training data set. So we have twenty time levels of
solution average pairs $\left(\bar{u}_{j}^{n},\bar{u}_{j}^{n+1}\right)$ used
in training the network. Spatial domain is either $D=[0,2\pi]$ or $D=[-1,5]$.
Mesh size of $\Delta x=\frac{2\pi}{100}$ for example 3.3.1 and $\Delta
x=\frac{6}{100}$ for other three examples are applied.

(a) Example 3.3.1: sine initial (b) Example 3.3.2: shock wave (c) Example
3.3.3: rarefaction wave (d) Example 3.3.4: rarefaction and shock interaction

Figure 5: Invisid Burgers’ equation neural network squared $L_{2}$ training
errors (2.11) at the last time level

###### Example 3.3.1

_Smooth initial sine wave_

We first consider the inviscid Burgers equation (3.7) associated with smooth
initial value of

$u(x,0)=\sin(x),$

over domain $D=[0,2\pi]$. Zero boundary conditions $u(0,t)=u(2\pi,t)=0$ and
its zero extension to out of domain ghost cells are applied. With $\Delta
t=0.1$ and $\Delta x=2\pi/100$ and characteristic speed less than one, the
network input vector is taken as

$\overrightarrow{V}_{j}^{in}=\Big{[}\bar{u}^{n}_{j-3},\bar{u}^{n}_{j-2},\bar{u}^{n}_{j-1},\bar{u}^{n}_{j},\bar{u}^{n}_{j+1},\bar{u}^{n}_{j+2},\bar{u}^{n}_{j+3}\Big{]}^{T}.$

Figure 6: neural network simulation for smooth $\sin(x)$ wave initial (Example
3.3.1)

This stencil includes the characteristic and the domain of dependence.
Training data pairs $\left(\bar{u}_{j}^{n},\bar{u}_{j}^{n+1}\right)$ are
obtained from highly accurate discontinuous Galerkin method. The network
itself contains 2 hidden layers with 8 neurons per layer. The network training
is conducted for up to $K=10^{5}$ iterations. We also output the squared
$L_{2}$ training error of (2.11) corresponding to iteration at the very last
time level $t=2$, see Fig 5.

Well trained network is applied to simulate solution evolution to $T=3$.
Before and after shock developed are illustrated in Fig 6. The cell-average
based neural network results are comparable to those obtained with classical
numerical methods. Shock evolution is sharply captured and no oscillation is
generated with the neural network method.

###### Example 3.3.2

_Single shock propagation_

In this example, we consider a Riemann problem of (3.7) with piece wise
constant initial

$u(x,0)=\left\\{\begin{array}[]{lr}1,~{}~{}~{}~{}x<0,\\\
0,~{}~{}~{}otherwise.\end{array}\right.$

Computational domain is $D=[-1,5]$. Dirichlet boundary conditions of
$u(-1,t)=1$, $u(5,t)=0$ and its out of domain extensions are applied. Cell
size $\Delta x=0.06$ and time step size $\Delta t=0.1$ are taken. The network
input vector of (2.6) is chosen as

$\overrightarrow{V}_{j}^{in}=\Big{[}\bar{u}^{n}_{j-4},\bar{u}^{n}_{j-3},\bar{u}^{n}_{j-2},\bar{u}^{n}_{j-1},\bar{u}^{n}_{j},\bar{u}^{n}_{j+1},\bar{u}^{n}_{j+2}\Big{]}^{T}.$

This is a biased choice which includes more points to the left. Again the
choice covers the domain of dependence. Training data are obtained from the
exact solution

$u(x,t)=\left\\{\begin{array}[]{lr}1,~{}~{}~{}x\leq\frac{1}{2}t,\\\
0,~{}~{}~{}otherwise.\end{array}\right.$ (3.9)

Figure 7: shock propagation (Example 3.3.2) with neural network method

The neural network picked contains 1 hidden layer of 8 neurons. Network
training is conducted for up to $K=10^{5}$ iterations. The squared $L_{2}$
training error of (2.11) to iteration at $t=2$ are output in Fig 5. After well
trained, the neural network is applied solving the Riemann problem to $T=8.0$.
Screen shots of $t=0$, $t=1.5$, $t=3.0$ and $t=6.0$ are shown in Figure 7. The
neural network method can accurately and sharply capture the shock evolution.
Notice the neural network simulation at $t=8$ is way later than the training
data set of $t=2$.

###### Example 3.3.3

_Rarefaction wave_

In this example, we consider a Riemann problem of (3.7) with piece wise
constant initial

$u(x,0)=\left\\{\begin{array}[]{lr}0,~{}~{}~{}x<0,\\\
1,~{}~{}~{}otherwise.\end{array}\right.$

Domain is $D=[-1,5]$. Dirichlet boundary condition $u(-1,t)=0$, $u(5,t)=1$ and
its out of domain extension are applied. Cell size $\Delta x=0.06$ and time
step size $\Delta t=0.1$ are chosen. Network input vector is taken as

$\overrightarrow{V}_{j}^{in}=\Big{[}\bar{u}^{n}_{j-2},\bar{u}^{n}_{j-1},\bar{u}^{n}_{j},\bar{u}^{n}_{j+1}\Big{]}^{T}.$

The choice still includes the domain of dependence. Training data are
generated from exact solution

$u(x,t)=\left\\{\begin{array}[]{lr}0,~{}~{}x<0,\\\ \frac{x}{t},~{}~{}0\leq
x\leq t,\\\ 1,~{}~{}~{}otherwise.\end{array}\right.$

Figure 8: rarefaction wave propagation (Example 3.3.3) with neural network
method

The network contains 2 hidden layers with 8 neurons per layer. Iteration steps
up to $K=10^{5}$ are applied. The squared $L_{2}$ training error of (2.11)
corresponding to iteration at $t=2$ are also included in Figure 5. After well
trained, the neural network solver is applied solving this Riemann problem up
to $T=4.0$. Screen shots are shown in Figure 8. The neural network method can
open up and well resolves the rarefaction wave evolution.

###### Example 3.3.4

_Interaction of rarefaction and shock waves_

In this example, we consider solving (3.7) with three piece wise constants
initial

$u(x,0)=\left\\{\begin{array}[]{lr}0,~{}~{}~{}x<0,\\\ 1,~{}~{}~{}0\leq x\leq
1,\\\ 0,~{}~{}~{}otherwise.\end{array}\right.$

Spatial domain is set as $D=[-1,5]$. Dirichlet boundary conditions
$u(-1,t)=u(5,t)=0$ and its out of domain extension are applied. Network input
vector is taken as

$\overrightarrow{V}_{j}^{in}=\Big{[}\bar{u}^{n}_{j-4},\bar{u}^{n}_{j-3},\bar{u}^{n}_{j-2},\bar{u}^{n}_{j-1},\bar{u}^{n}_{j},\bar{u}^{n}_{j+1},\bar{u}^{n}_{j+2}\Big{]}^{T}.$

Training data pairs $\left(\bar{u}_{j}^{n},\bar{u}_{j}^{n+1}\right)$ are
generated from the exact solution. Here the rarefaction wave travels faster
than the shock and it merges into the shock wave after a while. Before the
meet at $t=2$, the exact solution is given as

$u(x,t)=\left\\{\begin{array}[]{lr}0,~{}~{}~{}x<0,\\\
\frac{x}{t},~{}~{}~{}0\leq x<t,\\\ 1,~{}~{}~{}t\leq x\leq 1+\frac{t}{2},\\\
0,~{}~{}~{}otherwise.\end{array}\right.$

After $t>2$, the rarefaction wave runs into the shock wave. The solution is
composed with zero state to the left, followed by a rarefaction wave and
connected with a shock to the right with which the shock speed slows down as
time evolves. The shock speed is determined by the following formula

$\sigma^{\prime}(t)=\frac{\frac{1}{2}(\frac{\sigma(t)}{t})^{2}}{\frac{\sigma(t)}{t}}=\frac{\sigma(t)}{2t}.$

With initial value $\sigma(2)=2$, we obtain $\sigma(t)=\sqrt{2t}$. For $t>2$,
the exact solution is given by

$u(x,t)=\left\\{\begin{array}[]{lr}0,~{}~{}~{}x<0,\\\
\frac{x}{t},~{}~{}~{}0\leq x\leq\sqrt{2t},\\\
0,~{}~{}~{}\sqrt{2t}<x.\end{array}\right.$

Figure 9: rarefaction wave and shock wave interaction (Example 3.3.4) with
neural network method

The picked network contains 2 hidden layers with 8 neurons per layer. Training
is conducted up to $K=10^{5}$ iterations. Squared $L_{2}$ training error of
(2.11) corresponding to iteration at $t=2$ is listed in Figure 5. After well
trained, neural network method is applied with its solution screen shots
before and after waves interaction shown in Figure 9. Recall the training set
only include solution pairs up to $t=2$ for which the rarefaction wave has not
met the shock wave yet. The neural network method is capable of accurately
capturing shock propagation after the interaction. We mention that we also
check out the $L_{2}$ errors, which are around $O(10^{-3})$ for all four
examples, at a time when solutions all develop singularities. Recall the cell
size is around $\Delta x\approx 0.06$.

###  3.4 Nonlinear convection diffusion equation

In this section, we investigate the viscous Burgers’ equation

$u_{t}+\left(\frac{u^{2}}{2}\right)_{x}=\mu u_{xx},~{}~{}(x,t)\in D\times
R^{+},$ (3.10)

with $\mu=0.1$. Zero boundary conditions of $u(0,t)=u(2\pi,t)=0$ are applied.
Same time step size $\Delta t=0.1$ and cell size $\Delta x=\frac{2\pi}{100}$
as in Example 3.3.1 are taken. The cell-average based neural network method
can adapt large time step size $\Delta t$. We use this setting to simply
illustrate the effectiveness and efficiency of neural network method solving
nonlinear convection diffusion equations. The network input vector is taken
the same as that of Example 3.3.1

$\overrightarrow{V}_{j}^{in}=\Big{[}\bar{u}^{n}_{j-3},\bar{u}^{n}_{j-2},\bar{u}^{n}_{j-1},\bar{u}^{n}_{j},\bar{u}^{n}_{j+1},\bar{u}^{n}_{j+2},\bar{u}^{n}_{j+3}\Big{]}^{T}.$

Training data $\left\\{\bar{u}_{j}^{n},\bar{u}_{j}^{n+1}\right\\}$ are
generated from the refined mesh highly accurate discontinuous Galerkin method.
With time step $\Delta t=0.1$ and up to $t=2.0$, twenty time levels of
solution averages are applied in the training data set. The network structure
contains 2 hidden layers of each with 8 neurons. The network training is
conducted for up to $K=10^{5}$ iterations.

Figure 10: viscous Burgers’ equation (3.10) simulation with neural network
method

Solution evolution run by neural network method is shown in Figure 10. The
results are comparable to those obtained from classical numerical methods.

##  4 Conclusions

A finite volume or cell-average based neural network method is developed for
hyperbolic and parabolic partial differential equations. Classical numerical
methods and principles motivate and guide the design of machine learning
neural network methods. It is found the cell-average based neural network
method can be relieved from the explicit scheme CFL restriction and is able to
evolve the solution forward in time with large time step size. The method is
able to accurately capture shock and rarefaction waves. We believe this is the
right direction to explore neural networks solvers for partial differential
equations, which is built upon solution properties and successful numerical
methods.

## References

  * LeCun and Bengio [1995] Y. LeCun, Y. Bengio,  Convolutional networks for images, speech, and time-series,  The handbook of brain theory and neural networks (1995).
  * Bengio [2009] Y. Bengio,  Learning Deep Architectures for AI,  Found. Trends Mach. Learn. 2 (2009) 1–127.
  * Krizhevsky et al. [2012] A. Krizhevsky, I. Sutskever, G. Hinton,  Imagenet classification with deep convolutional neural networks,  Advances in Neural Information Processing Systems 25 (2012) 1097–1105.
  * LeCun et al. [2015] Y. LeCun, Y. Bengio, G. Hinton,  Deep learning,  Nature 521 (2015) 436–444.
  * E [2017] W. E,  A proposal on machine learning via dynamical systems,  Commun. Math. Stat. 5 (2017) 1–11.
  * Chaudhari et al. [2017] P. Chaudhari, A. Oberman, S. Osher, S. Soatto, G. Carlier, Deep relaxation: partial differential equations for optimizing deep neural networks, 2017\. arXiv:1704.04932.
  * Rudy et al. [2017] S. H. Rudy, S. L. Brunton, J. L. Proctor, J. N. Kutz,  Data-driven discovery of partial differential equations,  Science Advances 3 (2017) e1602614.
  * Chang et al. [2018] B. Chang, L. Meng, E. Haber, L. Ruthotto, D. Begert, E. Holtham,  Reversible architectures for arbitrarily deep residual neural networks,  in: Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, (AAAI-18), 2018, AAAI Press, 2018, pp. 2811–2818.
  * Long et al. [2019] Z. Long, Y. Lu, B. Dong,  PDE-Net 2.0: learning PDEs from data with a numeric-symbolic hybrid deep network,  J. Comput. Phys. 399 (2019) 108925, 17.
  * Ruthotto and Haber [2020] L. Ruthotto, E. Haber,  Deep neural networks motivated by partial differential equations,  J. Math. Imaging Vision 62 (2020) 352–364.
  * He and Xu [2019] J. He, J. Xu,  MgNet: a unified framework of multigrid and convolutional neural network,  Sci. China Math. 62 (2019) 1331–1354.
  * Cybenko [1989] G. Cybenko,  Approximation by superpositions of a sigmoidal function,  Mathematics of Control, Signals and Systems 2 (1989) 303–314.
  * Hornik et al. [1990] K. Hornik, M. Stinchcombe, H. White,  Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks,  Neural Networks 3 (1990) 551 – 560.
  * Barron [1993] A. R. Barron,  Universal approximation bounds for superpositions of a sigmoidal function,  IEEE Transactions on Information Theory 39 (1993) 930–945.
  * Leshno et al. [1993] M. Leshno, V. Y. Lin, A. Pinkus, S. Schocken,  Multilayer feedforward networks with a nonpolynomial activation function can approximate any function,  Neural Networks 6 (1993) 861–867.
  * Pinkus [1999] A. Pinkus,  Approximation theory of the mlp model in neural networks,  Acta Numer. 8 (1999) 143–195.
  * Lagaris et al. [1998] I. Lagaris, A. Likas, D. Fotiadis,  Artificial neural networks for solving ordinary and partial differential equations,  IEEE transactions on neural networks 95 (1998) 987–1000.
  * Rudd and Ferrari [2015] K. Rudd, S. Ferrari,  A constrained integration (cint) approach to solving partial differential equations using artificial neural networks,  Neurocomputing 155 (2015) 277–285.
  * Berg and Nyström [2018] J. Berg, K. Nyström,  A unified deep artificial neural network approach to partial differential equations in complex geometries,  Neurocomputing 317 (2018) 28–41.
  * Sirignano and Spiliopoulos [2018] J. Sirignano, K. Spiliopoulos,  DGM: a deep learning algorithm for solving partial differential equations,  J. Comput. Phys. 375 (2018) 1339–1364.
  * Raissi et al. [2017a] M. Raissi, P. Perdikaris, G. E. Karniadakis,  Physics informed deep learning (part I): data-driven solutions of nonlinear partial differential equations (2017a). arXiv:1711.10561.
  * Raissi et al. [2017b] M. Raissi, P. Perdikaris, G. E. Karniadakis,  Physics informed deep learning (part II): data-driven discovery of nonlinear partial differential equations,  arXiv abs/1711.10566 (2017b).
  * Raissi et al. [2019] M. Raissi, P. Perdikaris, G. E. Karniadakis,  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,  Journal of Computational Physics 378 (2019) 686–707.
  * Lu et al. [2021] L. Lu, X. Meng, Z. Mao, G. E. Karniadakis,  Deepxde: A deep learning library for solving differential equations,  SIAM Review 63 (2021) 208–228.
  * Dwivedi and Srinivasan [2020] V. Dwivedi, B. Srinivasan,  Physics informed extreme learning machine (PIELM)–a rapid method for the numerical solution of partial differential equations,  Neurocomputing 391 (2020) 96–118.
  * Jin et al. [2021] X. Jin, S. Cai, H. Li, G. E. Karniadakis,  Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations,  Journal of Computational Physics 426 (2021) 109951.
  * Zang et al. [2020] Y. Zang, G. Bao, X. Ye, H. Zhou,  Weak adversarial networks for high-dimensional partial differential equations,  Journal of Computational Physics 411 (2020) 109409.
  * Cai et al. [2021a] Z. Cai, J. Chen, M. Liu,  Least-squares relu neural network (lsnn) method for linear advection-reaction equation,  Journal of Computational Physics (2021a) 110514.
  * Cai et al. [2021b] Z. Cai, J. Chen, M. Liu, Least-squares relu neural network (lsnn) method for scalar nonlinear hyperbolic conservation law, 2021b. arXiv:2105.11627.
  * Kutyniok et al. [2021] G. Kutyniok, P. Petersen, M. Raslan, R. Schneider,  A theoretical analysis of deep neural networks and parametric pdes,  Constructive Approximation (2021). URL: https://doi.org/10.1007/s00365-021-09551-4.
  * Michoski et al. [2020] C. Michoski, M. Milosavljević, T. Oliver, D. R. Hatch,  Solving differential equations using deep neural networks,  Neurocomputing 399 (2020) 193–212.
  * Shin et al. [2020] Y. Shin, J. Darbon, G. Em Karniadakis,  On the convergence of physics informed neural networks for linear second-order elliptic and parabolic type pdes,  Communications in Computational Physics 28 (2020) 2042–2074.
  * Laakmann and Petersen [2021] F. Laakmann, P. Petersen,  Efficient approximation of solutions of parametric linear transport equations by ReLU DNNs,  Adv. Comput. Math. 47 (2021) Paper No. 11, 32.
  * Beck et al. [2019] C. Beck, W. E, A. Jentzen,  Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations,  J. Nonlinear Sci. 29 (2019) 1563–1619.
  * Chan-Wai-Nam et al. [2019] Q. Chan-Wai-Nam, J. Mikael, X. Warin,  Machine learning for semi linear PDEs,  J. Sci. Comput. 79 (2019) 1667–1712.
  * Pham et al. [2021] H. Pham, X. Warin, M. Germain,  Neural networks-based backward scheme for fully nonlinear pdes,  SN Partial Differ. Equ. Appl. 2 (2021).
  * Martin et al. [2020] H. Martin, J. Arnulf, K. Thomas, N. T. Anh,  A proof that rectified deep neural networks overcome the curse of dimensionality in the numerical approximation of semilinear heat equations,  SN Partial Differ. Equ. Appl. 1 (2020) 10.
  * Lye et al. [2020] K. O. Lye, S. Mishra, D. Ray,  Deep learning observables in computational fluid dynamics,  Journal of Computational Physics 410 (2020) 109339.
  * Fan et al. [2019] Y. Fan, L. Lin, L. Ying, L. Zepeda-Núñez,  A multiscale neural network based on hierarchical matrices,  Multiscale Model. Simul. 17 (2019) 1189–1213.
  * Khoo et al. [2020] Y. Khoo, J. Lu, L. Ying,  Solving parametric pde problems with artificial neural networks,  European Journal of Applied Mathematics (2020) 1–15.
  * Li et al. [2020] Y. Li, J. Lu, A. Mao,  Variational training of neural network approximations of solution maps for physical models,  J. Comput. Phys. 409 (2020) 109338.
  * Winovich et al. [2019] N. Winovich, K. Ramani, G. Lin,  ConvPDE-UQ: convolutional neural networks with quantified uncertainty for heterogeneous elliptic partial differential equations on varied domains,  J. Comput. Phys. 394 (2019) 263–279.
  * Wu and Xiu [2020] K. Wu, D. Xiu,  Data-driven deep learning of partial differential equations in modal space,  Journal of Computational Physics 408 (2020) 109307.
  * Qin et al. [2021] T. Qin, Z. Chen, J. D. Jakeman, D. Xiu,  Data-driven learning of nonautonomous systems,  SIAM Journal on Scientific Computing 43 (2021) A1607–A1624.
  * Ray and Hesthaven [2018] D. Ray, J. S. Hesthaven,  An artificial neural network as a troubled-cell indicator,  J. Comput. Phys. 367 (2018) 166–191.
  * Wang et al. [2020] Y. Wang, Z. Shen, Z. Long, B. Dong,  Learning to discretize: Solving 1d scalar conservation laws via deep reinforcement learning,  Communications in Computational Physics 28 (2020) 2158–2179.
  * Discacciati et al. [2020] N. Discacciati, J. S. Hesthaven, D. Ray,  Controlling oscillations in high-order discontinuous Galerkin schemes using artificial viscosity tuned by neural networks,  J. Comput. Phys. 409 (2020) 109304, 30.
  * Hsieh et al. [2019] J. Hsieh, S. Zhao, S. Eismann, L. Mirabella, S. Ermon,  Learning neural PDE solvers with convergence guarantees,  arXiv abs/1906.01200 (2019). URL: http://arxiv.org/abs/1906.01200. arXiv:1906.01200.
  * Sun et al. [2020] Z. Sun, S. Wang, L.-B. Chang, Y. Xing, D. Xiu,  Convolution neural network shock detector for numerical solution of conservation laws,  Communications in Computational Physics 28 (2020) 2075–2108.
  * Beck et al. [2020] A. D. Beck, J. Zeifang, A. Schwarz, D. G. Flad,  A neural network based shock detection and localization approach for discontinuous galerkin methods,  Journal of Computational Physics 423 (2020) 109824.
  * Yu and Shu [tted] X. Yu, C.-W. Shu,  Multi-layer perceptron estimator for the total variation bounded constant in limiters for discontinuous galerkin methods,  La Matematica: Official Journal of the Association for Women in Mathematics (submitted).