#
A Deep Learning-based Collocation Method for
Modeling Unknown PDEs from Sparse Observation

Priyabrata Saha  Saibal Mukhopadhyay

###### Abstract

Deep learning-based modeling of dynamical systems driven by partial
differential equations (PDEs) has become quite popular in recent years.
However, most of the existing deep learning-based methods either assume strong
physics prior, or depend on specific initial and boundary conditions, or
require data in dense regular grid making them inapt for modeling unknown PDEs
from sparsely-observed data. This paper presents a deep learning-based
collocation method for modeling dynamical systems driven by unknown PDEs when
data sites are sparsely distributed. The proposed method is spatial dimension-
independent, geometrically flexible, learns from sparsely-available data and
the learned model does not depend on any specific initial and boundary
conditions. We demonstrate our method in the forecasting task for two
dimensional wave equation and Burgers-Fisher equation in multiple geometries
with different boundary conditions.

††footnotetext: Authors are with School of Electrical and Computer
Engineering, Georgia Institute of Technology, Atlanta, USA.††footnotetext:
Correspondence: Priyabrata Saha<EMAIL_ADDRESS>
##  1 Introduction

Complex physical phenomena are often modeled with a set of differential
equations (DEs). Traditionally, model equations are constructed based on
physical laws and empirical observation. However, in many practical scenarios,
mathematical modeling is very challenging and/or can only describe the system
behavior partially. Furthermore, even if the involved differential equations
are known, their analytical solutions are difficult to obtain for most
practical problems and therefore, must be solved numerically [1]. Traditional
accurate numerical solvers are computationally intensive. The success of deep
learning in various sequence prediction tasks, for example natural language
and speech modeling [2, 3, 4], video prediction [5, 6, 7], has motivated many
researchers to approach both the aforementioned problems using deep learning.
Consequently, a number of deep learning models which incorporate knowledge
from physics or applied mathematics, have been introduced for modeling complex
dynamical systems involving partial differential equations (PDEs) [8, 9, 10,
11, 12, 13, 14, 15, 16, 17].

Typically deep learning methods requires huge amount of data to obtain
generalized solution. However, in many practical scenarios, data or
observation can be obtained only in few sites (for example, sensor
measurements). In this paper, we consider the problem of modeling unknown PDE-
driven spatiotemporal dynamics from sparsely-observed data. Given a set of $N$
distinct scattered spatial locations
$X=\\{\mathbf{x}_{1},\cdots,\mathbf{x}_{N}\\}\subset\Omega\subset\mathbb{R}^{d}$
and corresponding time series of measurements of a physical quantity at those
locations $\\{u(t,\mathbf{x}_{i})\in\mathbb{R}:\mathbf{x}_{i}\in X\
\text{and}\ t=t_{0},t_{0}+\Delta t,t_{0}+2\Delta t,\cdots\\}$, we want to
learn a model for the unknown PDE that governs the data. We assume the
underlying PDE belongs to a family of PDEs described by the following generic
form:

$\partial_{t}^{p}u=F(u,\partial_{\mathbf{x}}u,\partial_{\mathbf{x}}^{2}u,\cdots),\quad\mathbf{x}\in\Omega\subset\mathbb{R}^{d},\quad
t>0$ (1)

$\partial^{p}_{t}$ denotes the $p^{th}$ order partial differential operator
with respect to time and $\partial^{q}_{\mathbf{x}}$ denotes all possible
combination of $q^{th}$ order spatial (partial) differential operators. We
assume that the temporal order ($p$) of the underlying PDE is known a priori
while the form of the function $F$, nonlinear in general, and its arguments
are unknown. Though we observe the system at few discrete locations, we seek
to predict the evolution of $u$ at any generic location $\mathbf{x}\in\Omega$
or over the entire $\Omega$. We want the learned model to construct the
solution of the underlying PDE in time marching fashion, given any initial and
boundary conditions.

Our method integrates radial basis function (RBF) collocation method with deep
learning for modeling unknown PDEs from sparsely-observed data. RBF
collocation method is widely used as numerical solver for PDEs because of its
inherent meshfree nature that provides geometric flexibility with scattered
node layout [18]. Another advantage of RBF-based methods is that they can be
used to solve high dimensional problem since they convert problems involving
multiple space dimensions into virtually one-dimensional [19, 20]. We propose
to learn the application of the differential operators on a chosen RBF using
deep neural networks which makes the method independent of initial and
boundary conditions while also utilizing the meshfree, geometrically flexible
and spatial dimension-independent nature of RBF collocation method. By
choosing infinitely-smooth RBFs (e.g. Gaussian, multiquadric etc.) and tuning
their shape parameters, our method can be employed in scenarios where data
sites are sparsely distributed. We first show how RBF collocation method for
time-dependent linear PDEs can be integrated with deep neural network to
incorporate unknown linear PDEs. Next, we extend the architecture to model
unknown nonlinear PDEs of form (1) and finally to model coupled system of
multiple nonlinear PDEs. We evaluate our method in the forecasting task for
two dimensional wave equation and Burgers-Fisher equation, and demonstrate the
capability of transferring the learned model in different geometries with
different boundary conditions.

#### Related work.

Deep learning methods for learning or solving PDE-driven system can be broadly
categorized into two types. Fhe first kind of methods directly approximates
the PDE solution with deep neural networks (DNNs) as a function of space and
time [8, 13, 21]. These methods are meshfree, geometrically flexible and
suitable for high dimensional PDEs. However, since these methods directly
approximate the PDE solution as a function of space and time, the learned
model is dependent on the specific initial and boundary condition used in
training.

The second type of methods learns the differential operators using
convolutional neural networks (CNNs) [9, 11, 22, 17]. These methods are
independent of initial and boundary conditions, but they are not suitable for
high-dimensional PDEs and can only be used when data is available in a dense
regular grid. Furthermore, both types of methods often assume strong physics
prior either in model structure [10, 11, 12, 16] or in loss function [8, 13,
15] limiting their applicability in scenarios where the dynamics is unknown.

One promising approach for solving or learning PDEs using irregular scattered
data is to utilize graph neural networks (GNNs) [23]. Seo and Liu [24]
proposed differentiable physics-informed graph networks (DPGN) which
incorporates differentiable physics equations in GNN to model complex
dynamical systems. However, difference operators on graph can extract spatial
variations only when the nodes are close to each other and contribute to large
numerical error when data sites are sparsely distributed [25].

##  2 Background

###  2.1 Radial Basis Functions (RBFs) Interpolation

In standard scattered data interpolation problem we are generally given with a
set of $N$ distinct data points
$X=\\{\mathbf{x}_{1},\cdots,\mathbf{x}_{N}\\}\subset\Omega\subset\mathbb{R}^{d}$
and corresponding real values (measurements) $y_{i},i=1,\cdots,N$. The goal is
to find a (continuous) function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ that
satisfies the interpolation conditions

$f(\mathbf{x}_{i})=y_{i},\quad i=1,\cdots,N$ (2)

A common approach to solve this scattered data interpolation problem is to
consider the function $f(\mathbf{x})$ as a linear combination of basis
functions of certain class. Often these basis functions are formed using
radial functions. Value of a radial function
$\phi:\mathbb{R}^{d}\rightarrow\mathbb{R}$ at each point depends only on the
some arbitrary distance between that point and origin such that
$\phi(\mathbf{x})=\phi(\mathbb{\|\mathbf{x}\|})$ where $\|\cdot\|$ denotes
some norm in $\mathbf{R}^{d}$, usually the Euclidean norm. Effectively, the
function becomes a univariate function
$\phi:\mathbb{R}^{+}\rightarrow\mathbb{R}$.

In radial basis functions interpolation, the basis functions are formed using
the given datapoints $\mathbf{x}_{i},i=1,\cdots,N$ as centers or origins. In
terms of these basis functions $\phi(\|\cdot\|)$, the interpolant takes the
following form:

$f(\mathbf{x})=\sum_{j=1}^{N}c_{j}\phi(\|\mathbf{x}-\mathbf{x}_{j}\|),\quad\mathbf{x}\in\Omega$
(3)

where, $c_{j},j=1,\cdots,N$ are the unknown coefficients to be determined.
Plugging the interpolation conditions (2) into (3) leads to a system of linear
equations

$\mathbf{f}=\Phi\mathbf{c},$ (4)

where $\Phi$ is the interpolating matrix whose entries are given by
$\Phi_{ij}=\phi(\|\mathbf{x}_{i}-\mathbf{x}_{j}\|),\ i,j=1,\cdots,N$,
$\mathbf{c}$ is the coefficient vector $[c_{1},\cdots c_{N}]^{T}$ and
$\mathbf{f}=[f(\mathbf{x}_{1}),\cdots,f(\mathbf{x}_{N})]^{T}$.

Solution to the system (4) exists and is unique provided the interpolating
matrix $\Phi$ is non-singular. For a large class of radial functions including
(inverse) multiquadrics, Gaussian, the matrix $\Phi$ is non-singular if the
datapoints are all distinct [19]. Details regarding the invertibility of the
interpolating matrix for various choices of the radial functions can be found
in [26, 19].

###  2.2 Solving time-dependent PDEs using RBFs

Radial basis functions are widely used in meshfree methods for solving PDE
problems. Here we will describe how RBFs are used, particularly Kansa’s
unsymmetric collocation method [27], to solve time-dependent linear PDEs on
which we will base our method.

Consider the following time-dependent linear PDE:

$\partial_{t}u=Lu,\quad\mathbf{x}\in\Omega\subset\mathbb{R}^{d},\quad t>0,$
(5)

with initial and boundary conditions

$u(0,\mathbf{x})=u^{(0)}(\mathbf{x}),\quad\mathbf{x}\in\Omega,$ (6)
$u(t,\mathbf{x})=g(t,\mathbf{x})\quad\mathbf{x}\in\partial\Omega,\quad t>0,$
(7)

$L:\mathbb{R}\rightarrow\mathbb{R}$ is some linear (spatial) differential
operator, $\partial\Omega$ denotes the boundary of $\Omega$.

In order to obtain an approximate solution of the above PDE, consider
$X=\\{\mathbf{x}_{1},\cdots,\mathbf{x}_{N}\\}\subset\Omega$ be a set of $N$
collocation nodes. Furthermore, let us assume that out of the these $N$
collocation nodes, the last $b$ nodes are boundary nodes, i.e. the subset
$X_{B}=\\{\mathbf{x}_{N-b+1},\cdots,\mathbf{x}_{N}\\}\subset\partial\Omega$.
Using some radial basis functions $\phi(\|\cdot\|)$ centered at the
collocation nodes, an approximate solution to the PDE can be defined as
follows:

$\hat{u}(t,\mathbf{x})=\sum_{j=1}^{N}c_{j}(t)\phi(\|\mathbf{x}-\mathbf{x}_{j}\|),\quad\mathbf{x}\in\Omega,\quad
t\geq 0$ (8)

where $c_{j}(t),j=1,\cdots,N$ are the unknown time-dependent coefficients to
be determined at each time step. Plugging (8) into (4), (6), and (7) we get
the following ordinary differential equation (ODE)

$\sum_{j=1}^{N}\frac{dc_{j}(t)}{dt}\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)-c_{j}(t)L\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)=0,\quad\mathbf{x}\in\Omega,\quad
t>0,$ (9)

with initial and boundary conditions

$\sum_{j=1}^{N}c_{j}(0)\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)=u^{(0)}(\mathbf{x}),\quad\mathbf{x}\in\Omega,$
(10)
$\sum_{j=1}^{N}c_{j}(t)\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)=g(t,\mathbf{x}),\quad\mathbf{x}\in\partial\Omega,\quad
t>0,$ (11)

In (9), the spatial differential operator $L$ is applied to the radial basis
function. Approximating the time-derivatives $\frac{dc_{j}(t)}{dt}$ by first
order finite difference, we get

$\sum_{j=1}^{N}c_{j}(t+\Delta
t)\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)-c_{j}(t)\Big{(}\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)+\Delta
tL\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)\Big{)}=0,\quad\mathbf{x}\in\Omega,\quad
t\geq 0,$ (12)

where $\Delta t>0$ is the step-size in time. The coefficients $c_{j}(t)$ can
be determined by applying (11) and (12) on the boundary collocation nodes of
$X_{B}$ and interior collocation nodes of $X\setminus X_{B}$. In that case,
(10), (11) and (12) can be written as the following compact form:

$\displaystyle\Phi\mathbf{c}^{(t+\Delta t)}$
$\displaystyle=A\mathbf{c}^{(t)}+\mathbf{g}^{(t+\Delta t)},\quad t\geq 0$
$\displaystyle\Phi\mathbf{c}^{(0)}$ $\displaystyle=\mathbf{u}^{(0)}$ (13)

$\Phi$ is the RBF matrix whose entries are given by
$\Phi_{ij}=\phi(\|\mathbf{x}_{i}-\mathbf{x}_{j}\|)$, $\ i,j=1,\cdots,N$ and
$\mathbf{c}^{(t)}=[c_{1}(t),\cdots c_{N}(t)]^{T}$ is the coefficient vector at
time $t$. The elements of the matrix $A$ are given by

$\displaystyle
A_{ij}=\begin{cases}\phi(\|\mathbf{x}_{i}-\mathbf{x}_{j}\|)+\Delta
tL\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)|_{\mathbf{x}=\mathbf{x}_{i}},&if\quad
1\leq i\leq(N-b),\ 1\leq j\leq N\\\ 0,&otherwise\ (\text{i.e. }(N-b)<i\leq N,\
1\leq j\leq N)\end{cases}$ (14)

and,
$\mathbf{g}^{(t)}=[0,\cdots,0,g(t,\mathbf{x}_{N-b+1}),\cdots,g(t,\mathbf{x}_{N})]$
is the boundary condition vector at time $t$.
$\mathbf{u}^{(0)}=[u^{(0)}(\mathbf{x}_{1}),\cdots,u^{(0)}(\mathbf{x}_{N})]^{T}$
is the vector of initial values at the collocation nodes. The coefficient
vectors $\mathbf{c}^{(t)}$ are obtained by iteratively solving (13) which are
used to compute the numerical solution $\hat{u}(t,\mathbf{x})$ by (8).

##  3 Deep RBF Collocation

(a)

(b)

(c)

Figure 1: (a): Proposed model for linear PDE; $\Phi_{\sigma}$ is the RBF
interpolation matrix with trainable shape parameter $\sigma$; the unknown
linear differential operator $L$ is modeled using a feed-forward neural
network $L_{\theta}$ with parameter vector $\theta$. (b): Proposed model for
nonlinear PDE; $L_{\theta,h}$ is same as $L_{\theta}$ except it has $h$ output
neurons instead of just one; $F_{\vartheta}$ is the neural network
representation (with parameter vector $\vartheta$) of the unknown function
$F$. (c): Proposed model for a system of $M$ coupled PDEs.

We will first describe our approach for unknown linear PDE and then extend
that to nonlinear PDE of form (1). Finally, we will show how the proposed
model can be used to learn an unknown systems of coupled nonlinear PDEs.

###  3.1 Modeling unknown linear PDE

Consider the linear time-dependent PDE of (5) where the linear (spatial)
differential operator $L$ is unknown. Learning the application of $L$ directly
on $u$ requires fine-grid measurement for convolution operations. Rather, we
propose to learn the application of $L$ on a radial basis function $\phi$ that
allows scattered sparse measurements independent of any specific geometry.
Specifically, we propose to learn a feed-forward neural network $L_{\theta}$,
with parameter vector $\theta$, that approximates
$L\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)|_{\mathbf{x}=\mathbf{x}_{i}}$ given any
two spatial locations $\mathbf{x}_{i},\mathbf{x}_{j}$ and corresponding value
of the RBF $\phi(\|\mathbf{x}_{i}-\mathbf{x}_{j}\|)$. We train the neural
network $L_{\theta}$ using scattered time-series measurement of $u$ by
minimizing an objective function derived from (13). Let
$X=\\{\mathbf{x}_{1},\cdots,\mathbf{x}_{N}\\}\subset\Omega$ be the set of
measurement sites and
$\mathbf{u}^{(t)}=[u(t,\mathbf{x}_{1}),\cdots,u(t,\mathbf{x}_{N})]^{T}$ be the
measurement vector at time $t$. Since our goal is to learn a model for spatial
differential operator $L$ and boundary values are not effected by $L$, we only
consider internal measurement sites during training. For a radial basis
function $\phi$, the coefficient vector at time $t$, $\mathbf{c}^{(t)}$, can
be obtained by solving the linear system

$\Phi\mathbf{c}^{(t)}=\mathbf{u}^{(t)}$ (15)

where $\Phi$ is the RBF matrix whose entries are given by
$\Phi_{ij}=\phi(\|\mathbf{x}_{i}-\mathbf{x}_{j}\|)$, $\ i,j=1,\cdots,N$.
Replacing $L$ with $L_{\theta}$, ignoring boundary nodes, and plugging (15)
into (13) we get

$\mathbf{u}^{(t+\Delta t)}=\mathbf{u}^{(t)}+\Delta
tD_{\theta}\Phi^{-1}_{\sigma}\mathbf{u}^{(t)}=(I+\Delta
tD_{\theta}\Phi^{-1}_{\sigma})\mathbf{u}^{(t)}$ (16)

where $I$ is the identity matrix of order $N$, $D_{\theta}$ denotes the
approximated spatial derivative matrix, i.e.
$[D_{\theta}]_{ij}=L_{\theta}\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)|_{\mathbf{x}=\mathbf{x}_{i}}$,
$\ i,j=1,\cdots,N$ and, $\sigma$ is the trainable shape parameter of the RBF.
So, for a time-series of measurements
$\mathbf{u}^{(t)}\in\mathbb{R}^{n},t=t_{0},t_{0}+\Delta t,\cdots,t_{0}+K\Delta
t$, the neural network $L_{\theta}$ is trained with the following objective:

$\min_{\theta,\sigma}\frac{1}{K}\sum_{k=0}^{K-1}\big{\|}\mathbf{u}^{(t_{0}+(k+1)\Delta
t)}-\mathbf{u}^{(t_{0}+k\Delta t)}-\Delta
tD_{\theta}\Phi^{-1}_{\sigma}\mathbf{u}^{(t_{0}+k\Delta t)}\big{\|}^{2}$ (17)

Note, the term $D_{\theta}\Phi^{-1}\mathbf{u}^{(t)}$ in (16) basically
multiplies the approximated spatial derivative matrix $D_{\theta}$ with
coefficient vector $\mathbf{c}^{(t)}$ to obtain the radial approximation of
$L\mathbf{u}^{(t)}$. Once the model $L_{\theta}$ is trained, it can be used to
get the value of $L_{\theta}u(t,\mathbf{x})$ (Figure 1a) at any generic
location $\mathbf{x}\in\Omega$ and to finally predict $u(t,\mathbf{x})$.
During prediction phase, the boundary condition is applied in similar fashion
as in (13). Considering the last $b$ nodes in $X$ as boundary nodes and
incorporating corresponding boundary condition of (11) in (16), we get

$\mathbf{u}^{(t+\Delta t)}=H\mathbf{u}^{(t)}+\mathbf{g}^{(t+\Delta t)}$ (18)

where the rows of matrix $H$ are given by

$\displaystyle H_{i}=\begin{cases}[I+\Delta
tD_{\theta}\Phi_{\sigma}^{-1}]_{i}&if\quad 1\leq i\leq(N-b)\\\
\mathbf{0}&otherwise\end{cases}$ (19)

Once the model is trained, during prediction the matrix $H$ is solely
determined from the configuration of measurement sites. Under a fixed
Dirichlet boundary condition $g(t,\mathbf{x})=g(\mathbf{x})$, the iterative
scheme of (18) is basically a linear discrete-time dynamical system and its
convergence to a equilibrium point can be analyzed from the spectral
properties of $H$.
The linear discrete-time dynamical system $\mathbf{u}^{(t+\Delta
t)}=H\mathbf{u}^{(t)}+\mathbf{g}$ converges to a stable equilibrium point from
any arbitrary initialization if and only if $|\lambda_{H}|<1$, where
$\lambda_{H}$ are the eigenvalues of $H$.

###  3.2 Modeling unknown nonlinear PDE

In order to model an unknown nonlinear PDE of form (1), we use the proposed
model for linear PDE as the base module. The input arguments of the nonlinear
function $F$ in (1) are basically various unknown linear spatial derivatives
of the measurement variable $u$. Therefore, we first use the linear model to
generate a spatial derivative feature vector and then pass that through
another feed-forward neural network $F_{\vartheta}$, with parameter vector
$\vartheta$, that approximates the function $F$. The neural network
$L_{\theta}$ is modified to generate a vector of length $h$ (instead of a
scalar as in the case of linear PDE) whose entries basically represent spatial
derivative features of the RBF $\phi$. We denote this spatial derivative
feature extractor neural network as $L_{\theta,h}$. The output of
$L_{\theta,h}$ is multiplied with the coefficient vector $\mathbf{c}^{(t)}$ to
get the spatial derivative features of $u$, which is used by $F_{\vartheta}$
to obtain the final output. Figure 1b shows the overall model for the
nonlinear PDE. In this case, for a time-series of measurements (at internal
nodes) $\mathbf{u}^{(t)}\in\mathbb{R}^{N},t=t_{0},t_{0}+\Delta
t,\cdots,t_{0}+K\Delta t$, the neural networks $L_{\theta,h}$ and
$F_{\vartheta}$ are trained jointly with the following objective:

$\displaystyle\min_{\theta,\vartheta,\sigma}\frac{1}{K}\sum_{k=0}^{K-1}\big{\|}$
$\displaystyle\mathbf{u}^{(t_{0}+(k+1)\Delta t)}-\mathbf{u}^{(t_{0}+k\Delta
t)}-\Delta
tF_{\vartheta}\big{(}D_{\theta,h}\Phi^{-1}_{\sigma}\mathbf{u}^{(t_{0}+k\Delta
t)}\big{)}\big{\|}^{2},$ (20)

where $D_{\theta,h}\in\mathbb{R}^{h\times N\times N}$ is a 3D tensor
containing the $h$ spatial derivative matrices of the RBF. Note, the objective
(20) considers nonlinear PDE of form (1) where the temporal order is $p=1$. In
general, for $p\geq 1$, we use the finite difference approximation of
$\frac{\partial^{p}}{\partial t^{p}}$ and the objective is modified as

$\displaystyle\min_{\theta,\vartheta,\sigma}\frac{1}{K}\sum_{k=0}^{K-1}\big{\|}$
$\displaystyle\mathbf{u}^{(t_{0}+(k+1)\Delta
t)}-\textstyle\sum_{q=1}^{p}\binom{p}{q}(-1)^{q+1}\mathbf{u}^{(t_{0}+(k-q+1)\Delta
t)}-\Delta
tF_{\vartheta}\big{(}D_{\theta,h}\Phi^{-1}_{\sigma}\mathbf{u}^{(t_{0}+k\Delta
t)}\big{)}\big{\|}^{2}$ (21)

The training objective of (20) or (21) uses measurements from all sites to
predict the values of $F$ at those sites, but does not incorporate any
constraints to generalize $F_{\vartheta}$ at locations where measurement is
not available. Accuracy of RBF-based methods at any generic location
$\mathbf{x}\in\Omega$ strongly depends on the shape parameter $\sigma$ of the
RBF. Generally the shape parameter is chosen using leave-one-out cross
validation method [28, 29]. Note, as the spatial differential operators are
functions of $\sigma$, the parameters of $L_{\theta,h}$ are dependent on
$\sigma$. We adopt the leave-one-out strategy in our training algorithm. At
each training step, we randomly leave one measurement site and the model uses
the remaining $(N-1)$ sites for prediction. The left out measurement site is
used only for computing the loss function. The overall training procedure is
summarized in Algorithm 1.

1

Input: Training samples {measurement sites $X$,
measurement sequence $\mathbf{u}^{(t)},\ t=t_{0},t_{0}+\Delta
t,\cdots,t_{0}+K\Delta t$ },
a radial function $\phi$

2

3Arbitrarily initialize parameters $\theta,\vartheta,\sigma$

4while  _not converged_ do

5 Randomly select sample sequence {$X,\mathbf{u}^{(t)},t=t_{0},t_{0}+\Delta
t,\cdots,t_{0}+K\Delta t$}

6 Randomly choose $l$ from $\\{1,\cdots,N\\}$

7 Compute RBF matrix $\Phi_{\sigma}\in\mathbb{R}^{n\times(n-1)}$ such that
$[\Phi_{\sigma}]_{ij}=\phi_{\sigma}(\|\mathbf{x}_{i}-\mathbf{x}_{j}\|),\ 1\leq
i,j\leq N,\ j\neq l$

8 Forward pass $L_{\theta,h}$ to get $D_{\theta,h}\in\mathbb{R}^{h\times
n\times(n-1)}$ such that
$[D_{\theta,h}]_{rij}=L_{\theta,h}^{r}\phi(\|\mathbf{x}-\mathbf{x}_{j}\|)|_{\mathbf{x}=\mathbf{x}_{i}},\newline
1\leq r\leq h,\ 1\leq i,j\leq N,\ j\neq l$

9 Forward pass $F_{\vartheta}$ to get
$F_{\vartheta}\big{(}D_{\theta,h}\Phi^{-1}_{\sigma}\mathbf{u}^{(t)}_{-l}\big{)}\in\mathbb{R}^{N},\
t=t_{0},t_{0}+\Delta t,\cdots,t_{0}+K\Delta t$

10 Update $\theta,\vartheta,\sigma$ in the gradient descent direction of
objective (21)

11 end while

Output: Learned parameters $\theta,\vartheta,\sigma$

Algorithm 1 Training Procedure

Once the models $L_{\theta,h}$ and $F_{\vartheta}$ are trained, they can be
used to get the approximated value of $F$ (Figure 1b) at any generic location
$\mathbf{x}\in\Omega$ and to finally predict $u(t,\mathbf{x})$.

###  3.3 Modeling unknown system of coupled nonlinear PDEs

Consider a coupled system of $M$ nonlinear PDEs associated with $M$
measurement variables $u_{1},\cdots,u_{M}$. In order to learn such a system,
we require radial approximation for each of these $M$ measurement variables.
Therefore, given $N$ measurement sites, we consider a RBF coefficient matrix
$C^{(t)}\in\mathbb{R}^{N\times M}$ (at time $t$) where each column corresponds
to the coefficient vector of one measurement variable. Since the neural
network $L_{\theta,h}$ operates only on the spatial locations, its parameters
are shared across all measurement variables. The output of $L_{\theta,h}$,
i.e. the spatial derivative features of the RBF $\phi$, is multiplied with
coefficient matrix $C^{(t)}$ to obtain the spatial derivative features of
length $h$ for each of the measurement variables. These spatial derivative
features are concatenated and finally passed through the neural network
$F_{\vartheta}$ to get the final vector valued ($\in\mathbb{R}^{M}$) output.
Figure 1c shows the overall model for the coupled system of nonlinear PDEs.

##  4 Experimental Evaluation

We evaluate the proposed method for two dimensional wave equation and Burgers-
Fisher equation. The first one is a linear PDE, whereas the latter is a
coupled system of two nonlinear PDEs. Though we show results for two (spatial)
dimensional system, the method is scalable to higher dimensional system due to
RBF’s inherent dimension independent formulation.

###  4.1 Training and testing settings

For both experiments, we train the models with data from discrete measurement
sites on a square geometry $[-1,1]\times[-1,1]$. The measurement sites are
chosen using K-means clustering. For testing, we consider four different
settings:

  1. (i)

same square geometry with same number of measurement sites (though their
locations are different)

  2. (ii)

same square geometry with different number of measurement sites

  3. (iii)

a circular geometry with different boundary condition

  4. (iv)

an annular geometry with different boundary condition

Different geometries with example measurement site configurations are shown in
Figure 2.

Figure 2: (a): Training geometry, (b, c): two additional geometries used for
testing. Red markers show example measurement sites.

###  4.2 Model and optimization hyperparameters

We use Gaussian function
$\phi_{\sigma}(\|\mathbf{x}\|)=e^{-\|\mathbf{x}\|^{2}/2\sigma^{2}}$ as the
radial function, where $\sigma$ is the shape parameter.
For $L_{\theta,h}$ and $F_{\vartheta}$, we use fully connected networks of
layer configurations $(2d+1)-64-32-h$ and $Mh-128-64-32-M$, respectively, with
ReLU activations. $d$ is the dimension of the domain; $d=2$ for both the
experiments. $h$ denotes the length of spatial derivative feature vector; we
use $h=16$. $M$ is the number of PDEs present in the system; $M=1$ for wave
equation and $M=2$ for Burgers-Fisher equation.

We use $300$ sequences of length $200$ for training and validation; another
$50$ sequences are used for testing. We train the models for $200$ epochs in
mini-batches of $32$ using Adam optimizer [30] with an initial learning rate
of $0.001$. During prediction, at each timestep, the linear system
$\Phi\mathbf{c}^{(t)}=\mathbf{u}^{(t)}$ is solved using $\ell_{2}$-regularized
least square with regularization parameter $10^{-4}$ instead of direct
inversion to avoid blowing up coefficients due to noisy estimates. Models are
implemented, trained and tested in PyTorch [31].

###  4.3 Evaluation metrics

We use Signal-to-Noise Ratio (SNR) as metric to evaluate the prediction
quality. Depending on the scenario, we may seek prediction only on the
measurement sites $X$ or on the entire domain $\Omega$. Therefore, we define
two SNR metrics as follows,

$\displaystyle\text{SNR}_{X}(u,\hat{u},t)$
$\displaystyle=10\log_{10}\frac{\sum_{\mathbf{x}\in
X}u^{2}(t,\mathbf{x})}{\sum_{\mathbf{x}\in
X}(u(t,\mathbf{x})-\hat{u}(t,\mathbf{x}))^{2}}$ (22)
$\displaystyle\text{SNR}_{\Omega}(u,\hat{u},t)$
$\displaystyle=10\log_{10}\frac{\sum_{\mathbf{x}\in\Omega}u^{2}(t,\mathbf{x})}{\sum_{\mathbf{x}\in\Omega}(u(t,\mathbf{x})-\hat{u}(t,\mathbf{x}))^{2}}$
(23)

where $\hat{u}(t,\mathbf{x})$ denotes the predicted value of $u$ at time $t$
at location $\mathbf{x}$.

Figure 3: Quantitative analysis for wave equation. (a, b): SNR comparison in
settings (i) with different number of measurement sites; same number of
measurement sites is used in both training and testing. (d, e): SNR comparison
in settings (ii). (c): SNR comparison in settings (iii). (f): SNR comparison
in settings (iv). SNRX and SNRΩ respectively denote SNR only at the
measurement sites $X$ and on the entire domain $\Omega$.

Figure 4: Visualization of the true and predicted states (on entire domain
$\Omega$) of wave equation over time in different settings.

Figure 5: Quantitative analysis for Burgers-Fisher equation. (a, b): SNR
comparison in settings (i) with different number of measurement sites; same
number of measurement sites is used in both training and testing. (d, e): SNR
comparison in settings (ii). (c): SNR comparison in settings (iii). (f): SNR
comparison in settings (iv).

Figure 6: Visualization of the true and predicted states (on entire domain
$\Omega$) of one measurement variable of 2D Burgers-Fisher equation over time
in different settings.

###  4.4 Wave equation

We consider the two-dimensional wave equation given by

$\partial_{t}^{2}u=\mathsf{v}^{2}\nabla^{2}u,\quad\mathbf{x}\in\Omega\subseteq[-1,1]\times[-1,1],\quad
t\in[0,2],$ (24)

where $\mathsf{v}$ is the (constant) wave propagation speed and $\nabla^{2}$
denotes the differential operator Laplacian. Dataset is generated using finite
element method in PDE toolbox of MATLAB. We use $\mathsf{v}=0.1$, and the
following initial condition,

$u(0,\mathbf{x})=ae^{-\varepsilon\|\mathbf{x}-\mathbf{z}\|^{2}},$ (25)

where $a\sim\mathcal{U}(1,2)$, $\varepsilon\sim\mathcal{U}(10,100)$, and
$\mathbf{z}\sim\mathcal{U}(\Omega)$ are chosen randomly for each sequence.
$\mathcal{U}$ denotes the uniform distribution. For training dataset and test
settings (i) and (ii), we use the Dirichlet boundary condition
$u(t,\mathbf{x})=0,\mathbf{x}\in\partial\Omega,t>0$. For test settings (iii)
and (iv), following boundary condition is used,

$u(t,\mathbf{x})=0.2\
\text{sin}\Big{(}\text{tan}^{-1}\frac{x_{2}}{x_{1}}\Big{)},\quad\mathbf{x}=[x_{1},x_{2}]^{T}\in\partial\Omega,\quad
t>0$ (26)

Figure 4 shows the quantitative comparison of proposed method in different
settings in the task of forecasting. Increasing the number of measurement
sites improves SNR (Figure 4(a,b)). Comparing Figures 4a and 4b, it can be
seen that difference between prediction accuracy at measurement sites only
(SNRX) and prediction accuracy on the entire domain (SNRΩ) decreases with
time. Initially, the neural networks compensate for the spatial approximation
error of RBFs to some extent and provide better accuracy at the measurement
sites since they are trained only at those sites. However, as the prediction
quality of the neural networks degrades over time, both the SNRs become very
similar. Figure 4(d,e) show that the prediction quality is maintained even if
we use different number of measurement sites compared to the training
scenario. However, a larger change in number of measurement site would require
retuning the RBF shape parameter. Figure 4(c,f) show that the learned model
can be used in different settings (different geometries with different
boundary conditions). The model is trained in just one geometry with one
boundary condition resulting in accuracy drop when tested in other settings;
generalizability of the learned model can be improved if trained in multiple
settings. Qualitative comparison of predicted maps in different settings along
with ground truth is depicted in Figure 4.

###  4.5 Burgers-Fisher equation

We consider the two-dimensional coupled Burgers-Fisher equations given by

$\displaystyle\partial_{t}u_{1}+[u_{1},u_{2}]^{T}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\nabla
u_{1}$ $\displaystyle=\nu\nabla^{2}u_{1}+\alpha u_{1}(1-u_{1}),$
$\displaystyle\partial_{t}u_{2}+[u_{1},u_{2}]^{T}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\nabla
u_{2}$ $\displaystyle=\nu\nabla^{2}u_{2}+\alpha u_{2}(1-u_{2}),$
$\displaystyle\mathbf{x}\in\Omega\subseteq[-1,1]$
$\displaystyle\times[-1,1],\quad t\in[0,2],$ (27)

where $\nu$ and $\alpha$ are constant non-negative parameters, and $\nabla$
denotes the gradient operator. We use $[\cdot]^{T}$ to denote a column vector,
whereas $[\cdot]$ represents a closed interval. Dataset is generated using
finite element method in PDE toolbox of MATLAB. We use $\nu=0.1,\alpha=1$, and
the following initial condition,

$u_{m}(0,\mathbf{x})=\sum_{|\omega|,|\beta|<4}\lambda_{\omega,\beta}\
\text{cos}([\omega,\beta]^{T}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\mathbf{x})+\gamma_{\omega,\beta}\
\text{sin}([\omega,\beta]^{T}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\mathbf{x}),\quad
m=1,2,$ (28)

where $\lambda_{\omega,\beta},\gamma_{\omega,\beta}\sim\mathcal{N}(0,0.2)$ are
chosen randomly for each sequence. $\mathcal{N}$ denotes the normal
distribution, and
$\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}$
denotes dot product. For training dataset and test settings (i) and (ii), we
use the Dirichlet boundary condition $u_{m}(t,\mathbf{x})=0,m=1,2,\
\mathbf{x}\in\partial\Omega,\ t>0$. For test settings (iii) and (iv), we use
the same boundary condition as (26). Measurement in real environments are
often noisy; we add Gaussian noise $\mathcal{N}(0,0.01\times\text{SD}(u))$ to
the generated data, where $\text{SD}(u)$ is the standard deviation of the
clean training data.

Figure 6 shows the quantitative comparison of proposed method in different
settings in the task of forecasting. In this case, improvement in prediction
accuracy due to changes in number of measurement sites is relatively less
(Figure 6(a,b)). Convergence of the system to a smooth state (as opposed to
wave system) can be attributed to the saturation in RBF approximation with
relatively less number of measurement sites. Qualitative comparison of
predicted maps in different settings along with ground truth is depicted in
Figure 6.

##  5 Conclusion

By integrating RBF collocation method for solving PDEs with deep learning, we
developed a framework for modeling unknown PDEs from sparse observations. RBF
framework enables meshfree computation and makes the method viable for
learning high-dimensional PDEs. Learning the application of differential
operators in time-stepping fashion allows the learned model to be transferable
to different settings of initial and boundary conditions.

Although numerical experiments show promising results, the proposed method has
few limitations that we would like to address in future work. First, even
though the current model can deal with small changes in the number of
measurement sites without retraining, a larger change would require very
different value of the shape parameter. Since the parameters of the neural
networks depends on the shape parameter of the RBF, those also require
retuning to work in a setting where number of measurement sites is very
different from training setting. Developing a neural network model that is
independent of the RBF shape parameter would be an interesting problem to
solve.

Secondly, accuracy of the current method drops if we test in a different
setting involving different geometry and, different boundary and initial
conditions. Pure data-driven models do not perform well in general when
testing scenario is very different from training scenario. Incorporating
knowledge about the system or constraining the learned operators to maintain
physical properties can increase the generalizability of the models.

Another limitation of the current method is that it assumes the physical
parameter (e.g. diffusion coefficient, wave speed) of the underlying system is
uniform. However, in real scenarios, these physical parameters often vary
across space and time. Adapting the model with changes in physical parameters
of the system would be another interesting challenge to address.

## References

  * [1] N. Mai-Duy and T. Tran-Cong, “Numerical solution of differential equations using multiquadric radial basis function networks,” _Neural networks_ , vol. 14, no. 2, pp. 185–199, 2001.
  * [2] I. Sutskever, O. Vinyals, and Q. V. Le, “Sequence to sequence learning with neural networks,” in _Advances in neural information processing systems_ , 2014, pp. 3104–3112.
  * [3] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Advances in neural information processing systems_ , 2017, pp. 5998–6008.
  * [4] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, “Empirical evaluation of gated recurrent neural networks on sequence modeling,” _arXiv preprint arXiv:1412.3555_ , 2014.
  * [5] S. Xingjian, Z. Chen, H. Wang, D.-Y. Yeung, W.-K. Wong, and W.-c. Woo, “Convolutional lstm network: A machine learning approach for precipitation nowcasting,” in _Advances in neural information processing systems_ , 2015, pp. 802–810.
  * [6] C. Finn, I. Goodfellow, and S. Levine, “Unsupervised learning for physical interaction through video prediction,” in _Advances in neural information processing systems_ , 2016, pp. 64–72.
  * [7] N. Lee, W. Choi, P. Vernaza, C. B. Choy, P. H. Torr, and M. Chandraker, “Desire: Distant future prediction in dynamic scenes with interacting agents,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2017, pp. 336–345.
  * [8] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” _Journal of Computational Physics_ , vol. 378, pp. 686–707, 2019.
  * [9] Z. Long, Y. Lu, X. Ma, and B. Dong, “Pde-net: Learning pdes from data,” in _International Conference on Machine Learning_ , 2018, pp. 3208–3216.
  * [10] Y. Long, X. She, and S. Mukhopadhyay, “Hybridnet: Integrating model-based and data-driven learning to predict evolution of dynamical systems,” in _Conference on Robot Learning_ , 2018, pp. 551–560.
  * [11] E. de Bezenac, A. Pajot, and P. Gallinari, “Deep learning for physical processes: Incorporating prior scientific knowledge,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2019, no. 12, p. 124009, 2019\.
  * [12] P. Saha, S. Dash, and S. Mukhopadhyay, “Physics-incorporated convolutional recurrent neural networks for source identification and forecasting of dynamical systems,” _arXiv:2004.06243_ , 2020.
  * [13] J. Sirignano and K. Spiliopoulos, “Dgm: A deep learning algorithm for solving partial differential equations,” _Journal of computational physics_ , vol. 375, pp. 1339–1364, 2018.
  * [14] Y. Khoo, J. Lu, and L. Ying, “Solving parametric pde problems with artificial neural networks,” _European Journal of Applied Mathematics_ , p. 1–15, Jul 2020.
  * [15] N. Geneva and N. Zabaras, “Modeling the dynamics of pde systems with physics-constrained deep auto-regressive networks,” _Journal of Computational Physics_ , vol. 403, p. 109056, 2020.
  * [16] J. Sirignano, J. F. MacArt, and J. B. Freund, “Dpm: A deep learning pde augmentation method with application to large-eddy simulation,” _Journal of Computational Physics_ , vol. 423, p. 109811, 2020.
  * [17] Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner, “Learning data-driven discretizations for partial differential equations,” _Proceedings of the National Academy of Sciences_ , vol. 116, no. 31, pp. 15 344–15 349, 2019\.
  * [18] E. F. Bollig, N. Flyer, and G. Erlebacher, “Solution to pdes using radial basis function finite-differences (rbf-fd) on multiple gpus,” _Journal of Computational Physics_ , vol. 231, no. 21, pp. 7133–7151, 2012.
  * [19] G. E. Fasshauer, _Meshfree approximation methods with MATLAB_. World Scientific, 2007, vol. 6.
  * [20] Z. Wu and J. Liu, “Meshless method for numerical solution of pde using hermitian interpolation with radial basis,” in _Geometric Computation_. World Scientific, 2004, pp. 209–220.
  * [21] M. Raissi, “Deep hidden physics models: Deep learning of nonlinear partial differential equations,” _The Journal of Machine Learning Research_ , vol. 19, no. 1, pp. 932–955, 2018.
  * [22] Z. Long, Y. Lu, and B. Dong, “Pde-net 2.0: Learning pdes from data with a numeric-symbolic hybrid deep network,” _Journal of Computational Physics_ , vol. 399, p. 108925, 2019.
  * [23] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini, “The graph neural network model,” _IEEE Transactions on Neural Networks_ , vol. 20, no. 1, pp. 61–80, 2008.
  * [24] S. Seo and Y. Liu, “Differentiable physics-informed graph networks,” _arXiv preprint arXiv:1902.02950_ , 2019.
  * [25] S. Seo, C. Meng, and Y. Liu, “Physics-aware difference graph networks for sparsely-observed dynamics,” in _International Conference on Learning Representations_ , 2019.
  * [26] M. D. Buhmann, _Radial basis functions: theory and implementations_. Cambridge university press, 2003, vol. 12.
  * [27] E. J. Kansa, “Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics—ii solutions to parabolic, hyperbolic and elliptic partial differential equations,” _Computers & mathematics with applications_, vol. 19, no. 8-9, pp. 147–161, 1990.
  * [28] S. Rippa, “An algorithm for selecting a good value for the parameter c in radial basis function interpolation,” _Advances in Computational Mathematics_ , vol. 11, no. 2-3, pp. 193–210, 1999.
  * [29] G. E. Fasshauer and J. G. Zhang, “On choosing “optimal” shape parameters for rbf approximation,” _Numerical Algorithms_ , vol. 45, no. 1-4, pp. 345–368, 2007.
  * [30] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _3rd International Conference on Learning Representations, ICLR 2015_ , 2015\.
  * [31] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” in _Advances in neural information processing systems_ , 2019, pp. 8026–8037.