Year: 2020

Abstract:
Partial differential equations (PDEs) play a crucial role in studying a vast number of problems in science and engineering. Numerically solving nonlinear and/or highdimensional PDEs is frequently a challenging task. Inspired by the traditional finite difference and finite elements methods and emerging advancements in machine learning, we propose a sequence-to-sequence learning (Seq2Seq) framework called Neural-PDE, which allows one to automatically learn governing rules of any timedependent PDE system from existing data by using a bidirectional LSTM encoder, and predict the solutions in next n time steps. One critical feature of our proposed framework is that the Neural-PDE is able to simultaneously learn and simulate all variables of interest in a PDE system. We test the Neural-PDE by a range of examples, from one-dimensional PDEs to a multi-dimensional and nonlinear complex fluids model. The results show that the Neural-PDE is capable of learning the initial conditions, boundary conditions and differential operators defining the initial-boundary-value problem of a PDE system without the knowledge of the specific form of the PDE system. In our experiments, the Neural-PDE can efficiently extract the dynamics within 20 epochs training and produce accurate predictions. Furthermore, unlike the traditional machine learning approaches for learning PDEs, such as CNN and MLP, which require great quantity of parameters for model precision, the Neural-PDE shares parameters among all time steps, and thus considerably reduces computational complexity and leads to a fast learning algorithm.

heading: 1 INTRODUCTION
text: The research of time-dependent partial differential equations (PDEs) is regarded as one of the most important disciplines in applied mathematics. PDEs appear ubiquitously in a broad spectrum of fields including physics, biology, chemistry, and finance, to name a few. Despite their fundamental importance, most PDEs can not be solved analytically and have to rely on numerical solving methods. Developing efficient and accurate numerical schemes for solving PDEs, therefore, has been an active research area over the past few decades (Courant et al., 1967; Osher & Sethian, 1988; LeVeque; Cockburn et al., 2012; Thomas, 2013; Johnson, 2012). Still, devising stable and accurate schemes with acceptable computational cost is a difficult task, especially when nonlinear and(or) high-dimensional PDEs are considered. Additionally, PDE models emerged from science and engineering disciplines usually require huge empirical data for model calibration and validation, and determining the multidimensional parameters in such a PDE system poses another challenge (Peng et al., 2020).
Deep learning is considered to be the state-of-the-art tool in classification and prediction of nonlinear inputs, such as image, text, and speech (Litjens et al., 2017; Devlin et al., 2018; LeCun et al., 1998; Krizhevsky et al., 2012; Hinton et al., 2012). Recently, considerable efforts have been made to employ deep learning tools in designing data-driven methods for solving PDEs (Han et al., 2018; Long et al., 2018; Sirignano & Spiliopoulos, 2018; Raissi et al., 2019). Most of these approaches are based on fully-connected neural networks (FCNNs), convolutional neural networks(CNNs) and multilayer perceptron (MLP). These neural network structures usually require an increment of the layers to improve the predictive accuracy (Raissi et al., 2019), and subsequently lead to a more complicated model due to the additional parameters. Recurrent neural networks (RNNs) are one type of neural network architectures. RNNs predict the next time step value by using the input data from the current
and previous states and share parameters across all inputs. This idea (Sherstinsky, 2020) of using current and previous step states to calculate the state at the next time step is not unique to RNNs. In fact, it is ubiquitously used in numerical PDEs. Almost all time-stepping numerical methods applied to solve time-dependent PDEs, such as Euler’s, Crank-Nicolson, high-order Taylor and its variance Runge-Kutta (Ascher et al., 1997) time-stepping methods, update numerical solution by utilizing solution from previous steps.
This motivates us to think what would happen if we replace the previous step data in the neural network with numerical solution data to PDE supported on grids. It is possible that the neural network behaves like a time-stepping method, for example, forward Euler’s method yields the numerical solution at a new time point as the current state output (Chen et al., 2018). Since the numerical solution on each of the grid point (for finite difference) or grid cell (for finite element) computed at a set of contiguous time points can be treated as neural network input in the form of one time sequence of data, the deep learning framework can be trained to predict any time-dependent PDEs from the time series data supported on some grids if the bidirectional structure is applied (Huang et al., 2015; Schuster & Paliwal, 1997). In other words, the supervised training process can be regarded as a practice of the deep learning framework to learn the numerical solution from the input data, by learning the coefficients on neural network layers.
Long Short-Term Memory (LSTM) (Hochreiter & Schmidhuber, 1997) is a neural network built upon RNNs. Unlike vanilla RNNs, which suffer from losing long term information and high probability of gradient vanishing or exploding, LSTM has a specifically designed memory cell with a set of new gates such as input gate and forget gate. Equipped with these new gates which control the time to preserve and pass the information, LSTM is capable of learning long term dependencies without the danger of having gradient vanishing or exploding. In the past two decades, LSTM has been widely used in the field of natural language processing (NLP), such as machine translation, dialogue systems, question answering systems (Lipton et al., 2015).
Inspired by numerical PDE schemes and LSTM neural network, we propose a new deep learning framework, denoted as Neural-PDE. It simulates multi-dimensional governing laws, represented by time-dependent PDEs, from time series data generated on some grids and predicts the next n time steps data. The Neural-PDE is capable of intelligently processing related data from all spatial grids by using the bidirectional (Schuster & Paliwal, 1997) neural network, and thus guarantees the accuracy of the numerical solution and the feasibility in learning any time-dependent PDEs. The detailed structures of the Neural-PDE and data normalization are introduced in Section 3.
The rest of the paper is organized as follows. Section 2 briefly reviews finite difference method for solving PDEs. Section 3 contains detailed description of designing the Neural-PDE. In Section 4 and Appendix A of the paper, we apply the Neural-PDE to solve four different PDEs, including the 1-dimensional(1D) wave equation, the 2-dimensional(2D) heat equation, and two systems of PDEs: the invicid Burgers’ equations and a coupled Navier Stokes-Cahn Hilliard equations, which widely appear in multiscale modeling of complex fluid systems. We demonstrate the robustness of the Neural-PDE, which achieves convergence within 20 epochs with an admissible mean squared error, even when we add Gaussian noise in the input data.

heading: 2 PRELIMINARIES

heading: 2.1 TIME DEPENDENT PARTIAL DIFFERENTIAL EQUATIONS
text: A time-dependent partial differential equation is an equation of the form:
ut = f(x1, · · · , u, ∂u ∂x1 , · · · , ∂u ∂xn ,
∂2u
∂x1∂x1 , · · · , ∂
2u
∂x1∂xn , · · · , ∂
nu
∂x1 · · · ∂xn ) , (2.1.1)
where u = u(t, x1, ..., xn) is known, xi ∈ R are spatial variables, and the operator f maps R 7→ R. For example, consider the parabolic heat equation: ut = α2∆u, where u represents the temperature and f is the Laplacian operator ∆. Eq. (2.1.1) can be solved by finite difference methods, which is briefly reviewed below for the self-completeness of the paper.

heading: 2.2 FINITE DIFFERENCE METHOD
text: Consider using a finite difference method (FDM) to solve a two-dimensional second-order PDE of the form: ut = f(x, y, ux, uy, uxx, uyy), (x, y) ∈ Ω ⊂ R2, t ∈ R+ ∪ {0} , (2.2.1) with some proper boundary conditions. Let Ω be Ω = [xa, xb]× [ya, yb], and
uni,j = u(xi, yj , tn) (2.2.2)
where tn = nδt, 0 ≤ n ≤ N , and δt = TN for t ∈ [0, T ], and some large integer N . xi = iδx, 0 ≤ i ≤ Nx, δx = xa−xbNx for x ∈ [xa, xb]. yj = jδy, 0 ≤ j ≤ Ny, δy = ya−yb Ny
for y ∈ [ya, yb]. Nx and Ny are integers.
The central difference method approximates the spatial derivatives as follows (Thomas, 2013):
ux(xi, yj , t) = 1
2δx (ui+1,j − ui−1,j) +O(δx2) , (2.2.3)
uy(xi, yj , t) = 1
2δy (ui,j+1 − ui,j−1) +O(δy2) , (2.2.4)
uxx(xi, yj , t) = 1
δx2 (ui+1,j − 2ui,j + ui−1,j) +O(δx2) , (2.2.5)
uyy(xi, yj , t) = 1
δy2 (ui,j+1 − 2ui,j + ui,j−1) +O(δy2) . (2.2.6)
To this end, the explicit time-stepping scheme to update next step solution un+1 is given by:
uni,j ≈ Un+1i,j = U n i,j + δtf(xi, yj , U n i,j , U n i,j−1, U n i,j+1, U n i+1,j , U n i−1,j) , (2.2.7)
≡ F(xi, yj , δx, δy, δt, Uni,j , Uni,j−1, Uni,j+1, Uni+1,j , Uni−1,j) , (2.2.8)
Apparently, the finite difference method (2.2.7) for updating un+1 on a grid point relies on the previous time steps’ solutions, supported on the grid point and its neighbours. The scheme (2.2.7) updates un+1i,j using four points of un values (see Figure 1). Similarly, the finite element method (FEM) approximates the new solution by calculating the corresponded mesh cell coefficient (see Appendix), which is updated by its related nearby coefficients on the mesh. From this perspective, one may regard the numerical schemes for solving time-dependent PDEs as methods catching the information from neighbourhood data of interest.

heading: 3 PROPOSED METHOD

heading: 3.1 MATHEMATICAL MOTIVATION
text: Recurrent neural network including LSTM is an artificial neural network structure of the form (Lipton et al., 2015):
ht = σ(Whxxt + Whhht−1 + bh) ≡ σa(xt,ht−1) ≡ σb(x0,x1,x2, · · · ,xt) , (3.1.1) where xt ∈ Rd is the input data of the tth state and ht−1 ∈ Rh denotes the processed value in its previous state by the hidden layers. The output yt of the current state is updated by the current state value ht:
yt = σ(Whyht + by) (3.1.2)
≡ σc(ht) ≡ σd(x0,x1,x2, · · · ,xt) . (3.1.3)
Here Whx ∈ Rh×d, Whh ∈ Rh×h, Why ∈ Rh×h are the matrix of weights, vectors bh, by ∈ Rh are the coefficients of bias, and σ, σa, σb, σc, σd are corresponded activation and mapping functions.
With proper design of input and forget gate, LSTM can effectively yield a better control over the gradient flow and better preserve useful information from long-range dependencies (Graves & Schmidhuber, 2005).
Now consider a temporally continuous vector function u ∈ Rn given by an ordinary differential equation with the form:
du(t)
dt = g(u(t)) . (3.1.4)
Let un = u(t = nδt), a forward Euler’s method for solving u can be easily derived from the Taylor’s theorem which gives the following first-order accurate approximation of the time derivative:
dun
dt =
un+1 − un
δt +O(δt) . (3.1.5)
Then we have: du
dt = g(u)
(3.1.5)−−−−→ un+1 = un + δt g(un) +O(δt2)
→ ûn+1 = f1(ûn) = f1 ◦ f1 ◦ · · · f1(û0)︸ ︷︷ ︸ n
(3.1.6)
Here ûn ≈ u(nδt) is the numerical approximation and f1 ≡ un+δt g(un) : Rn → Rn. Combining equations (3.1.1) and (3.1.6) one may notice that the residual networks, recurrent neural network and also LSTM networks can be regarded as a numerical scheme for solving time-dependent differential equations if more layers are added and smaller time steps are taken. (Chen et al., 2018)
Canonical structure for such recurrent neural network usually calculate the current state value by its previous time step value ht−1 and current state input xt. Similarly, in numerical PDEs, the next step data at a grid point is updated from the previous (and current) values on its nearby grid points (see Eq. 2.2.7).
Thus, what if we replace the temporal input ht−1 and xt with spatial information? A simple sketch of the upwinding method for a 1d example of u(x, t):
ut + νux = 0 (3.1.7) will be:
un+1i = u n i − ν
δt δx (uni − uni−1) +O(δx, δt)→ ûn+1i = f2(û n i−1, û n i ) (3.1.8)
≡ fθ ( fη(xi,hi−1(u)) ) = fθ,η ( ûn0 , û n 1 , · · · , ûni−1, ûni ) = vn+1i (3.1.9)
xi = û n i , hi−1(û) = σ(û n i−1,hi−2(û)) ≡ fη(ûn0 , ûn1 , ûn2 , · · · , ûni−1). (3.1.10)
Here let vn+1i be the prediction of û n+1 i processed by neural network. We replace the temporal previous state ht−1with spacial grid value hi−1 and input the numerical solution ûni ≈ u(iδx, nδt) as current state value, which indicates the neural network could be seen as a forward Euler method for equation 3.1.7 (Lu et al., 2018). Function f2 ≡ ûni − ν δtδx (û n i − ûni−1) : R→ R and the function fθ represents the dynamics of the hidden layers in decoder with parameters θ, and fη specifies the dynamics of the LSTM layer (Hochreiter & Schmidhuber, 1997; Graves & Schmidhuber, 2005) in encoder withe parameters η. The function fθ,η simulates the dynamics of the Neural-PDE with paramaters θ and η. By applying Bidirectional neural network, all grid data are transferred and it enables LSTM to simulate the PDEs as :
vn+1i = fθ ( fη(hi+1(ˆ̂u), û n i ,hi−1(û)) ) (3.1.11)
hi+1(û) ≡ fη(ûni+1, ûni+2, ûni+3, · · · , ûnk ). (3.1.12) For a time-dependent PDE, if we map all our grid data into an input matrix which contains the information of δx, δt, then the neural network would regress such coefficients as constants and will learn and filter the physical rules from all the k mesh grids data as:
vn+1i = fθ,η ( ûn0 , û n 1 , û n 2 , · · · , ûnk ) (3.1.13)
The LSTM neural network is designed to overcome the vanishing gradient issue through hidden layers, therefore we use such recurrent structure to increase the stability of the numerical approach in deep learning. The highly nonlinear function fθ,η simulates the dynamics of updating rules for un+1i , which works in a way similar to a finite difference method (section 2.2) or a finite element method.
3.2 NEURAL-PDE
In particular, we use the bidirectional LSTM (Hochreiter & Schmidhuber, 1997; Graves & Schmidhuber, 2005) to better retain the state information from data on grid points which are neighbourhoods in the mesh but far away in input matrix.
The right frame of Figure 3 shows the overall design of the Neural-PDE. Denote the time series data at collocation points as aN1 ,a N 2 , · · · ,aNk with aNi = [û0i , û 1 i , · · · , ûNi ] at ith point. The superscript represents different time points. The Neural-PDE takes the past states {aN1 ,aN2 , · · · ,aNk } of all collocation points, and outputs the predicted future states {bM1 , bM2 , · · · , bMk }, where bMi = [v N+1 i , v N+2 i , · · · , v N+M i ] is the Neural-PDE prediction for the ith collocation point at time points from N + 1 to N +M . The data from time point 0 to N are the training data set.
The Neural-PDE is an encoder-decoder style sequence model that first maps the input data to a low dimensional latent space that
hi = −−−−→ LSTM(ai)⊕ ←−−−− LSTM(ai), (3.2.1)
where ⊕ denotes concatenation and hi is the latent embedding of point ai under the environment.
One then decoder, another bi-lstm with a dense layer: vi = (−−−−→ LSTM(hi)⊕ ←−−−− LSTM(hi) ) ·W, (3.2.2)
where W is the learnable weight matrix in the dense layer.
During training process, mean squared error (MSE) loss L is used as we typically don’t know the specific form of the PDE.
L = N+M∑ t=N+1 k∑ i=1 ||ûti − vti ||2 , (3.2.3)

heading: 3.3 DATA INITIALIZATION AND GRID POINT RESHAPE
text: In order to feed data into our sequence model framework, we map the PDE solution data onto a K × N matrix, where K ∈ Z+ is the dimension of the grid points and N ∈ Z+ is the length of the time series data on each grid point. There is no regularization for the input order of the grid points data in the matrix because of the bi-directional structure of the Neural-PDE. For example, a 2d
heat equation at some time t is reshaped into a 1d vector (See Fig. 2). Then the matrix is formed accordingly.
For a n-dimensional time-dependent partial differential equation with K collocation points, the input and output data for t ∈ (0, T ) will be of the form:
A(K,N) =  aN0 ... aN`
... aNK
 =  û00 û 1 0 · · · ûn0 · · · ûN0 ... ... . . . ... . . . ... û0` û 1 ` · · · ûn` · · · ûN` ... ... . . . ... . . . ...
û0K û 1 K · · · ûnK · · · ûNK
 (3.3.1)
B(K,M) =  bM0 ... bM`
... bMK
 =  vN+10 v N+2 0 · · · v N+m 0 · · · v N+M 0 ... ... . . . ... . . . ... vN+1` v N+2 ` · · · v N+m ` · · · v N+M k ... ... . . . ... . . . ...
vN+1K v N+2 K · · · v N+m K · · · v N+M K
 (3.3.2)
Here N = Tδt and each row ` represents the time series data at the ` th mesh grid, and M is the time length of the predicted data.
By adding Bidirectional LSTM encoder in the Neural-PDE, it will automatically extract the information from the time series data as:
B(K,M) = PDESolver(A(K,N)) = PDESolver(aN0 ,a N 1 , · · · aNi , · · · ,aNK) (3.3.3)

heading: 4 COMPUTER EXPERIMENTS
text: Since the Neural-PDE is a sequence to sequence learning framework which allows one to predict within any time period by the given data. One may test the Neural-PDE using different permutations of training and predicting time periods for its efficiency, robustness and accuracy. In the following examples, the whole dataset is randomly splitted in 80% for traning and 20% for testing. We will predict the next tp ∈ [31× δt, 40× δt] PDE solution by using its previous ttr ∈ [0, 30× δt] data as:
B(K, 10) = PDESolver(A(K, 30)) (4.0.1)
Table 1 summaries the experimental results of the Neural-PDE model on 4 different PDEs, which achieve extremely small MSEs from ∼ 10−5 to ∼ 10−7. Table 2 shows the comparison results of our proposed Neural-PDE with the state-of-the-art method Physically Informed Artificial Neural Networks (PINN) (Raissi et al., 2019) on two PDEs (1d Allen-Cahn and 1d Burgers’ equation). Neural-PDE is able to outperform PINN while having much less parameters, where PINN contains 4 hidden layers with 200 neurons per layer and Neural-PDE only consists of 3 layers (2 bi-lstm with 20 neurons per layer and 1 dense output layer with 10 neurons).
EXAMPLE: INVISCID BURGERS’ EQUATION
Inviscid Burgers’ equation is a classical nonlinear PDE in fluid dynamics. In this example, we consider an invicid Burgers’ equation which has the following hyperbolic form:
∂u ∂t + u ∂u ∂x + v ∂u ∂y = 0, ∂v ∂t + u ∂u ∂x + v ∂u ∂y = 0 (4.0.2)
Ω = [0, 1]× [0, 1], t ∈ [0, 1], (4.0.3)
and with initial and boundary conditions:
u(0.25 ≤ x ≤ 0.75, 0.25 ≤ y ≤ 0.75, t = 0) = 0.9 (4.0.4) v(0.25 ≤ x ≤ 0.75, 0.25 ≤ y ≤ 0.75, t = 0) = 0.5 (4.0.5) u(0, y, t) = u(1, y, t) = v(x, 0, t) = v(x, 1, t) = 0 (4.0.6)
The invicid Burgers’ equation is hard to deal with in numerical PDEs due to the discontinuities (shock waves) in the solutions. We use a upwinding finite difference scheme to create the training data and put the velocity u, v in to the input matrix. Let δx = δy = 10−2, δt = 10−3, our empirical results (see Figure 4) show that the Neural-PDE is able to learn the shock waves, boundary conditions and the rules of the equation, and predict u and v simultaneously with an overall MSE of 1.4018× 10−5. The heat maps of exact solution and predicted solution are shown in Figure 5.

heading: EXAMPLE: MULTISCALE MODELING: COUPLED CAHN–HILLIARD–NAVIER–STOKES SYSTEM
text: Finally, let’s consider the following 2d Cahn–Hilliard–Navier–Stokes system widely used for modeling complex fluids:
ut + u · ∇u = −∇p+ ν∆u− φ∇µ , (4.0.7) φt +∇ · (uφ) = M∆µ , (4.0.8)
µ = λ(−∆φ+ φ η2 (φ2 − 1)) (4.0.9)
∇ · u = 0 (4.0.10)
In this complicated example we will use the following initial condition:
φ(x, y, 0) = ( 1
2 − 50 tanh(f1 − 0.1)) + (
1 2 − 50 tanh(f2 − 0.1)), I.C. (4.0.11)
f1 = √ (x+ 0.12)2 + (y)2, f2 = √
(x− 0.12)2 + (y)2 (4.0.12) with x ∈ [−0.5, 0.5], y ∈ [−0.5, 0.5], t ∈ [0, 1], M = 0.1, ν = 0.01 (4.0.13)
This fluid system can be derived by the energetic variational approach (Forster, 2013). The complex fluids system has the following features: the micro-structures such as the molecular configurations, the interaction between different scales and the competition between multi-phase fluids (Hyon et al., 2010). Here u is the velocity and φ(x, y, t) ∈ [0, 1] denotes the volume fraction of one fluid phase. M is the diffusion coefficient and µ is the chemical potential of φ. Equation (4.0.10) indicates the incompressibility of the fluid. Solving such PDE system is notorious because of its high nonlinearity and multi-physical and coupled features. One may use the decoupled projection method (Guermond et al., 2006) to numerically solve it efficiently or an implicit method which however is computationally
expensive. Another challenge of deep learning in solving a system like this is how to process the data to improve the learning efficiency when the input matrix consists of variables such as φ ∈ [0, 1] with large magnitude value and variable of very small values such as p ∼ 10−5. For the Neural-PDE to better extract and learn the physical features of variables in different scales, we normalized the p data with a sigmoid function. Set δt = 5× 10−4, and here the training dataset is generated by a FEM solver FreeFem++ (Hecht, 2012) using a Crank-Nicolson finite element scheme. Our Neural-PDE prediction shows that the physical features of p and φ have been successfully captured with an overall MSE: 6.1631× 10−7 (see Figure 7).

heading: 5 CONCLUSIONS
text: In this paper, we proposed a novel sequence recurrent deep learning framework: Neural-PDE, which is capable of intelligently filtering and learning solutions of time-dependent PDEs. One key innovation of our method is that the time marching method from the numerical PDEs is applied in the deep learning framework, and the neural network is trained to explore the accurate numerical solutions for prediction.
The state-of-the-art researches have shown the promising power of deep learning in solving highdimensional nonlinear problems in engineering, biology and finance with efficiency in computation and accuracy in prediction. However, there are still unresolved issues in applying deep learning in PDEs. For instance, the stability and convergence of the numerical algorithms have been rigorously studied by applied mathematicians. Due to the high nonlinearity of the neural network system and the curse of dimensionality (Hutzenthaler et al., 2019), theorems guiding stability and convergence of solutions predicted by the neural network are to be revealed.
Lastly, it would be helpful and interesting if one can theoretically characterize a numerical scheme from the neural network coefficients and learn the forms or mechanics from the scheme and prediction. We leave these questions for further study.

heading: A APPENDIX

heading: A.1 FINITE ELEMENT METHOD
text: Finite element method (FEM) is a powerful numerical method in solving PDEs. Consider a 1D wave equation of u(x, t):
utt − v2uxx = f, x ∈ [a, b] ≡ Ω ⊂ R, t ∈ R+ ∪ {0} , (A.1.1) ux(a, t) = ux(b, t) = 0 . (A.1.2)
The function u is approximated by a function uh :
u(x, t) ≈ uh(x, t) = N∑ i=1 ai(t)ψi(x) (A.1.3)
(A.1.4) where ψi ∈ V , is the basis functions of some FEM space V , and ani denotes the coefficients. N denotes the degrees of freedom. Multiply the equation with an arbitrary test function ψj and integral over the whole domain we have:∫
Ω
uttψj dx+ v 2 ∫ Ω ∇u∇ψj dx = ∫ Ω fψj dx (A.1.5)
(A.1.6) and approximate u(x, t) by uh:
N∑ i ∂2ai(t) ∂t2 ∫ Ω
ψiψj dx︸ ︷︷ ︸ Mi,j
+v2 N∑ i ai(t) ∫ Ω
∇ψi∇ψj dx︸ ︷︷ ︸ Ai,j =
∫ Ω
fψj︸ ︷︷ ︸ b dx , (A.1.7)
≡MTatt + v2ATa = b . (A.1.8) Here M is the mass matrix and A is the stiffness matrix, a = (a1, .., aN )t is a N × 1 vector of the coefficients at time t. The central difference method for time discretization indicates that (Johnson, 2012):
an+1 = 2an − an−1 + M−1(b− v2ATan) , (A.1.9) un+1 ≈ un+1h = N∑ i an+1i ψi(x) . (A.1.10)

heading: A.2 LONG SHORT-TERM MEMORY
text: Long Short-Term Memory Networks (LSTM) (Hochreiter & Schmidhuber, 1997; Graves & Schmidhuber, 2005) are a class of artificial recurrent neural network (RNN) architecture that is commonly used for processing sequence data and can overcome the gradient vanishing issue in RNN. Similar to most RNNs (Mikolov et al., 2011), LSTM takes a sequence {x1,x2, · · · ,xt} as input and learns hidden vectors {h1,h2, · · · ,ht} for each corresponding input. In order to better retain long distance information, LSTM cells are specifically designed to update the hidden vectors. The computation process of the forward pass for each LSTM cell is defined as follows:
it = σ(W (x) i xt + W (h) i ht−1 + W (c) i ct−1 + bi) ,
ft = σ(W (x) f xt + W (h) f ht−1 + W (c) f ct−1 + bf ) ,
ct = ftct−1 + it tanh(W (x) c xt + W (h) c ht−1 + bc) ,
ot = σ(W (x) o xt + W (h) o ht−1 + W (c) o ct + bo),
ht = ot tanh(ct) ,
where σ is the logistic sigmoid function, Ws are weight matrices, bs are bias vectors, and subscripts i, f , o and c denote the input gate, forget gate, output gate and cell vectors respectively, all of which have the same size as hidden vector h.
This LSTM structure is used in the paper to simulate the numerical solutions of partial differential equations.

heading: A.3 EXAMPLES

heading: A.3.1 WAVE EQUATION
text: Consider the 1d wave equation:
utt = cuxx, x ∈ [0, 1], t ∈ [0, 2] , (A.3.1) u(x, 0) = sin(4πx) (A.3.2) u(0, t) = u(1, t) (A.3.3)
Let c = 116π2 and use the analytical solution given by the characteristics for the training and testing data:
u(x, t) = 1
2 (sin(4πx+ t) + sin(4πx− t)) (A.3.4)

heading: A.3.2 HEAT EQUATION
text: The heat equation describes how the motion or diffusion of a heat flow evolves over time. The Black–Scholes model (Black & Scholes, 1973) is also developed based on the physical laws behind the heat equation. Rather than the 1D case that maps the data into a matrix (??) with its original spatial locations, the high dimensional PDEs grids are mapped into matrix without regularization of the position, and the experimental results show that Neural-PDE is able to capture the valuable features regardless of the order of the mesh grids in the matrix. Let’s start with a 2D heat equation as follows:
ut = uxx + uyy, (A.3.5)
u(x, y, 0) = { 0.9, if (x− 1)2 + (y − 1)2 < 0.25 0.1, otherwise (A.3.6)
Ω = [0, 2]× [0, 2], t ∈ [0, 0.15] (A.3.7)

Review summaries:
This work proposes a sequence-to-sequence approach for learning the time evolution of PDEs. The method employs a bi-directional LSTM to predict solutions of a PDE-based formulation for a chosen number of time steps. By itself this is an interesting, and important goal, but the method does not seem to contain any novel components apart from demonstrating that LSTMs can be used to learn data from PDEs. The paper only compares to a simple form of PINNs, but not to a variety of other time forecasting algorithms available in the deep learning field (LSTM are just one of many methods used these days, a more state of the art one being e.g. transformers). In addition, the examples only contain single cases with relatively simple model equations.