# Solving Non-local Fokker-Planck Equations by Deep Learning

Senbao Jiang Corresponding Author<EMAIL_ADDRESS>Department of Applied
Mathematics, Illinois Institute of Technology  Xiaofan Li<EMAIL_ADDRESS>Department of Applied Mathematics, Illinois Institute of Technology

###### Abstract

Physics-informed neural networks (PiNNs) recently emerged as a powerful solver
for a large class of partial differential equations under various inital and
boundary conditions. In this paper, we propose _trapz-PiNNs_ , physics-
informed neural networks incoporated with a modified trapezoidal rule recently
developed for accurately evaluating fractional laplacian and solve the space-
fractional Fokker-Planck equations in 2D and 3D. We describe the modified
trapezoidal rule in detail and verify the second-order accuracy. We
demonstrate trapz-PiNNs have high expressive power through predicting solution
with low $\mathcal{L}^{2}$ relative error on a varitey of numerical examples.
We also use local metrics such as pointwise absolute and relative errors to
analyze where could be further improved. We present an effective method for
improving performance of trapz-PiNN on local metrics, provided that physical
observations of high-fidelity simulation of the true solution are available.
Besides the usual advantages of the deep learning solvers such as adaptivity
and mesh-independence, the trapz-PiNN is able to solve PDEs with fractional
laplacian with arbitrary $\alpha\in(0,2)$ and specializes on rectangular
domain. It also has potential to be generlized into higher dimensions.

##  1 Introduction

The Fokker-Planck equations (FPEs) describe the time evolution of probability
density functions of underlying stochastic dynamics [3]. If the driving noise
is Gaussian (brownian motions), the FPE is a parabolic partial differential
equation involving laplacian. The FPEs are widely used in studying stochatic
models in physical, chemical and biological systems. There are many cases that
the driving noises are non-Gaussian [2, 19, 18], such as isotropic
$\alpha$-stable Lévy motions, then the corresponding FPEs are non-local
parabolic PDE involving fractional laplacian. For a introduction to FPEs, see
[3].

We are insterested in the following non-local Fokker-Planck Equation defined
on $\Omega_{T}\coloneqq\Omega\times(0,T)$ where $\Omega\in\mathbb{R}^{n},\
n\geq 1$ is a rectangular domain

$\displaystyle\partial_{t}u(x,t)$ $\displaystyle=-\nabla\cdot(f\
u)+\frac{1}{2}\operatorname{Tr}\left\\{\nabla^{2}:(\sigma(x)\sigma(x)^{T}u(x,t))\right\\}-(-\Delta)^{\frac{\alpha}{2}}u(x,t),\quad(x,t)\in\Omega_{T},$
(1) $\displaystyle u(x,t)$ $\displaystyle=0,\quad x\not\in\Omega,\ t>0,$ (2)
$\displaystyle u(x,0)$ $\displaystyle=u_{0}(x).$ (3)

Here, $\alpha\in(0,2)$, $f(x)$ is a $C^{2}$ vector field in $\mathbb{R}^{n}$
and $\sigma(x)$ is a $C^{2}$ $n\times n$ matrix-valued map, the fractional
laplacian is defined by

$\displaystyle(-\Delta)^{\frac{\alpha}{2}}\varphi(x)$ $\displaystyle\coloneqq
c_{n,\alpha}\text{P.V.}\int_{\mathbb{R}^{n}\setminus\\{0\\}}\frac{\varphi(x)-\varphi(y)}{|x-y|^{n+\alpha}}\operatorname{d}\\!{y},$
(4) $\displaystyle\text{where}\quad c_{n,\alpha}$
$\displaystyle\coloneqq\frac{2^{\alpha}\Gamma(\frac{\alpha+n}{2})}{\pi^{\frac{n}{2}}|\Gamma(-\frac{\alpha}{2})|}.$
(5)

We sometimes denote Eq. 1 by $\partial_{t}u=\mathcal{L}u$ for convenience.
After rescaling, we can always assume $\Omega=(-1,1)^{n}$, a hypercube with
side length two. For more information about fractional laplacian, we refer to
[14].

Finite difference type numerical method for the non-local partial differential
equations involving fractional laplacian usually relies on the discretization
of fractional laplacian [26, 21, 22, 4, 5, 8]. Recently, a new accurate
numerical scheme has been developed for computing fractional laplacian in 2D
[11]. It is a modified trapezoidal rule with correction terms around
singularity in order to achieve high-order accuracy of the quadrature. The
correction weights are pre-computed. Once being stored, the modified
trapezoidal rule can be easily set up. Blended with finite difference method,
it has been applied to solving one and two dimensional non-local FPEs in [5,
8] respectively. It is not hard to generalize this rule to arbitrary $n$
dimensions. For full presentation of the modified trapezoidal rule in 2D and
its convergence analysis, please refer to [11].

Apart from classical methods, the surge of artificial intelligence influenced
numerical PDEs community. The Physics-informed Neural Networks (PiNNs) emerged
as powerful deep learning solvers for partial differential equations (PDEs)
[20, 9, 17], fractional PDEs [15, 16], or stochastic PDEs [25] with various
initial and/or boundary conditions. PiNN was formally introduced by Rassi et
al [17], in which the authors provided data-driven solution (forward problem)
and data-driven discovery (inverse problem) of some integer-order PDEs for
continuous time models and discrete time models. PiNN has been further
extended to solve nonlocal or fractional type PDEs by incorporating classcial
numerical methods to evaluate non-local or fractional operators, where the
automatic differentiation is not applicable. In [15], Pang et al proposed a
unified nonlocal operator encompassing both fractional laplacian and classical
laplacian, which is computed in spherical coordinates and Gauss-Legendre
quadrature rule. Combining this unified operator with PiNN, they are able to
solve non-local Poisson model and non-local turbulence. A parallel work by
Pang et al [16] combines PiNN with one-dimensional and multi-dimensional
Grunwald-Letnikov numerical schemes and Gauss-Legendre quadrature rule to
discretize the fractional laplacian and solves forward and inverse problem of
fractional advection-diffusion equations in one, two and three dimensions in
circular or spherical domains. Xu et al [23] verifies the PiNN accurately
solves the integer-order stationary FPEs in one, two and three dimenions.

In this work we propose the modified-trapezoidal-rule-incoporated Physics-
informed Neural Networks and focus on forward problem for the non-local FPEs.
We call _trapz-PiNNs_ for short, it is sufficient to only use the simplest
version of the modified trapezoidal rule: a trapezoidal rule plus one
correction term on singularity and provides second-order accuracy for
evaluating the fractional laplacian. The expressive power of trapz-PiNN is
demonstrated by accurately predicting the solution of nonlocal FPEs at time
$T_{\text{pred}}$ after training strictly before the time $T$ for some $T\leq
T_{\text{pred}}$. We verify the accuracy of the deep learning (DL) solution by
comparing it with reference solution obtained from finite difference method
(FDM) under $\mathcal{L}^{2}$ relative error. Apart from achieving good global
metric $\mathcal{L}^{2}$ relative error, we investigate the DL solution
profile from local metric such as pointwise absolute and relative errors. It
turns out that even achieving good $\mathcal{L}^{2}$ relative error, trapz-
PiNN predicts more accurately in the regions with large or moderate magnitude
than regions with small magnitude. If physical observation of solution data or
high-fidelity simulated data is available [15], we propse an effective loss
function integrating the extra data into trapz-PiNN framework to furthur
enhance the performance on the local metrics, provided physical observation of
solution data or high-fidelity simulated data is available.

It worths mentioning that trapz-PiNN is not limited to being the solver of the
non-local FPEs. It can be adapted to other non-local PDEs involving fractional
laplacian such as fractional reaction-diffusion equations, Cahn-Hilliard
equations etc. Besides the usual advantages of PiNNs such as independence of
mesh grids and easy adaptivity, the proposed trapz-PiNN has the following
characteristics:

  1. 1.

trapz-PiNN solves the non-local PDEs on rectangular spacial domains while the
previous relevent works emphasize on circular or spherical ones.

  2. 2.

The numerical scheme for fractional laplacian is valid for arbitrary
$\alpha\in(0,2)$, and so is trapz-PiNN.

  3. 3.

Although we only provide numerical examples in two or three dimensions in this
work, trapz-PiNN can be generalized to higher dimensional settings.

We describe some notations that appear throughout the paper. $n$ denotes
spatial dimensions, usually $2,3$ but could be higher. $\Delta t>0$ represents
the time step size and $h>0$ means space resolution. $|x|$ is the $l^{2}$ norm
of $x\in\mathbb{R}^{n}$. $\Omega=(-1,1)^{n}$ is a hypercube with side length
two and $\Omega_{T}=\Omega\times(0,T)$ for some $T>0$. To facilitate reading,
Table 1 records all common acronyms used in this paper.

Full word | acronyms
---|---
Partial Differential Equations | PDEs
Fokker-Planck Equations | FPEs
Ornstein–Uhlenbeck | O-U
Finite Difference Method | FDM
Neural Network | NN
Physics-informed Neural Networks | PiNNs
Mean Squared Error | MSE
Deep Learning | DL
Stochastic Gradient Descent | SGD
Table 1:  Common Acronyms

We organize this paper as follows. In Section 2 we discuss the modified
trapezoidal rule and demonstrate the second-order accuracy for approximating
the fractional laplacian. We elaborate the construction of the trapz-PiNNs in
Section 3. Numerical results are given in detail in Section 4, in particular,
we solve fractional heat equations in 2D and 3D and FPE with
Ornstein–Uhlenbeck potential in 2D. We also proposed a modified loss function
and compare the DL solution profiles computed by trapz-PiNNs equipped with the
original and modified loss functions. Section 5 summarizes the work and
present future directions.

##  2 Numerical Scheme for the Fractional Laplacian

###  2.1 The modified trapezoidal rule

Accurately computing the fractional laplacian is essential in constructing the
loss function in trapz-PiNNs. In this section, we introduce a simplest version
of the modified trapezoidal rule introduced in [11]. Let $\phi\in
C_{c}^{N}(\mathbb{R}^{n})$ with $N>\max\\{4-\alpha,n\\}$. The modified
trapezoidal rule is designed to numerically compute a class of weakly singular
integrals, in particular

$\displaystyle
I_{\alpha}^{(j)}=\int_{\mathbb{R}^{n}\setminus\\{0\\}}\phi(x)\frac{x_{j}^{2}}{|x|^{n+\alpha}}\operatorname{d}\\!{x},\quad
j=1,\cdots,n.$ (6)

It turns out that discretization of fractional laplacian boils down to
evaluating weakly singular integrals (6). See [8, section 3.3]. Let $f$ be a
Schwartz function or function with compact support on $\mathbb{R}^{n}$, we
define punctured-hole trapezoidal rule to be

$\displaystyle T_{h}[f]\coloneqq h^{n}\sum_{k\in\mathbb{Z}^{n},k\not=0}f(kh).$
(7)

Writing

$\displaystyle s_{j}(x)=\frac{x_{j}^{2}}{|x|^{n+\alpha}},$ (8)

the modified trapezoidal rule for $I_{\alpha}^{(j)}$ is given by

$\displaystyle Q_{h}^{(j)}[\phi s_{j}]=T_{h}[\phi
s_{j}]+h^{2-\alpha}\omega_{0}^{(j)}\phi(0).$ (9)

Here, $\omega_{0}^{(j)}$ is a correction weight that can be computed by the
following limit

$\displaystyle\omega_{0}^{(j)}=\lim_{h\to
0}\frac{1}{h^{2-\alpha}}\left(\int_{\mathbb{R}^{n}\setminus\\{0\\}}g(x)s_{j}(x)\operatorname{d}\\!{x}-T_{h}[g\cdot
s_{j}]\right),$ (10)

where $g$ is a radially symmetric Schwartz function. After a symmetry
argument, one can see that

$\displaystyle\omega_{0}^{(1)}=\cdots=\omega_{0}^{(n)}\eqqcolon\omega_{0}.$
(11)

It is a proven fact that this simple modified trapezoidal rule has order of
convergence $4-\alpha$ for weakly singular integral (6). For more information
about the full version of this numerical rule and related convergence
analysis, please refer to [11].

###  2.2 Numerical Fractional Laplacian

The modified trapezoidal rule (9) can be applied to fractional laplacian (4)
with second-order accuracy. Let $\varphi(x)\in C_{c}^{N}(\mathbb{R}^{n})$ with
$supp(\varphi)\subset\Omega$ and let $\Omega_{h}\coloneqq(x_{k})_{k}$ be any
uniform mesh grid for $\Omega$ with mesh size $h$, then the fractional
laplacian evaluated at $x_{k}$ is approximated by

$\displaystyle(-\Delta)^{\frac{\alpha}{2}}\varphi(x_{k})$
$\displaystyle\approx
h^{n}c_{n,\alpha}\sideset{}{{}^{\prime\prime}}{\sum}_{x_{l}\in\Omega_{h}\setminus\\{x_{k}\\}}\frac{\varphi(x_{k})-\varphi(x_{l})}{|x_{k}-x_{l}|^{n+\alpha}}-\frac{1}{2}c_{n,\alpha}\omega_{0}h^{2-\alpha}\sum_{j=1}^{n}\partial_{jj}^{2}\varphi(x_{k})$
$\displaystyle+c_{n,\alpha}\varphi(x_{k})I_{\alpha}(x_{k}),$ (12)

where $\sum_{x_{l}}^{{}^{\prime\prime}}$ indicates the summand is divided by
$2^{m}$ if grid point $x_{l}$ lies on $m$ edges and

$\displaystyle
I_{\alpha}(x)\coloneqq\int_{\Omega^{c}}\frac{\operatorname{d}\\!{y}}{|y-x|^{n+\alpha}}$
(13)

can be evaluated analytically. In particular, the analytical formula for Eq.
13 in 2D is available in [8]. We denote the discretized fractional laplacian
(the RHS of Section 2.2) by $(-\Delta)^{\nicefrac{{\alpha}}{{2}}}_{h}$.

We now verify the second-order accuracy of Section 2.2 in 2D, while the
verification in other dimensions are similar. Let

$\displaystyle\phi(x)$
$\displaystyle\coloneqq(1+x_{1}+2x_{1}^{2})(1+x_{2}^{2})(1-x_{1}^{2})_{+}^{5}(1-x_{2}^{2})_{+}^{5}\in
C_{c}^{4}(\mathbb{R}^{2})$ (14)

so that $supp(\phi)=[-1,1]^{2}$. In 2D the correction weight $\omega_{0}$ is
given by

$\displaystyle\omega_{0}=\left\\{\begin{array}[]{cc}0.960844610589965,&\alpha=0.5\\\
1.950132460000978,&\alpha=1\\\
5.038779739396576,&\alpha=1.5\end{array}\right..$ (18)

We compute $(-\Delta)_{h}^{\nicefrac{{\alpha}}{{2}}}\phi(x)$ at
$x=(0.25,-0.125),(0.875,0.25),(0.375,-0.625)$ for each
$h=\frac{1}{2^{3}},\frac{1}{2^{4}}\cdots,\frac{1}{2^{12}}$. Since the
analytical value of $(-\Delta)^{\nicefrac{{\alpha}}{{2}}}\phi(x)$ is not
available, we examine the absolute value of the difference

$\displaystyle\text{Diff}(h)=|(-\Delta)_{h}^{\nicefrac{{\alpha}}{{2}}}\phi(x)-(-\Delta)_{\frac{h}{2}}^{\nicefrac{{\alpha}}{{2}}}\phi(x)|,$

and the log-log plot between $\text{Diff}(h)$ and the mesh size $h$ to
determine the order of accuracy. Fig. 1 confirms the second-order accuracy of
$(-\Delta)_{h}^{\nicefrac{{\alpha}}{{2}}}$ for $\alpha=0.5,1,1.5$
respectively.

(a) $\alpha=0.5$

(b) $\alpha=1$

(c) $\alpha=1.5$

Figure 1: The accuracy of the approximation
$(-\Delta)_{h}^{\nicefrac{{h}}{{2}}}$ for the fractional laplacian by Section
2.2

##  3 Structure of Physics-informed Neural Networks

In this section we elaborate the structure of the trapz-PiNNs. It consists of
three main parts: training data, neural network solver, loss function and
training algorithm.

Training data:

We use the grid points $\mathcal{T}\coloneqq\\{(x_{j},t_{k}):1\leq j\leq
J^{n},1\leq k\leq K\\}$, where $J=\frac{2}{h}$ and $K=\frac{T}{\Delta t}$, as
training data, $(x_{j})_{j}$, $(t_{k})_{k}$ are uniform meshes for $\Omega$
and $(0,T)$ with mesh sizes $h$ and $\Delta t$, respectively. Denote
$(x_{j})_{j}$ by $\Omega_{h}$. Moreover, we can use validation data to monitor
the training process and avoid overfitting phenomenon [24] and the validation
data can be set as
$\mathcal{V}\coloneqq\Omega_{h}\times\\{T_{\text{pred}}\\}$.

Neural Network solver:

We use the fully-connected feedforward neural network (NN) in this paper,
which is the foundation for all variants of neural networks [24]. It is also
the core part of our NN solver Eq. 21. A NN $u_{nn}$ with input dimension
$d\in\mathbb{N}$, depth $N+1\in\mathbb{N}$ with $N$’s hidden layers, width
$M\in\mathbb{N}$ and output dimension $d^{\prime}\in\mathbb{N}$ is of the form

$\displaystyle u_{nn}(x;\Theta)=A_{N+1}\circ\sigma\circ A_{N}\circ\sigma\circ
A_{N-1}\circ\cdots\circ\sigma\circ A_{1}(x),$ (19)

where

  1. 1.

$A_{1}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{M}$,
$A_{i}:\mathbb{R}^{M}\rightarrow\mathbb{R}^{M},\ 2\leq i\leq N$ and
$A_{N+1}:\mathbb{R}^{M}\rightarrow\mathbb{R}^{d^{\prime}}$ are affine
transformations with matrix representations

$\displaystyle A_{i}(x)$ $\displaystyle=W_{i}x+b_{i},\quad i=1,\dots,N+1,$

  2. 2.

$\sigma:\mathbb{R}\rightarrow\mathbb{R}$ is a non-linear activation function.
With some abuse of notation, we follow the convention that $\sigma$ has the
vectorized form $\sigma(x)=(\sigma(x_{1}),\cdots,\sigma(x_{n})),\ \forall
x\in\mathbb{R}^{n},\ \forall n\in\mathbb{N}$.

  3. 3.

The (learnable) parameters $\Theta$ of NN in Eq. 19 are

$\displaystyle\Theta\coloneqq\\{(W_{i},b_{i}):i=1,\cdots,N+1\\}\cong\mathbb{R}^{(N-1)M^{2}+(N+d+d^{\prime})M+d^{\prime}}.$
(20)

We call an NN has _shape_ $(M,N+1)$ if NN has width $M$ and depth $N+1$. To
solve the FPE Eqs. 1 to 3, it is convenient to write the NN solver as an
ansatz form

$\displaystyle\hat{u}(x,t;\Theta)=\rho(t)\tau(x)u_{nn}(x,t;\Theta)+u_{0}(x),$
(21)

so that $\hat{u}$ satisfies the initial-boundary condition Eqs. 3, LABEL: and
2. Here, $\rho(t)$ is a non-negative, bounded, strictly-increasing
$C^{\infty}([0,\infty))$ function such that $\rho(0)=0$ and $\tau(x)$ is a
$C^{\infty}(\Omega)$ function with $\tau|_{\Omega^{c}}=0$. The input and
output dimensions of neural network $u_{nn}(t,x;\Theta)$ are $n+1$ and $1$
while the width and depth vary. A good choice of functions $\rho(t),\ \tau(x)$
deserves careful consideration. We give a heuristic explanation here. The
Universal Approximation Theorem [10, 12] guarantees that NN with arbitrary
width bounded depth or bounded width arbitrary depth uniformly converges to
continuous function on compact set in $\mathbb{R}^{n}$ up to arbitrary degree
of precision. By the ansatz form Eq. 21, the NN tries to learn

$\displaystyle
u_{nn}(t,x;\Theta)\underset{\Theta}{\longrightarrow}\frac{u(t,x)-u_{0}(x)}{\rho(t)\tau(x)}\eqqcolon\varphi(t,x),$
(22)

where we assume $u\in C(\Omega_{T})$ is the true solution and $u_{0}$ is the
initial condition. Since $\Omega_{T}$ is open, we prefer $\varphi$ can be
continuously extended to the compact set $\overline{\Omega_{T}}$ so that
$\varphi\in C(\overline{\Omega_{T}})$. This is possible if and only if $\rho$
and $\tau$ are chosen such that $\varphi$ has well-defined limiting behaviors
as $x\to\partial\Omega$ and/or $t\to 0+$.

Loss function and training:

The loss function $\mathbf{L}(\Theta)$ is the classic Mean Squared Error
(MSE). We adopt the mini-batch training in this work, the loss function is
therefore evaluated on a randomly selected training batch
$\mathcal{B}\subset\mathcal{T}$ of fixed size in each epoch, i.e.

$\displaystyle\mathbf{L}_{\mathcal{B}}(\Theta)=\frac{1}{|\mathcal{B}|}\sum_{(x_{j},t_{k})\in\mathcal{B}}|\partial_{t}\hat{u}(x_{j},t_{k};\Theta)-\mathcal{L}_{h}\hat{u}(x_{j},t_{k};\Theta)|^{2}.$
(23)

$\mathcal{L}_{h}$ denotes the discretized operator of $\mathcal{L}$ at Eq. 1.
In the RHS of Eq. 23, (partial) derivatives with integer order can be computed
directly by automatic differentiation [7], fractional laplacian is evaluated
by modified trapezoidal rule Section 2.2. The grid points for the trapezoidal
rule coincide with training data $\mathcal{T}$. The NN solver at Eq. 21
already satisfied the initial and boundary conditions. Otherwise, proper loss
functions with respect to initial and boundary conditions have to be added
into Eq. 23. Once we set up the training data, NN solver and loss function, we
are looking for best possible parameters $\Theta^{*}$ that minimize the loss
function, i.e.

$\displaystyle\Theta^{*}=\underset{\Theta}{\arg\min}\ \mathbf{L}(\Theta).$
(24)

To achieve this, we use ADAM, a popular stochastic gradient descent (SGD)
algorithm [13]. Several hyperparameters are required to be initialized before
running the algorithm. They are NN shape, batch size, learning rate, epoch,
time and space resolutions. In addition, a initial guess $\Theta_{0}$ is
provided by Xavier Initialization scheme [6].

We summarize this section into the following Algorithm 1.

Algorithm 1 Training algorithm for trapz-PiNN

NN shape, batch size $B$, learning rate $\eta$, epoch $K$, space resolution
$h$, time resolution $\Delta t$ and initial parameter $\Theta_{0}$.

Initialize $\hat{u}(\cdot;\Theta_{0})$, $\mathcal{T}$ and $\mathcal{V}$.

for $i=0,\cdots,K-1$ do

1\. Randomly select $\mathcal{B}\in\mathcal{T}$ with $|\mathcal{B}|=B$

2\. Compute
$\displaystyle\nabla_{\Theta}\mathbf{L}_{\mathcal{B}}(\Theta)|_{\Theta=\Theta_{i}}=\nabla_{\Theta}\bigg{\\{}\frac{1}{|\mathcal{B}|}\sum_{(x_{j},t_{k})\in\mathcal{B}}|\partial_{t}\hat{u}(x_{j},t_{k};\Theta)-\mathcal{L}_{h}\hat{u}(x_{j},t_{k};\Theta)|^{2}\bigg{\\}}\bigg{|}_{\Theta=\Theta_{i}}$

3\. Update
$\Theta_{i+1}=\Theta_{i}-\eta\nabla_{\Theta}\mathbf{L}_{\mathcal{B}}(\Theta)|_{\Theta=\Theta_{i}}$
$\triangleright$ This is vanilla SGD, we use ADAM in practice

4\. Randomly select $\mathcal{B}^{\prime}\in\mathcal{V}$ with
$|\mathcal{B}^{\prime}|=B$

5\. Record the training error $\mathbf{L}_{\mathcal{B}}(\Theta_{i+1})$ and the
validation error $\mathbf{L}_{\mathcal{B}^{\prime}}(\Theta_{i+1})$.
$\triangleright$ For monitoring the training process

end for

##  4 Numerical Experiments

In this section we demonstrate the expressive power of trapz-PiNNs. In each
numerical examples, the trapz-PiNNs are trained on $\Omega_{T}$ with $T=0.2$
for each equation and predict the solutions at multiple time steps starting
from $t=0.2$. Moreover, we verify the accuracy of DL solution by comparing
with the reference solutions obtained from FDM together with modified
trapezoidal rule Eq. 9 and determining the $\mathcal{L}^{2}$ relative error.
We also investigate the DL solution profile from the pointwise absolute and
relative errors.

Let $\hat{u},u_{\text{ref}}$ be the DL solution and reference solution,
respectively. Recall that $\Omega_{h}=(x_{j})_{j}$ is uniform mesh for
$\Omega$ with mesh size $h$. Given a sample time $t$, the $\mathcal{L}^{2}$
relative error at $t$ is defined by

$\displaystyle\varepsilon=\sqrt{\frac{\sum_{j}|\hat{u}(x_{j},t)-u_{\text{ref}}(x_{j},t)|^{2}}{\sum_{j}|u_{\text{ref}}(x_{j},t)|^{2}}}$
(25)

and for each $x_{j}\in\Omega_{h}$, the absolute and relative errors at $x_{j}$
are defined by

$\displaystyle\epsilon_{\text{abs}}(x_{j})$
$\displaystyle=|\hat{u}(x_{j},t)-u_{\text{ref}}(x_{j},t)|,$ (26)
$\displaystyle\epsilon_{\text{re}}(x_{j})$
$\displaystyle=\left|\frac{\epsilon_{\text{abs}}(x_{j},t)}{u_{\text{ref}}(x_{j},t)}\right|.$
(27)

The activation function $\sigma$, the functions $\rho(t)$ and $\tau(x)$ in the
ansatz (21) are fixed throughout

$\displaystyle\sigma(x)$
$\displaystyle=\frac{1}{1+\exp(-x)},\quad\text{sigmoid function}$ (28)
$\displaystyle\rho(t)$ $\displaystyle=1-\exp(-t),$ (29) $\displaystyle\tau(x)$
$\displaystyle=\left(\prod_{k=1}^{n}(1-x_{k}^{2})_{+}\right)^{\frac{1}{2}},$
(30)

and the hyperparameters will be specified in each example. It worths
mentioning that $\tau(x)$ is chosen a posterior to the initial condition.

To construct FDM reference solution for FPE Eqs. 1 to 3, we use grid points
$\mathcal{T}$ where, if not specified, space resolution $h=\frac{1}{32}$ and
time step $\Delta t=\frac{h^{2}}{4}$. If any, first-order partial derivatives
$\partial_{x_{j}}\\{\cdot\\},j=1,\cdots,n$ are evaluated by classical upwind
scheme, second-order partial derivatives
$\partial_{x_{i}x_{j}}^{2}\\{\cdot\\},i,j=1,\cdots,n$ are discretized by
central differences. Fractional laplacian is discretized by modified
trapezoidal rule Eq. 9. For time evolution, we adopt third-order total-
variation diminishing Runge-Kutta scheme, i.e., given an ODE
$\frac{du}{dt}=R(u)$,

$\displaystyle U^{(1)}$ $\displaystyle=U^{n}+\Delta tR(U^{n}),$ $\displaystyle
U^{(2)}$ $\displaystyle=\frac{3}{4}U^{n}+\frac{1}{4}U^{(1)}+\frac{1}{4}\Delta
tR(U^{(1)}),$ $\displaystyle U^{(3)}$
$\displaystyle=\frac{1}{3}U^{n}+\frac{2}{3}U^{(2)}+\frac{2}{3}\Delta
tR(U^{(2)}),$

where $U^{n}$ denotes numerical solution of $u$ evaluated at time $t=t_{n}$.

As a common practice, we adopt $\mathcal{L}^{2}$ relative error
$\varepsilon=0.01$ as benchmark. DL solution with $\mathcal{L}^{2}$ relative
error $\varepsilon<0.01$ is considered accurate in general. We call this
accuracy by $\varepsilon$-accuracy. Another standard is attaining relative
error $\epsilon_{\text{re}}<1\%$ uniformly over domain. We stress that this is
a very strict critera for DL solvers in general and is hard to achieve in
reality. One primary reason is that the DL solvers are optimized with respect
to MSE, which itself is an indicator of global accuracy. To author’s best
knowledge, without utilizing external information such as incorporating high-
fidelty simulated solution into loss function, there is no known loss function
aiming for such high accuracy of (maximum) relative error.

All numerical experiments are coded in Python with machine learning package
JAX [1] and performed on a Google Colab Pro+ account and the source codes are
available at https://github.com/sjiang23.

###  4.1 Fractional Heat Equation

In both 2D and 3D, we consider the following fractional heat equation with the
initial and boundary conditions:

$\displaystyle\partial_{t}u(x,t)$
$\displaystyle=-(-\Delta)^{\frac{\alpha}{2}}u(x,t),$ (31) $\displaystyle
u(x,0)$
$\displaystyle=c_{n}\left(\prod_{j=1}^{n}(1-x_{j}^{2})_{+}\right)^{4},$ (32)
$\displaystyle u(x,t)$ $\displaystyle=0,\quad x\not\in(-1,1)^{n}.$ (33)

Here $n=2,3$, $c_{2}=(\frac{315}{216})^{2}$ and $c_{3}=(\frac{693}{512})^{3}$,
where $c_{n}$ are chosen so that
$\int_{\mathbb{R}^{n}}u(x,0)\operatorname{d}\\!{x}=1$. It is not necessary to
choose $u(x,0)$ to be probability density and other choices are also feasible.
We solve Eqs. 31 to 33 for $\alpha=0.5,1,1.5$ in 2D and $\alpha=1$ in 3D.

####  4.1.1 2D case

In 2D, the hyperparameters for trapz-PiNNs are chosen as the NN shape
$(20,5)$, the learning rate $10^{-3}$, the batch size $64$, the epoch $2\times
10^{5}$, the space resolution $\nicefrac{{1}}{{32}}$ and the time resolution
$10^{-2}$. Fig. 2(a) records the $\mathcal{L}^{2}$ relative errors
$\varepsilon$ defined in Eq. 25 at times steps
$T+\nicefrac{{k}}{{100}},k=0,\cdots,25$. Fig. 3 presents the FDM reference
solutions, the DL solutions, the absolute error and relative error surfaces of
the DL solutions, evaluated at time $t=0.225$ for $\alpha=0.5,1,1.5$. A
detailed study of FDM numerical solution profile for Eqs. 31 to 33 in 2D can
be found at [8].

From Fig. 2(a), we observe that $\varepsilon$-accuracy is achieved starting
from time $t=0.2$ and is preserved into some time range of future for all
values of $\alpha$, in particular, the $\varepsilon$-accuracies stay less than
0.01 for longer period of time when the value of $\alpha$ is smaller. A
possible contributing factor to this behavior is that the solution of larger
$\alpha$ departs further away from the initial condtions as shown in Figs.
3(a), LABEL:, 3(e), LABEL: and 3(i). It is not surprising that the
$\mathcal{L}^{2}$ relative errors for all values of $\alpha$ are strictly
increasing with respect to time $t$, i.e., the closer to the training space
$\Omega_{T}$, the more accurate the DL solution and vice versa. The
$\mathcal{L}^{2}$ relative error $\varepsilon$ increases at fastest pace in
time when $\alpha=1.5$, comparing with those for $\alpha=0.5,1$. We also note
that trapz-PiNN with the fixed NN shape $(20,5)$ is versatile enough to solve
three different cases simultaneously while achieving reasonable accuracy.

Figs. 3(a), LABEL:, 3(e), LABEL: and 3(i) show the FDM solutions to fractional
heat equations Eqs. 31 to 33 at time $t=0.225$ for $\alpha=0.5,1,1.5$
respectively, while Figs. 3(b), LABEL:, 3(f), LABEL: and 3(j) show the
corresponding DL solutions. It is hard to recognize any difference between the
FDM solutions and the DL solutions from the contour plots. Maximum absolute
errors of the DL solutions are of the order $10^{-3}$ for all $\alpha$. When
$\alpha=0.5$, Fig. 3(c) shows the absolute error at the corners of the domain
$\Omega$ are much larger than that at other regions, Fig. 3(d) shows that that
the a majority of $\Omega$ have relative errors strictly below $10^{-2}$, in
fact, $62.9\%$ of the total area of $\Omega$ have relative errors
$\epsilon_{\text{re}}<1\%$. We also note that the relative errors close to
corners are on the order of $10^{-1}$ and accounts for $5.7\%$ of the area in
total. In constrast, Figs. 3(g), LABEL: and 3(k) show some points in the
central area of $\Omega$ have the largest absolute errors for $\alpha=1,1.5$
respectively, while Figs. 3(h), LABEL: and 3(l) show that these points have
very low relative errors. Furthermore, $82.7\%$ and $71.1\%$ of $\Omega$ have
relative errors $\epsilon_{\text{re}}<1\%$ . Similar to the case of
$\alpha=0.5$, relative errors above $10^{-1}$ appear at region close to
corners and boundaries, making up $0.3\%$ and $2.2\%$ in total area for
$\alpha=1,1.5$ respectively.

In general, high $\varepsilon$-accuracy indicates small relative errors for
majority of the domain $\Omega$. As can be seen from the reference solutions
shown in Figs. 3(a), LABEL:, 3(e), LABEL: and 3(i), the DL solutions at the
regions with high relative errors have very small magnitude, compared with the
interior region where the solutions have larger magnitude and small relative
errors. It implies that, if true solution profile has different scales, trapz-
PiNN predicts the region with moderate or large magnitude better than the
region with small magnitude.

(a) Fractional heat equation in 2D

(b) Fokker-Planck equation with O-U potential in 2D

Figure 2: The $\mathcal{L}^{2}$ relative errors $\varepsilon$ at time steps
$0.2+\nicefrac{{k}}{{100}},k=0,\cdots,25$ for the DL solutions of (a)
fractional heat equations in 2D (Eqs. 31 to 33) and (b) FPEs with O-U
potential in 2D (Eqs. 34 to 36) when $\alpha=0.5,1,1.5$.

(a) FDM solution
$\alpha=0.5$

(b) DL solution
$\alpha=0.5$

(c) $\epsilon_{\text{abs}}$ for $\alpha=0.5$

(d) $\log_{10}(\epsilon_{\text{re}})$ for $\alpha=0.5$

(e) FDM solution
$\alpha=1$

(f) DL solution
$\alpha=1$

(g) $\epsilon_{\text{abs}}$ for $\alpha=1$

(h) $\log_{10}(\epsilon_{\text{re}})$ for $\alpha=1$

(i) FDM solution
$\alpha=1.5$

(j) DL solution
$\alpha=1.5$

(k) $\epsilon_{\text{abs}}$ for $\alpha=1.5$

(l) $\log_{10}(\epsilon_{\text{re}})$ for $\alpha=1.5$

Figure 3: The contour plots of the FDM solutions, the DL solutions, the
absolute and relative errors of the DL solutions for the fractional heat
equations in 2D Eqs. 31 to 33 sampled at $t=0.225$ for $\alpha=0.5,1,1.5$.

####  4.1.2 3D case

In 3D, the hyperparameters for trapz-PiNN are chosen as the NN shape $(20,7)$,
the learning rate $0.001$, the epoch $3\times 10^{5}$, the batch size $100$,
the space resolution $\nicefrac{{1}}{{20}}$ and the time resolution $0.01$.
The $\mathcal{L}^{2}$ relative errors $\varepsilon$ at
$t=0.2+\nicefrac{{k}}{{100}},k=0,\cdots,5$ are recorded at Table 2. Fig. 4
presents the contour plots of the FDM and DL solutions, the absolute and
relative errors of the DL solutions evaluated at $t=0.21$ for the cross-
sections $x_{1}=-0.1$ and $x_{2}=0.85$, respecitvely. Due to the limitation of
RAM, we use moderate space resolution $\frac{1}{20}$ for both FDM and DL
solutions.

From Table 2, we observe that the $\varepsilon$-accuracy ($\varepsilon<0.01$)
is achieved starting from $t=0.2$ and preserved before $t=0.25$. Similar to
the 2D case, the FDM solutions shown in Figs. 4(a), LABEL: and 4(e) and the
corresponding DL solutions shown in Figs. 4(b), LABEL: and 4(f) are
indistinguishable. Fig. 4(c) shows the maximum absolute error
$\epsilon_{\text{abs}}$ attained at some central points in the cross-section
$x_{1}=-0.1$. Fig. 4(d) shows that the majority of the cross-section have
relative error below or around $10^{-2}$. Fig. 4(g) also shows the absolute
error peaks at some central points in the cross-section $x_{2}=0.85$. Fig.
4(h) shows the pointwise relative error $\epsilon_{\text{re}}$ at a
significant part of the area close to the boundary of the cross-section
$x_{2}=0.85$ is on the order of $10^{-1}$ while the other areas is below
$10^{-1.5}\approx 3.2\%$. Furthermore, $76.1\%$ of $\Omega$ have the relative
error $\epsilon_{\text{re}}<3\%$ while the regions with relative error
$\epsilon_{\text{re}}>10\%$ occupy $5.2\%$ in $\Omega$.

We see that moderate space resolution $h$ reduces the size of training data
$\mathcal{T}$ and therefore the range of $\varepsilon$-accuracy, compared with
the 2D counter-part. From Figs. 4(d), LABEL: and 4(h), we find that the cross-
section $x_{1}=-0.1$ has better overall accuracy than that of $x_{2}=0.85$,
possibly since the latter one is closer to boundary of the cube $\Omega$ in
$\mathbb{R}^{3}$. We also know from the FDM solutions shown in Figs. 4(a),
LABEL: and 4(e) that the magnitude of solution profile for the cross-section
$x_{2}=0.85$ is small. It follows again that trapz-PiNN predicts better in the
region with moderate or large magnitude than region with small magnitude.

Table 2: The $\mathcal{L}^{2}$ relative errors $\varepsilon$ at time steps $0.2+\nicefrac{{k}}{{100}},k=0,\cdots,5$ for the DL solutions of fractional heat equations in 3D (Eqs. 31 to 33) with $\alpha=1$ $T_{\text{pred}}$ | $t=0.2$ | $t=0.21$ | $t=0.22$ | $t=0.23$ | $t=0.24$ | $t=0.25$
---|---|---|---|---|---|---
$\mathcal{L}^{2}$ relative errors | 6.856e-3 | 7.311e-3 | 7.838e-3 | 8.458e-3 | 9.196e-3 | 1.001e-2

(a) FDM solution
$x_{1}=-0.1$

(b) DL solution
$x_{1}=-0.1$

(c) $\epsilon_{\text{abs}}$ for $x_{1}=-0.1$

(d) $\log_{10}(\epsilon_{\text{re}})$ for $x_{1}=-0.1$

(e) FDM solution
$x_{2}=0.85$

(f) DL solution
$x_{2}=0.85$

(g) $\epsilon_{\text{abs}}$ for $x_{2}=0.85$

(h) $\log_{10}(\epsilon_{\text{re}})$ for $x_{2}=0.85$

Figure 4: The contour plots of the FDM solutions, the DL solutions, the
absolute and relative errors of the cross-sections $x_{1}=-0.1$ and
$x_{2}=0.85$ of the DL solutions for the fractional heat equations in 3D Eqs.
31 to 33 at $t=0.21$ when $\alpha=1$.

###  4.2 FPEs with Ornstein–Uhlenbeck Potential in 2D

We consider the FPE with O-U potential in 2D:

$\displaystyle\partial_{t}u(x,t)$
$\displaystyle=-(\partial_{1}(x_{1}u(x,t))-\partial_{2}(x_{2}u(x,t)))-(-\Delta)^{\frac{\alpha}{2}}u(x,t),$
(34) $\displaystyle u(x,0)$
$\displaystyle=\left(\frac{315}{216}\right)^{2}((1-x_{1}^{2})(1-x_{2}^{2}))_{+}^{4},$
(35) $\displaystyle u(x,t)$ $\displaystyle=0,\quad x\not\in(-1,1)^{2},$ (36)

with $\alpha=0.5,1,1.5$. The solution $u$ to the Fokker-Planck equation Eqs.
34 to 36 is the probability density function corresponding to the stochastic
differential equations in 2D

$\displaystyle\operatorname{d}\\!{X}_{t}=-X_{t}\operatorname{d}\\!{t}+\operatorname{d}\\!{L}^{\alpha}_{t},\quad
X_{0}=X,$ (37)

where $L_{t}^{\alpha}$ is the $\alpha$-stable Lévy process and random variable
$X$ has distribution as $u(x,0)$, with stochastic process vanishes outside
domain $\Omega$. Without the $\alpha$-stable noise, the determinstic dynamical
system drives the process to the unique stable point, the origin $(0,0)$.
Hyperparameters for trapz-PiNNs are chosen as the NN shape $(20,6)$, the
learning rate $10^{-3}$, the epoch $2\times 10^{5}$, the batch size $64$, the
space resolution $\nicefrac{{1}}{{50}}$ and the time resolution $5\times
10^{-3}$. Fig. 2(b) records the $\mathcal{L}^{2}$ relative errors
$\varepsilon$ at $t=0.2+\nicefrac{{k}}{{100}},k=0,\cdots,10$. Fig. 5 presents
the FDM and DL solutions and the absolute and relative error surfaces of the
DL solutions evaluated at $t=0.2$ for $\alpha=0.5,1,1.5$ respecitvely. For
more information on FDM numerical solution profile for equations Eqs. 34 to 36
in 2D, we refer to [8].

From Fig. 2(b), we observe that the $\varepsilon$-accuracy
($\varepsilon<0.01$) is achieved only at time $t=0.2$ for $\alpha=0.5$ under
the effect of the O-U potential, on the other hand, the $\varepsilon$-accuracy
is preserved up to time $t=0.25$ for $\alpha=1.5$ and up to $t=0.3$ for
$\alpha=1$. $\mathcal{L}^{2}$ relative errors are increasing with respect to
time $t$ for all values of $\alpha$. The rate of increase is highest for
$\alpha=1.5$.

Figs. 5(a), LABEL:, 5(e), LABEL: and 5(i) show the FDM solutions to Fokker-
Planck equation with O-U potential in 2D Eqs. 34 to 36 at time $t=0.2$ for
$\alpha=0.5,1,1.5$ respectively, while Figs. 5(b), LABEL:, 5(f), LABEL: and
5(j) are the corresponding DL solutions. They are hard to tell apart. When
$\alpha=0.5$, Fig. 5(c) shows maximum absolute error is of the order $10^{-2}$
and is attained at some interior points in the domain, Fig. 5(d) shows these
peak points have the relative errors $\epsilon_{\text{re}}$ around
$10^{-1.25}\approx 5.6\%$. Together with the FDM solution profile shown in
Fig. 5(a), we see that the area with the relative errors
$\epsilon_{\text{re}}$ greater than $3\%$ concertrates in the area with small
magnitude (deep blue area) of the solution. When $\alpha=0.5,1$, Figs. 5(g),
LABEL: and 5(k) show the maximum absolute errors are below $5\times 10^{-3}$.
Together with the FDM solution profiles shown in Figs. 5(e), LABEL: and 5(i),
Figs. 5(h), LABEL: and 5(l) show the points attaining maximum absolute error
have relative error below $10^{-1.5}\approx 3.2\%$. Furthurmore, for
$\alpha=1\text{ and }1.5$, more than $76.9\%$ and $80\%$ of the total area has
the relative error below $3\%$ respectively and the area with large relative
errors locates at the area with the solution of small magnitude. They account
for $8.6\%$ and $2.1\%$ of the total area, respectively.

Under the effect of O-U potential, $\varepsilon$-accuracy and range of
$\varepsilon$-accuracy are no longer similar to the case without O-U potential
in Section 4.1.1. O-U potential drives the solution profile more higher at the
origin and for $\alpha<1$, it drives the profile near the boundary of the
domain close to zero [8]. This may explain the trapz-PiNN can predict
accurately the solution for $\alpha=0.5$ at $t=0.2$ only but gives much better
predictions for the cases of $\alpha=1$ and $1.5$.

(a) FDM solution $\alpha=0.5$

(b) DL solution $\alpha=0.5$

(c) $\epsilon_{\text{abs}}$ for $\alpha=0.5$

(d) $\log_{10}(\epsilon_{\text{re}})$ for $\alpha=0.5$

(e) FDM solution $\alpha=1$

(f) DL solution $\alpha=1$

(g) $\epsilon_{\text{abs}}$ for $\alpha=1$

(h) $\log_{10}(\epsilon_{\text{re}})$ for $\alpha=1$

(i) FDM solution $\alpha=1.5$

(j) DL solution $\alpha=1.5$

(k) $\epsilon_{\text{abs}}$ for $\alpha=1.5$

(l) $\log_{10}(\epsilon_{\text{re}})$ for $\alpha=1.5$

Figure 5: The contour plots of the FDM solutions, the DL solutions, the
absolute and relative errors of DL solutions for FPE with O-U potential in 2D
Eqs. 34 to 36 sampled at $t=0.2$ for $\alpha=0.5,1,1.5$.

###  4.3 A modified loss function

From the numerical results in Sections 4.1, LABEL: and 4.2, we find that the
trapz-PiNN is able to predict solutions with good overall accuracy and more
accurate in the region with large or moderate magnitude than the region with
small magnitude, when the true solution profile has different scales. The loss
function MSE defined by Eq. 23 is a metric of global accuracy, rather than a
pointwise one. If physical observations of solutions or high-fidelity
simulated data is available, we propose an effective loss function to address
above-mentioned issue.

Let $U(x,t)$ denote the physical observation or the high-fidelity simulated
data at $(x,t)\in\Omega_{T}$ and $\mathcal{B}$ be a training batch. The new
loss function is defined by

$\displaystyle\mathbf{L}_{\mathcal{B}}^{\text{new}}(\Theta)$
$\displaystyle\coloneqq\frac{\lambda_{1}}{|\mathcal{B}|}\sum_{(x_{j},t_{k})\in\mathcal{B}}|\partial_{t}\hat{u}(x_{j},t_{k};\Theta)-\mathcal{L}_{h}\hat{u}(x_{j},t_{k};\Theta)|^{2}$
$\displaystyle+\frac{\lambda_{2}}{|\mathcal{B}|}\sum_{(x_{j},t_{k})\in\mathcal{B}}\left(\frac{|\hat{u}(x_{j},t_{k};\Theta)-U(x_{j},t_{k})|+\delta}{|U(x_{j},t_{k})|+\delta}\right)^{2},$
(38)

where $\lambda_{1},\lambda_{2}>0$ are constants, $\delta=10^{-6}$ if $U$
vanishes at some points in $\Omega_{T}$, otherwise $\delta=0$. One can see
that the new loss function $\mathbf{L}_{\mathcal{B}}^{\text{new}}$ is a
weighted sum of the loss function $\mathbf{L}_{\mathcal{B}}$ in Eq. 23 and the
‘mean-squared’ pointwise relative errors. To avoid repetitive presentation of
the results in the same nature, we only study the Fokker-Planck equation with
O-U potential when $\alpha=0.5$ in this subsection. We compare the DL
solutions computed by trapz-PiNN with two loss functions through pointwise
absolute and relative errors. To have a fair comparison, we remain using the
same set of hyperparameters chosen in Section 4.2. We refer trapz-PiNNs
equipped with the loss functions $\mathbf{L}_{\mathcal{B}}$ at Eq. 23 and
$\mathbf{L}_{\mathcal{B}}^{\text{new}}$ at Section 4.3 as the original and the
new trapz-PiNN, respectively. The high-fidelity simulated data are the FDM
solutions evaluated at $t=0.01+\nicefrac{{k}}{{100}},\ k=0,\cdots,19$. Fig. 6
shows the comparison between the DL solutions predicted at $t=0.2$ by the
original and new trapz-PiNNs.

From Figs. 6(a), LABEL: and 6(b), we can see that the new trapz-PiNN reduces
the maximum absolute error $\epsilon_{\text{abs}}$ from $0.01$ to $0.006$.
Moreover, as shown in Figs. 6(c), LABEL: and 6(d), the pointwise relative
errors $\epsilon_{\text{re}}$ for the new trapz-PiNN decreases to $3\%$ at
majority of the total area from the order of $10\%$ corresponding to the
original trapz-PiNN. More precisely, $61.3\%$ of the total area has the
pointwise relative error $\epsilon_{\text{re}}$ below $3\%$. In comparison,
only $40.7\%$ of the area achieves such accuracy for the original trapz-PiNN.
The $\mathcal{L}^{2}$ relative error $\varepsilon$ for the DL solution
predicted by the new trapz-PiNN is $5.8\times 10^{-3}$ while the original
trapz-PiNN attains $\varepsilon=9.9\times 10^{-3}$. The range of
$\varepsilon$-accuracy of the new trapz-PiNN extended from $t=0.2$ to at least
beyond $t=0.21$, as the $\mathcal{L}^{2}$ relative error at $t=0.21$ is
$\varepsilon=6.1\times 10^{-3}$.

The new loss function improves the performance of the trapz-PiNN on pointwise
absolute and relative errors and extends the range of $\varepsilon$-accuracy.
Smaller maximum relative error can be achieved if larger NN shape and longer
epoch are chosen. If physical observations or high-fidelity simulation is
available, using the new loss function can achieve higher global and local
accuracies.

(a) $\epsilon_{\text{abs}}$ \- original

(b) $\epsilon_{\text{abs}}$ \- new

(c) $\log_{10}(\epsilon_{\text{re}})$ \- original

(d) $\log_{10}(\epsilon_{\text{re}})$ \- new

Figure 6: The contour plots of the absolute and relative errors of the DL
solutions computed by the original and the new trapz-PiNN, for FPE with O-U
potential in 2D Eqs. 34 to 36 sampled at $t=0.2$ for $\alpha=0.5$.

##  5 Conclusion

We propose trapz-PiNN, a new physics-informed neural network, based on a
recently developed modified trapezoidal rule, to solve non-local Fokker-Planck
equations involving fractional laplacian. We have presented the simplest
version of the modified trapezoidal rule in $\mathbb{R}^{n}$ and have verified
second order accuracy for computing fractional laplacian in 2D. We have
demonstrated trapz-PiNNs have high expressive power through numerical examples
on fractional heat equations in 2D and 3D and Fokker-Planck eqution with O-U
potential in 2D. The DL solutions that has low $\mathcal{L}^{2}$ relative
error ($\varepsilon$-accuracy) in general garantuee a small pointwise relative
error for areas with high or moderate magnitude in $\Omega$. We also observe
that trapz-PiNNs have some range of $\varepsilon$-accuracy for almost all
cases we studied. If physical oberservation or high-fidelity simulation of
true solution is available, we propose an effective loss function integrating
the extra information so that trapz-PiNN improves the performance on local
metrics such as pointwise absolute and relative errors.

There are numerous questions remain to be investigated in future. Without
resorting to using more brute force training or adopting external information,
we will focus on designing an effective loss function to control the relative
errors from region with vanishing magnitude, especially for PDEs with multi-
scale solution profile. Due to RAM restriction, we only use moderate space
resolution in solving fractional heat equation in 3D at Section 4.1.2. Low-
memory fast algorithm for numerical fractional laplacian by modified
trapezoidal rule has been developed and implemented in [21, 8]. One can
incorporate this fast algorithm into deep learning algorithm to alleviate the
“curse of dimensionality”. We are also interested in developing efficient
trapz-PiNNs for Fokker-Planck equation Eq. 1 with discontinuous initial
condition at boundary or with natural boundary (unbounded) condition.

##  6 Acknowledgement

The authors thanks Yiwei Wang for very helpful discussion about this work.

## References

  * [1] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang.  JAX: composable transformations of Python+NumPy programs, 2018\.
  * [2] P. D. Ditlevsen.  Observation of alpha-stable noise induced millennial climate changes from an ice-core record.  In AGU Fall Meeting Abstracts, volume 2005, pages NG23B–0097, December 2005.
  * [3] J. Duan.  An Introduction to Stochastic Dynamics.  Cambridge Texts in Applied Mathematics. Cambridge University Press, 2015\.
  * [4] Siwei Duo and Yanzhi Zhang.  Accurate numerical methods for two and three dimensional integral fractional laplacian with applications.  Computer Methods in Applied Mechanics and Engineering, 355:639–662, 2019.
  * [5] Ting Gao, Jinqiao Duan, and Xiaofan Li.  Fokker–Planck equations for stochastic dynamical systems with symmetric Lévy motions.  Applied Mathematics and Computation, 278(Supplement C):1 – 20, 2016\.
  * [6] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pages 249–256. JMLR Workshop and Conference Proceedings, 2010.
  * [7] A. Griewank and A. Walther.  Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Second Edition.  Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104), 2008.
  * [8] Hansen Ha.  Non-Gaussian Stochastic Dynamics with Deterministic Numerical Tools.  PhD thesis, 2018.  Copyright - Database copyright ProQuest LLC; ProQuest does not claim copyright in the individual underlying works; Last updated - 2018-11-22.
  * [9] Jiequn Han, Arnulf Jentzen, and Weinan E.  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences, 115(34):8505–8510, 2018.
  * [10] Kurt Hornik.  Approximation capabilities of multilayer feedforward networks.  Neural Networks, 4(2):251–257, 1991.
  * [11] Senbao Jiang and Xiaofan Li.  Arbitrarily high-order trapezoidal rules for functions with fractional singularities in two dimensions.  Applied Mathematics and Computation, 429:127236, 2022.
  * [12] Patrick Kidger and Terry Lyons.  Universal approximation with deep narrow networks.  In Conference on learning theory, pages 2306–2327. PMLR, 2020.
  * [13] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [14] Anna Lischke, Guofei Pang, Mamikon Gulian, Fangying Song, Christian Glusa, Xiaoning Zheng, Zhiping Mao, Wei Cai, Mark M. Meerschaert, Mark Ainsworth, and George Em Karniadakis.  What is the fractional laplacian? a comparative review with new results.  Journal of Computational Physics, 404:109009, 2020.
  * [15] G. Pang, M. D’Elia, M. Parks, and G.E. Karniadakis.  nPINNs: Nonlocal physics-informed neural networks for a parametrized nonlocal universal laplacian operator. algorithms and applications.  Journal of Computational Physics, 422:109760, 2020.
  * [16] Guofei Pang, Lu Lu, and George Karniadakis.  fPINNs: Fractional physics-informed neural networks.  SIAM Journal on Scientific Computing, 41(4):A2603–A2626, 2019.
  * [17] M. Raissi, P. Perdikaris, and G.E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [18] James A Roberts, Tjeerd W Boonstra, and Michael Breakspear.  The heavy tail of the human brain.  Current opinion in neurobiology, 31:164–172, 2015.
  * [19] Xiaorong Shen, Hao Zhang, Yong Xu, and Suxia Meng.  Observation of alpha-stable noise in the laser gyroscope data.  IEEE Sensors Journal, 16(7):1998–2003, 2016.
  * [20] Justin Sirignano and Konstantinos Spiliopoulos.  Dgm: A deep learning algorithm for solving partial differential equations.  Journal of Computational Physics, 375:1339–1364, 2018.
  * [21] Hong Wang and Treena S. Basu.  A fast finite difference method for two-dimensional space-fractional diffusion equations.  SIAM Journal on Scientific Computing, 34(5):A2444–A2458, 2012.
  * [22] Hong Wang and Ning Du.  A fast finite difference method for three-dimensional time-dependent space-fractional diffusion equations and its efficient implementation.  Journal of Computational Physics, 253:50–63, 2013.
  * [23] Yong Xu, Hao Zhang, Yongge Li, Kuang Zhou, Qi Liu, and Jürgen Kurths.  Solving Fokker-Planck equation using deep learning.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(1):013133, 2020.
  * [24] Aston Zhang, Zachary C. Lipton, Mu Li, and Alexander J. Smola.  Dive into deep learning.  arXiv preprint arXiv:2106.11342, 2021.
  * [25] Dongkun Zhang, Ling Guo, and George Em Karniadakis.  Learning in modal space: Solving time-dependent stochastic pdes using physics-informed neural networks.  SIAM Journal on Scientific Computing, 42(2):A639–A665, 2020.
  * [26] Meng Zhao, Hong Wang, and Aijie Cheng.  A fast finite difference method for three-dimensional time-dependent space-fractional diffusion equations with fractional derivative boundary conditions.  Journal of Scientific Computing, 74(2):1009–1033, 2018.