PDE-Net 2.0: Learning PDEs from Data with A Numeric-Symbolic Hybrid Deep Network

Zichao Long, Yiping Lu, Bin Dong

Introduction

Differential equations, especially partial differential equations(PDEs), play a prominent role in many disciplines to describe the governing physical laws underlying a given system of interest. Traditionally, PDEs are derived mathematically or physically based on some basic principles, e.g. from Schrödinger’s equations in quantum mechanics to molecular dynamic models, from Boltzmann equations to Navier-Stokes equations, etc. However, the mechanisms behind many complex systems in modern applications (such as many problems in multiphase flow, neuroscience, finance, biological science, etc.) are still generally unclear, and the governing equations of these systems are commonly obtained by empirical formulas . With the recent rapid development of sensors, computational power, and data storage in the last decade, huge quantities of data can now be easily collected, stored and processed. Such vast quantity of data offers new opportunities for data-driven discovery of (potentially new) physical laws. Then, one may ask the following interesting and intriguing question: can we learn a PDE model to approximate the observed complex dynamic data?

Recent work greatly advanced the progress of PDE identification from observed data. However, SINDy requires to build a sufficiently large dictionary which may lead to high memory load and computation cost, especially when the number of model variables is large. Furthermore, the existing methods based on SINDy treat spatial and temporal information of the data separately and does not take full advantage of the temporal dependence of the PDE model. Although the framework presented by is able to learn hidden physical laws using less data than the approach based on SINDy, the explicit form of the PDEs is assumed to be known except for a few scalar learnable parameters. The approach of is specifically designed for advection-diffusion equations, and cannot be readily extended to other types of equations. Therefore, extracting governing equations from data in a less restrictive setting remains a great challenge.

The main objective of this paper is to design a transparent deep neural network to uncover hidden PDE models from observed complex dynamic data with minor prior knowledge on the mechanisms of the dynamics, and to perform accurate predictions at the same time. The reason we emphasize on both model recovery and prediction is because: 1) the ability to conduct accurate long-term prediction is an important indicator of accuracy of the learned PDE model (the more accurate is the prediction, the more confident we have on the underlying recovered PDE model); 2) the trained neural network can be readily used in applications and does not need to be re-trained when initial conditions are altered. Our inspiration comes from the latest development of deep learning techniques in computer vision. An interesting fact is that some popular networks in computer vision, such as ResNet, have close relationship with ODEs/PDEs and can be naturally merged with traditional computational mathematics in various tasks . However, existing deep networks designed in deep learning mostly emphasis on expressive power and prediction accuracy. These networks are not transparent enough to be able to reveal the underlying PDE models, although they may perfectly fit the observed data and perform accurate predictions. Therefore, we need to carefully design the network by combining knowledge from deep learning and numerical PDEs.

2 Our Approach

The proposed deep neural network is an upgraded version of our original PDE-Net . The main difference is the use of a symbolic network to approximate the nonlinear response function, which significantly relaxes the requirement on the prior knowledge on the PDEs to be recovered. During training, we no longer need to assume the general type of the PDE (e.g. convection, diffusion, etc.) is known. Furthermore, due to the lack of prior knowledge on the general type of the unknown PDE models, more carefully designed constraints on the convolution filters as well as the parameters of the symbolic network are introduced. We refer to this upgraded network as PDE-Net 2.0.

Assume that the PDE to be recovered takes the following generic form

PDE-Net 2.0 is designed as a feed-forward network by discretizing the above PDE using forward Euler in time and finite difference in space. The forward Euler approximation of temporal derivative makes PDE-Net 2.0 ResNet-like , and the finite difference is realized by convolutions with trainable kernels (or filters). The nonlinear response function FF is approximated by a symbolic neural network, which shall be referred to as SymNetSymNet. All the parameters of the SymNetSymNet and the convolution kernels are jointly learned from data. To grant full transparency to the PDE-Net 2.0, proper constraints are enforced on the SymNetSymNet and the filters. Full details on the architecture and constraints will be presented in Section 2.

3 Relation with Model Reduction

Data-driven discovery of hidden physical laws and model reduction have a lot in common. Both of them concern on representing observed data using relatively simple models. The main difference is that, model reduction emphasis more on numerical precision rather than acquiring the analytic form of the model.

It is common practice in model reduction to use a function approximator to express the unknown terms in the reduced models, such as approximating subgrid stress for large-eddy simulation or approximating interatomic forces for coarse-grained molecular dynamic systems. Our work may serve as an alternative approach to model reduction and help with analyzing the reduced models.

4 Novelty

The particular novelties of our approach are that we impose appropriate constraints on the learnable filters and use a properly designed symbolic neural network to approximate the response function FF. Using learnable filters makes the PDE-Net 2.0 more flexible, and enables more powerful approximation of unknown dynamics and longer time prediction (see numerical experiments in Section 3 and Section 4). Furthermore, the constraints on the learnable filters and the use of a deep symbolic neural network enable us to uncover the analytic form of FF with minor prior knowledge on the dynamic, which is the main advantage of PDE-Net 2.0 over the original PDE-Net. In addition, the composite representation by the symbolic network is more efficient and flexible than SINDy. Therefore, the proposed PDE-Net 2.0 is distinct from the existing learning based methods to discover PDEs from data.

PDE-Net 2.0: Architecture, Constraints and Training

Inspired by the dynamic system perspective of deep neural networks , we consider forward Euler as the temporal discretization of the evolution PDE (1), and unroll the discrete dynamics to a feed-forward network. One may consider more sophisticated temporal discretization which naturally leads to different network architectures . For simplicity, we focus on forward Euler in this paper.

Here, the operators DijD_{ij} are convolution operators with the underlying filters denoted by qijq_{ij}, i.e. Diju=qij⊛uD_{ij}u=q_{ij}\circledast u. The operators D10D_{10}, D01D_{01}, D11D_{11}, etc. approximate differential operators, i.e. Diju≈∂i+ju∂ix∂jyD_{ij}u\approx\frac{\partial^{i+j}u}{\partial^{i}x\partial^{j}y}. In particular, the operators D00D_{00} is a certain averaging operator. The purpose of introducing the average operators in stead of simply using the identity is to improve the expressive power of the network and enables it to capture more complex dynamics.

Other than the assumption that the observed dynamics is governed by a PDE of the form (1), we assume that the highest order of the PDE is less than some positive integer. Then, we can assume that FF is a function of mm variables with known mm. The task of approximating FF in (1) is equivalent to a multivariate regression problem. In order to be able to identify the analytic form of FF, we use a symbolic neural network denote by SymNetmkSymNet_{m}^{k} to approximate FF, where kk denotes the depth of the network. Note that, if FF is a vector function, we use multiple SymNetmkSymNet_{m}^{k} to approximate the components of FF separately.

Combining the aforementioned approximation of differential operators and the nonlinear response function, we obtain an approximation framework (2) which will be referred to as a δt\delta t-block (see Figure 1). Details of these two components can be found later in Section 2.2 and Section 2.3.

1.2 Multiple δ​t𝛿𝑡\delta t-Blocks:

One δt\delta t-block only guarantees the accuracy of one-step dynamics, which does not take error accumulation into consideration. In order to facilitate a long-term prediction, we stack multiple δt\delta t-blocks into a deep network, and call this network the PDE-Net 2.0 (see Figure 2). The importance of stacking multiple δt\delta t-blocks will be demonstrated by our numerical experiments in Section 3 and 4.

2 Convolutions and Differentiations

In the original PDE-Net , the learnable filters are properly constrained so that we can easily identify their correspondence to differential operators. PDE-Net 2.0 adopts the same constrains as the original version of the PDE-Net. For completeness, we shall review the related notions and concepts, and provide more details.

A profound relationship between convolutions and differentiations was presented by , where the authors discussed the connection between the order of sum rules of filters and the orders of differential operators. Note that the definition of convolution we use follows the convention of deep learning, which is defined as

This is essentially correlation instead of convolution in the mathematics convention. Note that if ff is of finite size, we use periodic boundary condition.

The order of sum rules is closely related to the order of vanishing moments in wavelet theory . We first recall the definition of the order of sum rules.

In practical implementation, the filters are normally finite and can be understood as matrices. For an N×NN\times N filter qq (NN being an odd number), assuming the indices of qq start from −N−12-\frac{N-1}{2}, (3) can be written in the following simpler form

The following proposition from links the orders of sum rules with orders of differential operators.

If, in addition, qq has total sum rules of order K\{∣α∣+1}K\backslash\{|\alpha|+1\} for some K>∣α∣K>|\alpha|, then

According to Proposition 2.1, an α\alphath order differential operator can be approximated by the convolution of a filter with α\alpha order of sum rules. Furthermore, according to (5), one can obtain a high order approximation of a given differential operator if the corresponding filter has an order of total sum rules with K>∣α∣+k,k⩾1K>|\alpha|+k,k\geqslant 1. For example, consider filter

It has a sum rules of order (1,0)(1,0), and a total sum rules of order 3\{2}3\backslash\{2\}. Thus, up to a constant and a proper scaling, qq corresponds to a discretization of ∂∂x\frac{\partial}{\partial x} with second order accuracy.

For an N×NN\times N filter qq, define the moment matrix of qq as

From (8) we can see that filter qq can be designed to approximate any differential operator with prescribed order of accuracy by imposing constraints on M(q)M(q).

For example, if we want to approximate ∂u∂x\frac{\partial u}{\partial x} (up to a constant) by convolution q⊛uq\circledast u where qq is a 3×33\times 3 filter, we can consider the following constrains on M(q)M(q):

Here, ⋆\star means no constraint on the corresponding entry. The constraints described by the moment matrix on the left of (9) guarantee the approximation accuracy is at least first order, and the one on the right guarantees an approximation of at least second order. In particular, when all entries of M(q)M(q) are constrained, e.g.

the corresponding filter can be uniquely determined. In the PDE-Net 2.0, all filters are learned subjected to partial constraints on their associated moment matrices, with at least second order accuracy.

3 Design of S​y​m​N​e​t𝑆𝑦𝑚𝑁𝑒𝑡SymNet: a Symbolic Neural Network

The symbolic neural network SymNetmkSymNet_{m}^{k} of the PDE-Net 2.0 is introduced to approximate the multivariate nonlinear response function FF of (1). Neural networks have recently been proven effective in approximating multivariate functions in various scenarios . For the problem we have in hand, we not only require the network to have good expressive power, but also good transparency so that the analytic form of FF can be readily inferred after training. Our design of SymNetmkSymNet_{m}^{k} is motivated by EQL/EQL÷ proposed by .

The SymNetmkSymNet_{m}^{k}, as illustrated in Figure 3, is a network that takes an mm dimensional vector as input and has kk hidden layers.

Figure 3 shows the symbolic neural network with two hidden layers, i.e. SymNetm2SymNet_{m}^{2}, where ff is a dyadic operation unit, e.g. multiplication or division. In this paper, we focus on multiplication, i.e. we take f(a,b)=a×bf(a,b)=a\times b. Different from EQL/EQL÷, each hidden layer of the SymNetmkSymNet_{m}^{k} directly takes the outputs from the preceding layer as inputs, rather than a linear combination of them. Furthermore, it adds one additional variable (i.e. f(⋅,⋅)f(\cdot,\cdot)) at each hidden layer. To better understand SymNetmkSymNet_{m}^{k}, we present an example in Algorithm 1 showing how SymNet62SymNet_{6}^{2} is constructed. In particular, when bi=0,i=1,2,3b^{i}=0,i=1,2,3,

and W3=(0,0,0,0,0,0,−1,−1)W^{3}=(0,0,0,0,0,0,-1,-1), then SymNet62(u,ux,uy,v,vx,vy)=−uux−vuySymNet_{6}^{2}(u,u_{x},u_{y},v,v_{x},v_{y})=-uu_{x}-vu_{y} which is the right-hand-side of ∂u∂t\frac{\partial u}{\partial t} of the Burgers’ equation without viscosity.

The SymNetmkSymNet_{m}^{k} can represent all polynomials of variables (x1,x2,…,xm)(x_{1},x_{2},\ldots,x_{m}) with the total number of multiplications not exceeding kk. If needed, one can add more operations to the SymNetmkSymNet_{m}^{k} to increase the capacity of the network.

Now, we show that SymNetmkSymNet_{m}^{k} is more compact than the dictionaries of SINDy. For that, we first introduce some notions.

Define the set of all polynomials of mm variables (x1,⋯ ,xm)(x_{1},\cdots,x_{m}) with the total number of multiplications not exceeding kk as Pk[x1,⋯ ,xm]\mathcal{P}^{k}[x_{1},\cdots,x_{m}]. Here, the total number of multiplications of Pk[x1,⋯ ,xm]\mathcal{P}^{k}[x_{1},\cdots,x_{m}] is counted as follows:

For any monomial of degree kk, if k≥2k\geq 2, then the number of multiplications of the monomial is counted as k−1k-1. When k=1k=1 or , the count is 0.

For any polynomial PP, the total number of multiplications is counted as the sum of the number of multiplications of its monomials.

For example, ∑i=1mxi+∑i=1kxixi+1\sum_{i=1}^{m}x_{i}+\sum_{i=1}^{k}x_{i}x_{i+1} and ∏i=1k+1xi\prod_{i=1}^{k+1}x_{i} with k<mk<m are all members of Pk[x1,⋯ ,xm]\mathcal{P}^{k}[x_{1},\cdots,x_{m}]. The elements in Pk[x1,⋯ ,xm]\mathcal{P}^{k}[x_{1},\cdots,x_{m}] are of simple forms when kk is relatively small. The following proposition shows that SymNetmkSymNet_{m}^{k} can represent all polynomials of variables (x1,x2,…,xm)(x_{1},x_{2},\ldots,x_{m}) with the total number of multiplications not exceeding kk. Note that the actual capacity of SymNetmkSymNet_{m}^{k} is larger than Pk[x1,⋯ ,xm]\mathcal{P}^{k}[x_{1},\cdots,x_{m}], i.e. Pk[x1,⋯ ,xm]\mathcal{P}^{k}[x_{1},\cdots,x_{m}] is a subset of the set of functions that SymNetmkSymNet_{m}^{k} can represent.

For any P∈Pk[x1,⋯ ,xm]P\in\mathcal{P}^{k}[x_{1},\cdots,x_{m}], there exists a set of parameters for SymNetmkSymNet_{m}^{k} such that

Proof: We prove this proposition by induction. When k=1k=1, the conclusion obviously holds. Suppose the conclusion holds for kk. For any polynomial P∈Pk+1[x1,⋯ ,xm]P\in\mathcal{P}^{k+1}[x_{1},\cdots,x_{m}], we only need to consider the cases when PP has a total number of multiplications greater than 1.

We take any monomial of PP that has degree greater than 1, which we suppose take the form x1x2⋅Ax_{1}x_{2}\cdot A where AA is a monomial of variable (x1,…,xm)(x_{1},\ldots,x_{m}). Then, PP can be written as P=x1x2A+QP=x_{1}x_{2}A+Q. Define new variable xm+1=x1x2x_{m+1}=x_{1}x_{2}. Then, we have

By the induction hypothesis, there exists a set of parameters such that P=SymNetm+1k(x1,⋯ ,xm+1)P=SymNet_{m+1}^{k}(x_{1},\cdots,x_{m+1}).

We take the linear transform between the input layer and the first hidden layer of SymNetmk+1SymNet_{m}^{k+1} as

Then, the output of the first hidden layer is x1,x2,⋯ ,xm,x1x2x_{1},x_{2},\cdots,x_{m},x_{1}x_{2}. If we use it as the input of SymNetm+1kSymNet_{m+1}^{k}, we have

Let P∈Pk[x1,⋯ ,xm]P\in\mathcal{P}^{k}[x_{1},\cdots,x_{m}] and suppose PP have monomials of degree ≤l\leq l.

The memory load of SymNetmkSymNet_{m}^{k} that approximates PP is O(m+k)O(m+k). The number of flops for evaluating SymNetmkSymNet_{m}^{k} is O(k(m+k))O(k(m+k)).

Constructing a dictionary with all possible polynomials of degree ll requires a memory load of (m+ll)\binom{m+l}{l}, and evaluation of a linear combination of dictionary members requires O((m+ll))O(\binom{m+l}{l}) flops.

We use the following example to show the advantage of SymNetSymNet over SINDy. Consider two variables u,vu,v and all of their derivatives of order ≤2\leq 2:

Suppose the polynomial to be approximated is P=−uux−vuy+uxxP=-uu_{x}-vu_{y}+u_{xx}. For k=l=3k=l=3, the size of the dictionary of SINDy is (153)=455\binom{15}{3}=455 and the computation of linear combination of the elements requires 909909 flops. The memory load of SymNet123SymNet_{12}^{3}, however, is 15 and an evaluation of the network requires 180 flops. Therefore, SymNetSymNet can significantly reduce memory load and computation cost when input data is large. Note that when kk is large and ll small, SymNetmkSymNet_{m}^{k} is worse than SINDy. However, for system identification problems, we normally wish to obtain a compact representation (i.e. smaller kk). Thus, SymNetSymNet takes full advantage of this prior knowledge and can significantly save on memory and computation cost which is crucial in the training of the PDE-Net 2.0.

4 Loss Function and Regularization

We adopt the following loss function to train the proposed PDE-Net 2.0:

where the hyper-parameters λ1\lambda_{1} and λ2\lambda_{2} are chosen as λ1=0.001\lambda_{1}=0.001 and λ2=0.005\lambda_{2}=0.005. Now, we present details on each of the term of the loss function and introduce pseudo-upwind as an additional constraint on PDE-Net 2.0.

Consider the data set {Uj(ti,⋅):1≤i≤n,1≤j≤N}\{U_{j}(t_{i},\cdot):1\leq i\leq n,1\leq j\leq N\}, where nn is the number of δt\delta t-blocks and NN is the total number of samples. The index jj indicates the jj-th solution path with a certain initial condition of the unknown dynamics. Note that one can split a long solution path into multiple shorter ones. We would like to train the PDE-Net 2.0 with nn δt\delta t-blocks. For a given n≥1n\geq 1, every pair of the data {Uj(t0,⋅),Uj(ti,⋅)}\{U_{j}(t_{0},\cdot),U_{j}(t_{i},\cdot)\}, for each jj and i≤ni\leq n, is a training sample, where Uj(t0,⋅)U_{j}(t_{0},\cdot) is the input and Uj(ti,⋅)U_{j}(t_{i},\cdot) is the label that we need to match with the output from the network. For that, we define the data approximation term LdataL^{data} as:

where qij{q_{ij}} are the filters of PDE-Net 2.0 and M(q)M(q) is the moment matrix of qq. We use this loss function to regularize the learnable filters to reduce overfitting. In our numerical experiments, we will use s=0.01s=0.01.

We set s=0.001s=0.001 in our numerical experiments.

4.3 Pseudo-upwind

In numerical PDEs, to ensure stability of a numerical scheme, we need to design conservation schemes or use upwind schemes . This is also important for PDE-Net 2.0 during inferencing. However, the challenge we face is that we do not know apriori the form or the type of the PDE. Therefore, we introduce a method called pseudo-upwind to help with maintaining stability of the PDE-Net 2.0.

Given a 2D filter qq, define the flipping operators flipx(q)\text{flip}_{x}(q) and flipy(q)\text{flip}_{y}(q) as

In each δt\delta t-block of the PDE-Net 2.0, before we apply convolution with a filter, we first use SymNetSymNet to determine whether we should use the filter or flip it first. We use the following univariate PDE as an example to demonstrate our idea. Given PDE ut=F(u,⋯ )u_{t}=F(u,\cdots), suppose the input of a δt\delta t-block is uu. The algorithm of pseudo-upwind is described by the following Algorithm 2.

Note that the algorithm does not exactly enforce upwind in general. This is why we call it pseudo-upwind. We further note that:

Given a PDE of the form ut=G(u)ux+H(u)uy+λ(uxx+uyy)u_{t}=G(u)u_{x}+H(u)u_{y}+\lambda(u_{xx}+u_{yy}), we can use G(u)G(u) and H(u)H(u) to determine whether we should flip a filter or not.

5 Initialization and training

We use layer-wise training to train the PDE-Net 2.0. We start with training the PDE-Net 2.0 on the first δ\deltat-block with a batch of data, and then use the results of the first δ\deltat-block as the initialization and restart training on the first two δ\deltat-blocks with another batch. Repeat this procedure until we complete all nn δ\deltat-blocks. Note that all the parameters in each of the δ\deltat-block are shared across layers. In addition, we add a warm-up step before the training of the first δ\deltat-block by fixing filters and setting regularization term to be 0 (i.e. λ1=λ2=0\lambda_{1}=\lambda_{2}=0). The warm-up step is to obtain a good initial guess of the parameters of SymNetSymNet.

To demonstrate the necessity of having learnable filters, we will compare the PDE-Net 2.0 containing learnable filters with the PDE-Net 2.0 having fixed filters. To differentiate the two cases, we shall call the PDE-Net 2.0 with fixed filters the “Frozen-PDE-Net 2.0”. Note that for Frozen-PDE-Net 2.0, the filters are fixed to be the initial values we choose to train the regular PDE-Net 2.0. This is a natural choice since when we know apriori that the PDE is Burgers’ equation, it would be a stable finite difference scheme. However, intuitively speaking, freezing any finite difference approximations of the differential operators during training of PDE-Net 2.0 is not ideal, because you cannot possibly know which numerical scheme to use without knowing the form of the PDE. Therefore, for inverse problem, it is better to learn both the PDE model and the discretization of the PDE model simultaneously. This assertion is supported by our emperical comparisons between frozen and regular PDE-Net 2.0 in Table 1 and 3.

Numerical Studies: Burgers’ Equation

Burgers’ equation is a fundamental partial differential equation in many areas such as fluid mechanics and traffic flow modeling. It has a lot in common with the Navier-Stokes equation, e.g. the same type of advective nonlinearity and the presence of viscosity.

In this section we consider a 2-dimensional Burger’s equation with periodic boundary condition on Ω=[0,2π]×[0,2π]\Omega=[0,2\pi]\times[0,2\pi],

with (t,x,y)∈×Ω,(t,x,y)\in\times\Omega, where ν=0.05\nu=0.05

The training data is generated by a finite difference scheme on a 128×128128\times 128 mesh and then restricted to a 32×3232\times 32 mesh. The temporal discretization is 2nd order Runge-Kutta with time step δt=11600\delta t=\frac{1}{1600}, the spatial discretization uses a 2nd order upwind scheme for ∇\nabla and the central difference scheme for Δ\Delta. The initial value u0(x,y),v0(x,y)u_{0}(x,y),v_{0}(x,y) takes the following form,

where w0(x,y)=∑∣k∣,∣l∣≤4λk,lcos⁡(kx+ly)+γk,lsin⁡(kx+ly)w_{0}(x,y)=\sum_{|k|,|l|\leq 4}\lambda_{k,l}\cos(kx+ly)+\gamma_{k,l}\sin(kx+ly), λk,l,γk,l∼N(0,1),c∼U(−2,2)\lambda_{k,l},\gamma_{k,l}\sim\mathcal{N}(0,1),c\sim\mathcal{U}(-2,2). Here, N(0,1),U(−2,2)\mathcal{N}(0,1),\mathcal{U}(-2,2) represents the standard normal distribution and uniform distribution on $$ respectively. We also add noise to the generated data:

where M=max⁡x,y,t{U(t,x,y)}M=\max_{x,y,t}\{U(t,x,y)\}, W∼N(0,1)W\sim\mathcal{N}(0,1).

Suppose we know a priori that the order of the underlying PDE is no more than 2, we can use two SymNet125SymNet_{12}^{5} to approximate the right-hand-side nonlinear response function of (10) component-wise. Let U=(u,v)⊤U=(u,v)^{\top}. We denote the two SymNet125SymNet_{12}^{5} as NetuNet_{u} and NetvNet_{v} respectively. Then, each δt\delta t-block of the PDE-Net 2.0 can be written as

where {Dij:0≤i+j≤2}\{D_{ij}:0\leq i+j\leq 2\} are convolution operators.

During training and testing, the data is generated on-the-fly. The size of the filters that will be used is 5×55\times 5. The total number of parameters in NetuNet_{u} and NetvNet_{v} is 336, and the number of trainable parameters in moment matrices is 105 for 5×55\times 5 filters (6×5×5−456\times 5\times 5-45 (constraint on moment)). During training, we use BFGS, instead of SGD, to optimize the parameters. We use 28 data samples per batch to train each δt\delta t-block and we only construct the PDE-Net 2.0 up to 9 layers, which requires totally 420 data samples during the whole training procedure.

2 Results and discussions

We first demonstrate the ability of the trained PDE-Net 2.0 to recover the analytic form of the unknown PDE model. We use the symbolic math tool in python to obtain the analytic form of SymNetSymNet. Results are summarized in Table 1. As one can see from Table 1 that we can recover the terms of the Burgers’ equation with good accuracy, and using learnable filters helps with the identification of the PDE model. Furthermore, the terms that are not included in the Burgers’ equation all have relatively small weights in the SymNetSymNet (see Figure 7).

This subsection demonstrates the importance of enforcing sparsity on the SymNetSymNet and using pseudo-upwind. As we can see from Figure 7 that having sparsity constraints on the SymNetSymNet helps with suppressing the weights on the terms that do not exist in the Burgers’ equation. Furthermore, Figure 8 and Figure 9 show that having sparsity constraint on the SymNetSymNet or using pseudo-upwind can significantly reduce prediction errors.

Numerical Studies: Diffusion Equation

Diffusion phenomenon has been studied in many applications in physics e.g. the collective motion of micro-particles in materials due to random movement of each particle, or modeling the distribution of temperature in a given region over time.

Consider the 2-dimensional heat equation with periodic boundary condition on Ω=[0,2π]×[0,2π]\Omega=[0,2\pi]\times[0,2\pi]

where c=0.1c=0.1. The training data of the heat equation is generated by 2nd order Runge-Kutta in time with δt=11600\delta t=\frac{1}{1600}, and central difference scheme in space on a 128×128128\times 128 mesh. We then restrict the data to a 32×3232\times 32 mesh. The initial value u0(x,y)u_{0}(x,y) is also generated from (11).

2 Results and discussions

The demonstration on the ability of the trained PDE-Net 2.0 to identify the PDE model is given in Table 2. As one can see from Table 2 that we can recover the terms of the heat equation with good accuracy. Furthermore, all the terms that are not included in the heat equation have much smaller weights in the SymNetSymNet.

We also demonstrate the ability of the trained PDE-Net 2.0 in prediction. The testing method is exactly the same as the method described in Section 3. Comparisons between PDE-Net 2.0 and Frozen-PDE-Net 2.0 are shown in Figure 10, where we can clearly see the advantage of learning the filters. Visualization of the predicted dynamics is given in Figure 11. All these results show that the learned PDE-Net 2.0 performs well in prediction.

Numerical Studies: Convection Diffusion Equation with A Reactive Source

Convection diffusion systems are mathematical models which correspond to the transferring of some physical quantities such as energy or materials due to diffusion and convection. Specifically, a convection diffusion system with a reactive source can be used to model a large range of chemical systems in which the transferring of materials competes with productions of materials induced by several chemical reactions.

Consider a 2-dimensional convection diffusion equation with a reactive source and the periodic boundary condition on Ω=[0,2π]×[0,2π]\Omega=[0,2\pi]\times[0,2\pi]:

where (t,x,y)∈[0,1.5]×Ω,(t,x,y)\in[0,1.5]\times\Omega, and ν=0.1,β=1\nu=0.1,\beta=1. Training data is generated the same way as what we did for Burgers’ equation in Section 3. A 2nd order Runge-Kutta with time step δt=1/10000\delta t=1/10000 is adopted for temporal discretization. We choose a 2nd order upwind scheme for the convection terms and the central difference scheme for Δ\Delta on a 128×128128\times 128 mesh. We then restrict the data to a 32×3232\times 32 mesh. Noise is added the same way as the Burgers’ equation. The initial values u0(x,y),v0(x,y)u_{0}(x,y),v_{0}(x,y) are also generated from (11).

2 Results and discussions

The capability of the trained PDE-Net 2.0 to identify the underlying PDE model is demonstrated in Table 3. As one can see that we can recover the terms of the reaction convection diffusion equation with good accuracy. Furthermore, all the terms that are not included in this equation have relatively small weights in the SymNetSymNet.

We also demonstrate the ability of the trained PDE-Net 2.0 in prediction. The testing method is exactly the same as the method described in Section 3. Comparisons between PDE-Net 2.0 and Frozen-PDE-Net 2.0 are shown in Figure 12. Visualization of the predicted dynamics and errors maps are given in Figure 13 and Figure 14. Similar to what we observed in Section 3, we can clearly see the benefit from learning discretizations. PDE-Net 2.0 obtains more accurate estimations of the coefficients for the nonlinear convection terms (i.e. the term −uux−vuy-uu_{x}-vu_{y} in Table 3) and makes more accurate predictions (Figure 12) than Frozen-PDE-Net 2.0.

Conclusions and Future Work

In this paper, we proposed a numeric-symbolic hybrid deep network, called PDE-Net 2.0, for PDE model recovery from observed dynamic data. PDE-Net 2.0 is able to recover the analytic form of the PDE model with minor assumptions on the mechanisms of the observed dynamics. For example, it is able to recover the analytic form of Burgers’ equation with good confidence without any prior knowledge on the type of the equation. Therefore, PDE-Net 2.0 has the potential to uncover potentially new PDEs from observed data. Furthermore, after training, the network can perform accurate long-term prediction without re-training for new initial conditions. The limitations and possible future extensions of the current version of PDE-Net 2.0 is twofold: 1) having only addition and multiplication in the SymNetSymNet may still be too restrictive, and one may include division as an additional operation in SymNetSymNet to further improve its expressive power; 2) using forward Euler as temporal discretization is the most straightforward treatment, while a more sophisticated temporal scheme may help with the model recovery and prediction. Both of these worth further exploration. Furthermore, we would like to apply the network to real biological dynamic data. We will further investigate the reliability of the network and explore the possibility to uncover new dynamical principles that have meaningful scientific explanations.

Acknowledgments

Zichao Long is supported by The Elite Program of Computational and Applied Mathematics for PhD Candidates of Peking University. Yiping Lu is supported by the Elite Undergraduate Training Program of the School of Mathematical Sciences at Peking University. Bin Dong is supported in part by Beijing Natural Science Foundation (No. 180001) and Beijing Academy of Artificial Intelligence (BAAI).

References