Recent advances in machine learning in addition to new data recordings and sensor technologies have the potential to revolutionize our understanding of the physical world in modern application areas such as neuroscience, epidemiology, finance, and dynamic network analysis where first-principles derivations may be intractable [1]. In particular, many concepts from statistical learning can be integrated with classical methods in applied mathematics to help us discover sufficiently sophisticated and accurate mathematical models of complex dynamical systems directly from data. This integration of nonlinear dynamics and machine learning opens the door for principled methods for model construction, predictive modeling, nonlinear control, and reinforcement learning strategies. The literature on data-driven discovery of dynamical systems [2] is vast and encompasses equation-free modeling [3], ar-tificial neural networks [4, 5, 6, 7], nonlinear regression [8], empirical dynamic modeling [9, 10], modeling emergent behavior [11], automated inference of dynamics [12, 13, 14], normal form identification in climate [15], nonlinear Laplacian spectral analysis [16], modeling emergent behavior [11], Koopman analysis [17, 18, 19, 20], automated inference of dynamics [12, 13, 14], and symbolic regression [21, 22]. More recently, sparsity [23] has been used to determine the governing dynamical system [24, 25, 26, 27, 28, 29, 30, 31, 32, 33].
Less well studied is how to discover the underlying physical laws expressed by partial differential equations from scattered data collected in space and time. Inspired by recent developments in physics-informed deep learning [34, 35], we construct structured nonlinear regression models that can uncover the dynamic dependencies in a given set of spatio-temporal dataset, and return a closed form model that can be subsequently used to forecast future states. In contrast to recent approaches to systems identification [24, 1], here we do not need to have direct access or approximations to temporal or spatial derivatives. Moreover, we are using a richer class of function approximators to represent the nonlinear dynamics and consequently we do not have to commit to a particular family of basis functions. Specifically, we consider nonlinear partial differential equations of the general form
where N is a nonlinear function of time t, space x, solution u and its derivatives.Here, the subscripts denote partial differentiation in either time t or space
Given a set of scattered and potentially noisy observations of the solution u, we are interested in learning the nonlinear function N and consequently discovering the hidden laws of physics that govern the evolution of the observed data.
For instance, let us assume that we would like to discover the Burger’s equation [36] in one space dimension . Although not pursued in the current work, a viable approach [1] to tackle this problem is to create a dictionary of possible terms and write the following expansion
Given the aforementioned large collection of candidate terms for constructing the partial differential equation, one could then use sparse regression techniques [1] to determine the coefficients and consequently the right-hand-side terms that are contributing to the dynamics. A huge advantage of this approach is the interpretability of the learned equations. However, there are two major drawbacks with this method.
First, it relies on numerical differentiation to compute the derivatives involved in equation (1). Derivatives are taken either using finite differences for clean data or with polynomial interpolation in the presence of noise. Numerical approximations of derivatives are inherently ill-conditioned and unstable [37] even in the absence of noise in the data. This is due to the introduction of truncation and round-off errors inflicted by the limited precision of computations and the chosen value of the step size for finite differencing. Thus, this approach requires far more data points than library functions. This need for using a large number of points lies more in the numerical evaluation of derivatives than in supplying sufficient data for the regression.
Second, in applying the algorithm [1] outlined above we assume that the chosen library is sufficiently rich to have a sparse representation of the time dynamics of the dataset. However, when applying this approach to a dataset where the dynamics are in fact unknown it is not unlikely that the basis chosen above is insufficient. Specially, in higher dimensions (i.e., for input x or output u) the required number of terms to include in the library increases exponentially. Moreover, an additional issue with this approach is that it can only estimate parameters appearing as coefficients. For example, this method cannot estimate parameters of a partial differential equation (e.g., the sine-Gordon equation) involving a term like sin(unknown parameter, even if we include sines and cosines in the dictionary of possible terms.
One could avoid the first drawback concerning numerical differentiation by assigning prior distributions in the forms of Gaussian processes [38, 39, 40, 41] or neural networks [35, 34] to the unknown solution u. Derivatives of the prior on u can now be evaluated at machine precision using symbolic or automatic differentiation [37]. This removes the requirement for having or generating data on derivatives of the solution u. This is enabling as it allows us to work with noisy observations of the solution u, scattered in space and time. Moreover, this approach requires far fewer data points than the method proposed in [1] simply because, as explained above, the need for using a large number of data points was due to the numerical evaluation of derivatives.
The second drawback can be addressed in a similar fashion by approximating the nonlinear function N (see equation (1)) with a neural network. Representing the nonlinear function N by a deep neural network is the novelty of the current work. Deep neural networks are a richer family of function approximators and consequently we do not have to commit to a particular class of basis functions such as polynomials or sines and cosines. This expressiveness comes at the cost of losing interpretability of the learned dynamics. However, there is nothing hindering the use of a particular class of basis functions in order obtain more interpretable equations [35].
We proceed by approximating both the solution u and the nonlinear function N with two deep neural networksand define a deep hidden physics model f to be given by
We obtain the derivatives of the neural network u with respect to time t and space x by applying the chain rule for differentiating compositions of functions using automatic differentiation [37]. It is worth emphasizing that automatic differentiation is different from, and in several aspects superior to, numerical or symbolic differentiation; two commonly encountered techniques of computing derivatives. In its most basic description [37], automatic differentiation relies on the fact that all numerical computations are ultimately compositions of a finite set of elementary operations for which derivatives are known. Combining the derivatives of the constituent operations through the chain rule gives the derivative of the overall composition. This allows accurate evaluation of derivatives at machine precision with ideal asymptotic efficiency and only a small constant factor of overhead. In particular, to compute the derivatives involved in definition (2) we rely on Tensorflow [42] which is a popular and relatively well documented open source software library for automatic differentiation and deep learning computations.
Parameters of the neural networks u and N can be learned by minimizing the sum of squared errors
where denote the training data on
tries to fit the data by adjusting the parameters of the neural network u while the term
learns the parameters of the network N by trying to satisfy the partial differential equation (1) at the collocation points (
Training the parameters of the neural networks u and N can be performed simultaneously by minimizing the sum of squared error (3) or in a sequential fashion by training u first and N second.
How can we make sure that the algorithm presented above results in an acceptable function N? One answer would be to solve the learned equations and compare the resulting solution to the solution of the exact partial differential equation. However, it should be pointed out that the learned function N is a black-box function; i.e., we do not know its functional form. Consequently, none of the classical partial differential equation solvers such as finite differences, finite elements or spectral methods are applicable here. Therefore, to solve the learned equations we have no other choice than to resort to modern black-box solvers such as physic informed neural networks (PINNs) introduced in [34]. The steps involved in PINNs as solverssimilar to equations (1), (2), and (3) with the nonlinear function N being known and the data residing on the boundary of the domain.
The proposed framework provides a universal treatment of nonlinear partial differential equations of fundamentally different nature. This generality will be demonstrated by applying the algorithm to a wide range of canonical problems spanning a number of scientific domains including the Burgers’, Korteweg-de Vries (KdV), Kuramoto-Sivashinsky, nonlinear Schr¨odinger, and Navier-Stokes equations. These examples are motivated by the pioneering work of [1]. All data and codes used in this manuscript are publicly available on GitHub at https://github.com/maziarraissi/ DeepHPMs.
3.1. Burgers’ equation
Let us start with the Burgers’ equation arising in various areas of engineering and applied mathematics, including fluid mechanics, nonlinear acoustics, gas dynamics, and traffic flow [36]. In one space dimension the Burgers’
equation reads as
To obtain a set of training and test data we simulate the Burger’s equation (4) using conventional spectral methods. Specifically, starting from an initial condition 8] and assuming periodic boundary conditions, we integrate equation (4) up to the final time t = 10. We use the Chebfun package [43] with a spectral Fourier discretization with 256 modes and a fourth-order explicit Runge-Kutta temporal integrator with time-step size 10
The solution is saved every ∆t = 0.05 to give us a total of 201 snapshots. Out of this data-set, we generate a smaller training subset, scattered in space and time, by randomly sub-sampling 10000 data points from time t = 0 to t = 6.7. We call the portion of the domain from time t = 0 to t = 6.7 the training portion. The rest of the domain from time t = 6.7 to the final time t = 10 will be referred to as the test portion. Using this terminology, we are in fact sub-sampling from the original dataset only in the training portion of the domain. Given the training data, we are interested in learning N as a function of the solution u and its derivatives up to the 2nd order
We represent the solution u by a 5-layer deep neural network with 50 neurons per hidden layer. Furthermore, we let N to be a neural network with 2 hidden layers and 100 neurons per hidden layer. These two networks are trained by minimizing the sum of squared errors loss of equation (3). To illustrate the effectiveness of our approach, we solve the learned partial differential equation (5), along with periodic boundary conditions and the same initial condition as the one used to generate the original dataset, using the PINNs algorithm [34]. The original dataset alongside the resulting solution of the learned partial differential equation are depicted in figure 1. This fig-ure indicates that our algorithm is able to accurately identify the underlying partial differential equation with a relative -error of 4.78e-03. It should be highlighted that the training data are collected in roughly two-thirds of the domain between times t = 0 and t = 6.7. The algorithm is thus extrapolating from time t = 6.7 onwards. The relative
-error on the training portion of
Figure 1: Burgers equation: A solution to the Burger’s equation (left panel) is compared to the corresponding solution of the learned partial differential equation (right panel). The identified system correctly captures the form of the dynamics and accurately reproduces the solution with a relative -error of 4.78e-03. It should be emphasized that the training data are collected in roughly two-thirds of the domain between times t = 0 and t = 6.7 represented by the white vertical lines. The algorithm is thus extrapolating from time t = 6.7 onwards. The relative
-error on the training portion of the domain is 3.89e-03.
the domain is 3.89e-03.
Furthermore, we performed a systematic study of the reported results in figure 1 with respect to noise levels in the data by keeping the total number of training observations as well as the neural network architectures fixed to the settings described above. In particular, we added white noise with magnitude equal to one, two, and five percent of the standard deviation of the solution function. The results of this study are summarized in table 1. The key observation here is that less noise in the data enhances the performance of the algorithm. Our experience so far indicates that the negative consequences of more noise in the data can be remedied to some extent by obtaining more data. Another fundamental point to make is that in general the choice of the neural network architectures, i.e., activation functions and number of layers/neurons, is crucial. However, in many cases of practical interest this choice still remains an art and systematic studies with respect to the neural network architectures often fail to reveal consistent patterns [34, 35, 4]. We usually observe some variability and non monotonic trends in systematic studies pertaining to the network architectures. In this regard, there exist a series of open questions mandating further investigations. For instance, how common techniques such as batch normalization, drop out, and regularization [44] could enhance the robustness of the proposed algorithm with respect to the neural network architectures as well as noise in the data.
To further scrutinize the performance of the algorithm, let us change the initial condition to ) and solve the Burgers’ equation (4) using a classical partial differential equation solver. In particular, the new data-set [1] contains 101 time snapshots of a solution to the Burgers’ equation (4) with a Gaussian initial condition, propagating into a traveling wave. The snapshots are ∆t = 0.1 apart and stretch from time t = 0 to t = 10. The spatial discretization of each snapshot involves a uniform grid with 256 cells. We compare the resulting solution to the one obtained by solving the learned partial differential equation (5) using the PINNs algorithm [34]. It is worth emphasizing that the algorithm is trained on the dataset depicted in figure 1 and is being tested on a totally different dataset as shown in figure 2. The surprising result reported in figure 2 is a strong indication that the algorithm is capable of accurately identifying the underlying partial differential equation. The algorithm hasn’t seen even a single observation of the dataset shown in figure 2 during model training and is yet capable of achieving a relatively accurate approximation of the true solution. The identified system reproduces the solution to the Burgers’ equation with a relative
7.33e-02.
However, the aforementioned result seems to be sensitive to making the nonlinear function N of equation (5) depend on either time t or space x, or both. For instance, if we look for equations of the form the relative
-error between the exact and the learned solutions corresponding to figure 2 increases to 4.25e-01. Similarly, if we look for equations of the form
), the relative
-error increases to 2.58e-01. Moreover, if we make the nonlinear function N both time and space dependent
Table 1: -error between solutions of the Burgers’ equation and the learned partial differential equation as a function of noise corruption level in the data. Here, the total number of training data as well as the neural network architectures are kept fixed.
Figure 2: Burgers equation: A solution to the Burger’s equation (left panel) is compared to the corresponding solution of the learned partial differential equation (right panel). The identified system correctly captures the form of the dynamics and accurately reproduces the solution with a relative -error of 7.33e-02. It should be highlighted that the algorithm is trained on a dataset totally different from the one used in this figure.
the relative -error increases even further to 1.46e+00. A possible explanation for this behavior could be that some of the dynamics in the training data is now being explained by time t or space x. This makes the algorithm over-fit to the training data and consequently harder for it to generalize to datasets it has never seen before. Based on our experience more data seems to help resolve this issue. Another interesting observation is that if we train on the dataset presented in figure 2 and test on the one depicted in figure 1, i.e., swap the roles of these two datasets, the algorithm fails to generalize to the new test data. This observation suggests that the dynamics (see figure 1) generated by
8) as the initial condition is a reacher one compared to the dynamics (see figure 2) generated by exp(
also evident from a visual comparison of the datasets given in figures 1 and 2.
Let us now take a closer look at equation (5) and ask: what would happen if we included derivatives of higher order than two in our formulation? As will be demonstrated in the following, the algorithm proposed in the current work is capable of handling such cases, however, this is a fundamental question worthy of a moment of reflection. The choice of the order of the partial differential equation (5) determines the form and the number of boundary conditions needed to end up with a well-posed problem. For instance, in the one dimensional setting of equation (5) including a third order derivative requires three boundary conditions, namely 8), assuming period boundary conditions. Consequently, in cases of practical interest, the available information on the boundary of the domain could help us determine the order of the partial differential equation we are try to identify. With this in mind, let us study the robustness of our framework with respect to the highest order of the derivatives included in equation (5). As for the boundary conditions, we use
8) when solving the identified partial differential equation regardless of the initial choice of its order. The results are summarized in table 2. The first column of table 2 demonstrates that a single first order derivative is clearly not enough to capture the second order dynamics of the Burgers’ equation. Moreover, the method seems to be generally robust with respect to the number and order of derivatives included in the nonlinear function N. Therefore, in addition to any information residing on the domain boundary, studies such as table 2, albeit for training or validation datasets, could help us choose the best order for the underlying partial differential equation. In this case, table 2 suggests the order of the equation to be two.
In addition, it must be stated that including higher order derivatives comes at the cost of reducing the speed of the algorithm due to the growth in the complexity of the resulting computational graph for the corresponding deep hidden physics model (see equation (2)). Also, another drawback is that higher order derivatives are usually less accurate specially if one uses single precision floating-point system (float32) rather than double precision (float64). It is true that automatic differentiation enables us to evaluate derivatives at machine precision, however, for float32 the machine epsilon is approximately 1.19e-07. For improved performance in terms of speed of the algorithm and constrained by usual GPU (graphics processing unit) platforms we often end up using float32 which boils down to committing an error of approximately 1.19e-07 while computing the required derivatives. The aforementioned issues do not seem to be too serious since computer infrastructure (both hardware and software) for deep learning is constantly improving.
3.2. The KdV equation
As a mathematical model of waves on shallow water surfaces one could consider the Korteweg-de Vries (KdV) equation. The KdV equation reads as
To obtain a set of training data we simulate the KdV equation (6) using conventional spectral methods. In particular, we start from an initial condition 20] and assume periodic boundary conditions. We integrate equation (6) up to the final time t = 40. We use the Chebfun package [43] with a spectral Fourier discretization with 512 modes and a fourth-order explicit Runge-Kutta temporal integrator with time-step size 10
. The solution is saved every ∆t = 0.2 to give us a total of 201 snapshots. Out of this data-set, we generate a smaller training subset, scattered in space and time, by randomly sub-sampling 10000 data points from time t = 0 to t = 26.8. In other words, we are sub-sampling from the original dataset only in the training portion of the domain from time t = 0 to t = 26.8. Given the training data, we are interested in learning N as a function of the solution u and its derivatives up to the 3rd order
We represent the solution u by a 5-layer deep neural network with 50 neurons per hidden layer. Furthermore, we let N to be a neural network with 2 hidden layers and 100 neurons per hidden layer. These two networks are trained by minimizing the sum of squared errors loss of equation (3). To illustrate the ef-
A detailed study of the choice of the order is provided in section 3.1 for the Burgers’ equation.
Table 2: -error between solutions of the Burgers’ equation and the learned partial differential equation as a function of the highest order of spatial derivatives included in our formulation. For instance, the case corresponding to the 3rd order means that we are looking for a nonlinear function N such that
). Here, the total number of training data as well as the neural network architectures are kept fixed and the data are assumed to be noiseless.
Figure 3: The KdV equation: A solution to the KdV equation (left panel) is compared to the corresponding solution of the learned partial differential equation (right panel). The identified system correctly captures the form of the dynamics and accurately reproduces the solution with a relative -error of 6.28e-02. It should be emphasized that the training data are collected in roughly two-thirds of the domain between times t = 0 and t = 26.8 represented by the white vertical lines. The algorithm is thus extrapolating from time t = 26.8 onwards. The relative
-error on the training portion of the domain is 3.78e-02.
fectiveness of our approach, we solve the learned partial differential equation (7) using the PINNs algorithm [34]. We assume periodic boundary conditions and the same initial condition as the one used to generate the original dataset. The resulting solution of the learned partial differential equation as well as the exact solution of the KdV equation are depicted in figure 3. This figure indicates that our algorithm is capable of accurately identifying the underlying partial differential equation with a relative -error of 6.28e-02. It should be highlighted that the training data are collected in roughly two-thirds of the domain between times t = 0 and t = 26.8. The algorithm is thus extrapolating from time t = 26.8 onwards. The corresponding relative
-error on the training portion of the domain is 3.78e-02.
To test the algorithm even further, let us change the initial condition to cos(20) and solve the KdV (6) using the conventional spectral method outlined above. We compare the resulting solution to the one obtained by solving the learned partial differential equation (5) using the PINNs algorithm [34]. It is worth emphasizing that the algorithm is trained on the dataset depicted in figure 3 and is being tested on a different dataset as shown in figure 4. The surprising result reported in figure 4 strongly indicates that the algorithm is accurately learning the underlying partial differ-
Figure 4: The KdV equation: A solution to the KdV equation (left panel) is compared to the corresponding solution of the learned partial differential equation (right panel). The identified system correctly captures the form of the dynamics and accurately reproduces the solution with a relative -error of 3.44e-02. It should be highlighted that the algorithm is trained on a dataset totally different from the one shown in this figure.
ential equation; i.e., the nonlinear function N. The algorithm hasn’t seen the dataset shown in figure 4 during model training and is yet capable of achieving a relatively accurate approximation of the true solution. To be precise, the identified system reproduces the solution to the KdV equation with a relative -error of 3.44e-02.
3.3. Kuramoto-Sivashinsky equation
As a canonical model of a pattern forming system with spatio-temporal chaotic behavior we consider the Kuramoto-Sivashinsky equation. In one space dimension the Kuramoto-Sivashinsky equation reads as
We generate a dataset containing a direct numerical solution of the KuramotoSivashinsky (8) equation with 512 spatial points and 251 snapshots. To be precise, assuming periodic boundary conditions, we start from the initial condition 10] and integrate equation (8) up to the final time t = 50. We employ the Chebfun package [43] with a spectral Fourier discretization with 512 modes and a fourth-order explicit Runge-Kutta temporal integrator with time-step size 10
. The snapshots are saved every ∆t = 0.2. From this dataset, we create a smaller training subset, scattered in space and time, by randomly sub-sampling 10000 data points from time t = 0 to the final time t = 50.0. Given the resulting training
Figure 5: Kuramoto-Sivashinsky equation: A solution to the Kuramoto-Sivashinsky equation (left panel) is compared to the corresponding solution of the learned partial differential equation (right panel). The identified system correctly captures the form of the dynamics and reproduces the solution with a relative -error of 7.63e-02.
data, we are interested in learning N as a function of the solution u and its derivatives up to the 4rd order
We let the solution u to be represented by a 5-layer deep neural network with 50 neurons per hidden layer. Furthermore, we approximate the nonlinear function N by a neural network with 2 hidden layers and 100 neurons per hidden layer. These two networks are trained by minimizing the sum of squared errors loss of equation (3). To demonstrate the effectiveness of our approach, we solve the learned partial differential equation (9) using the PINNs algorithm [34]. We assume the same initial and boundary conditions as the ones used to generate the original dataset. The resulting solution of the learned partial differential equation alongside the exact solution of the Kuramoto-Sivashinsky equation are depicted in figure 5. This figure indicates that our algorithm is capable of identifying the underlying partial differential equation with a relative -error of 7.63e-02.
However, it must be mentioned that we are avoiding the regimes where the Kuramoto-Sivashinsky equation becomes chaotic. For instance, by changing
Figure 6: Kuramoto-Sivashinsky equation: A solution to the Kuramoto-Sivashinsky equation. This problem remains stubbornly unsolved in the face of the algorithm proposed in the current work.
the domain to [0] and the initial condition to cos(x/16)(1 + sin(x/16)) and by integrating the Kuramoto-Sivashinsky equation (8) up to the final time t = 100, one could end up with a relatively complicated solution as depicted in figure 6. This problem remains stubbornly unsolved in the face of the algorithm proposed in the current work as well as the PINNs framework introduced in [34, 35]. It is not difficult for a plain vanilla neural network to approximate the function depicted in figure 6. However, as we compute the derivatives required in equation (2), minimizing the loss function (3) becomes a challenge. Understanding what makes this problem hard to solve should be at the core of future extensions of this line of research.
3.4. Nonlinear Schr¨odinger equation
The one-dimensional nonlinear Schr¨odinger equation is a classical field equation that is used to study nonlinear wave propagation in optical fibers and/or waveguides, Bose-Einstein condensates, and plasma waves. This example aims to highlight the ability of our framework to handle complexvalued solutions as well as different types of nonlinearities in the governing partial differential equations. The nonlinear Schr¨odinger equation is given by
Let u denote the real part of the imaginary part. Then, the nonlinear Schr¨odinger equation can be equivalently written as a system of partial
differential equations
In order to assess the performance of our method, we simulate equation (10) using conventional spectral methods to create a high-resolution data set. Specifically, starting from an initial state ) and assuming periodic boundary conditions
integrate equation (10) up to a final time
2 using the Chebfun package [43]. We are in fact using a spectral Fourier discretization with 256 modes and a fourth-order explicit Runge-Kutta temporal integrator with time-step ∆
. The solution is saved approximately every ∆t = 0.0031 to give us a total of 501 snapshots. Out of this data-set, we generate a smaller training subset, scattered in space and time, by randomly sub-sampling 10000 data points from time
2. Given the resulting training data, we are interested in learning two nonlinear functions
as functions of the solutions u, v and their derivatives up to the 2nd
order
We represent the solutions u and v by two 5-layer deep neural networks with 50 neurons per hidden layer. Furthermore, we let to be two neural networks with 2 hidden layers and 100 neurons per hidden layer. These four networks are trained by minimizing the sum of squared errors loss of equation (3). To illustrate the effectiveness of our approach, we solve the learned partial differential equation (12), along with periodic boundary conditions and the same initial condition as the one used to generate the original dataset, using the PINNs algorithm [34]. The original dataset (in absolute values, i.e.,
) alongside the resulting solution (also in absolute values) of the learned partial differential equation are depicted in figure 7. This figure indicates that our algorithm is able to accurately identify the underlying partial differential equation with a relative
-error of 6.28e-03.
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
Figure 7: Absolute value of a solution to the nonlinear Schr¨odinger equation (left panel) is compared to the absolute value of the corresponding solution of the learned partial differential equation (right panel). The identified system correctly captures the form of the dynamics and accurately reproduces the absolute value of the solution with a relative
-error of 6.28e-03.
3.5. Navier-Stokes equation
Let us consider the Navier-Stokes equation in two dimensionsexplicitly by
where w denotes the vorticity, u the x-component of the velocity field, and v the y-component. To generate a training dataset for this problem we follow the exact same instructions as the ones provided in [45, 1]. Specifically, we simulate the Navier-Stokes equations describing the two-dimensional fluid flow past a circular cylinder at Reynolds number 100 using the Immersed Boundary Projection Method [46, 47]. This approach utilizes a multi-domain scheme with four nested domains, each successive grid being twice as large as the previous one. Length and time are nondimensionalized so that the cylinder has unit diameter and the flow has unit velocity. Data is collected on the finest domain with dimensions 9 4 at a grid resolution of 449
The flow solver uses a 3rd-order Runge Kutta integration scheme with a time step of t = 0.02, which has been verified to yield well-resolved and converged flow fields. After simulation converges to steady periodic vortex shedding, 151 flow snapshots are saved every ∆t = 0.02. We use a small portion of
Figure 8: Navier-Stokes equation: A snapshot of the vorticity field of a solution to the Navier-Stokes equations for the fluid flow past a cylinder. The dashed red box in this panel specifies the sampling region.
the resulting data-set for model training. In particular, we subsample 50000 data points, scattered in space and time, in the rectangular region (dashed red box) downstream of the cylinder as shown in figure 8.
Given the training data, we are interested in learning N as a function of the stream-wise u and transverse v velocity components in addition to the vorticity w and its derivatives up to the 2nd order
We represent the solution w by a 5-layer deep neural network with 200 neurons per hidden layer. Furthermore, we let N to be a neural network with 2 hidden layers and 100 neurons per hidden layer. These two networks are trained by minimizing the sum of squared errors loss
Figure 9: Navier-Stokes equation: A randomly picked snapshot of a solution to the NavierStokes equation (left panel) is compared to the corresponding snapshot of the solution of the learned partial differential equation (right panel). The identified system correctly captures the form of the dynamics and accurately reproduces the solution with a relative -error of 5.79e-03 in space and time.
given by
To illustrate the effectiveness of our approach, we solve the learned partial differential equation (14), in the region specified in figure 8 by the dashed red box, using the PINNs algorithm [34]. We use the exact solution to provide us with the required Dirichlet boundary conditions as well as the initial condition needed to solve the leaned partial differential equation (14). A randomly picked snapshot of the vorticity field in the original dataset alongside the corresponding snapshot of the solution of the learned partial differential equation are depicted in figure 9. This figure indicates that our algorithm is able to accurately identify the underlying partial differential equation with a relative -error of 5.79e-03 in space and time.
We have presented a deep learning approach for extracting nonlinear partial differential equations from spatio-temporal datasets. The proposed algorithm leverages recent developments in automatic differentiation to construct efficient algorithms for learning infinite dimensional dynamical systems using deep neural networks. In order to validate the performance of our approach we had no other choice than to rely on black-box solvers (see e.g., [35]). This signifies the importance of developing general purpose partial differen-tial equation solvers. Developing these types of solvers is still in its infancy and more collaborative work is needed to bring them to the maturity level of conventional methods such as finite elements, finite differences, and spectral methods which have been around for more than half a century or so.
There exist a series of open questions mandating further investigations. For instance, many real-world partial differential equations depend on parameters and, when the parameters are varied, they may go through bifurcations (e.g., the Reynold number for the Navier-Stokes equations). Here, the goal would be to collect data from the underlying dynamics corresponding to various parameter values, and infer the parameterized partial differ-ential equation. Another exciting avenue of future research would be to apply convolutional architectures [44] for mitigating the complexity associated with partial differential equations with very high-dimensional inputs. These types of equations appear routinely in dynamic programming, optimal control, or reinforcement learning. Moreover, a quick look at the list of nonlinear partial differential equations on Wikipedia reveals that many of these equations take the form specified in equation (1). However, a handful of them do not take this form, including the Boussinesq type equation framework outlined in the current work to incorporate all such cases. In the end, it is not always clear what measurements of a dynamical system to take. Even if we did know, collecting these measurements might be prohibitively expensive. It is well-known that time-delay coordinates of a single variable can act as additional variables. It might be interesting to investigating this idea for the infinite dimensional setting of partial differential equations.
This work received support by the DARPA EQUiPS grant N66001-15-2-4055 and the AFOSR grant FA9550-17-1-0013. All data and codes used in this manuscript are publicly available on GitHub at https://github.com/ maziarraissi/DeepHPMs.
[1] S. H. Rudy, S. L. Brunton, J. L. Proctor, J. N. Kutz, Data-driven discovery of partial differential equations, Science Advances 3 (2017).
[2] J. P. Crutchfield, B. S. McNamara, Equations of motion from a data series, Complex systems 1 (1987) 121.
[3] I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidid, O. Run- borg, C. Theodoropoulos, et al., Equation-free, coarse-grained multi-scale computation: Enabling mocroscopic simulators to perform systemlevel analysis, Communications in Mathematical Sciences 1 (2003) 715– 762.
[4] M. Raissi, P. Perdikaris, G. E. Karniadakis, Multistep neural networks for data-driven discovery of nonlinear dynamical systems, arXiv preprint arXiv:1801.01236 (2018).
[5] R. Gonzalez-Garcia, R. Rico-Martinez, I. Kevrekidis, Identification of distributed parameter systems: A neural net based approach, Computers & chemical engineering 22 (1998) S965–S968.
[6] J. Anderson, I. Kevrekidis, R. Rico-Martinez, A comparison of recurrent training algorithms for time series analysis and system identification, Computers & chemical engineering 20 (1996) S751–S756.
[7] R. Rico-Martinez, K. Krischer, I. Kevrekidis, M. Kube, J. Hudson, Discrete-vs. continuous-time nonlinear signal processing of cu electrodissolution data, Chemical Engineering Communications 118 (1992) 25–48.
[8] H. U. Voss, P. Kolodner, M. Abel, J. Kurths, Amplitude equations from spatiotemporal binary-fluid convection data, Physical review letters 83 (1999) 3422.
[9] G. Sugihara, R. May, H. Ye, C.-h. Hsieh, E. Deyle, M. Fogarty, S. Munch, Detecting causality in complex ecosystems, science 338 (2012) 496–500.
[10] H. Ye, R. J. Beamish, S. M. Glaser, S. C. Grant, C.-h. Hsieh, L. J. Richards, J. T. Schnute, G. Sugihara, Equation-free mechanistic ecosystem forecasting using empirical dynamic modeling, Proceedings of the National Academy of Sciences 112 (2015) E1569–E1576.
[11] A. J. Roberts, Model emergent dynamics in complex systems, SIAM, 2014.
[12] M. D. Schmidt, R. R. Vallabhajosyula, J. W. Jenkins, J. E. Hood, A. S. Soni, J. P. Wikswo, H. Lipson, Automated refinement and inference of analytical models for metabolic networks, Physical biology 8 (2011) 055011.
[13] B. C. Daniels, I. Nemenman, Automated adaptive inference of phenomenological dynamical models, Nature communications 6 (2015).
[14] B. C. Daniels, I. Nemenman, Efficient inference of parsimonious phe- nomenological models of cellular dynamics using s-systems and alternating regression, PloS one 10 (2015) e0119821.
[15] A. J. Majda, C. Franzke, D. Crommelin, Normal forms for reduced stochastic climate models, Proceedings of the National Academy of Sciences 106 (2009) 3649–3653.
[16] D. Giannakis, A. J. Majda, Nonlinear laplacian spectral analysis for time series with intermittency and low-frequency variability, Proceedings of the National Academy of Sciences 109 (2012) 2222–2227.
[17] I. Mezi´c, Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynamics 41 (2005) 309–325.
[18] M. Budiˇsi´c, R. Mohr, I. Mezi´c, Applied koopmanism a, Chaos: An Interdisciplinary Journal of Nonlinear Science 22 (2012) 047510.
[19] I. Mezi´c, Analysis of fluid flows via spectral properties of the koopman operator, Annual Review of Fluid Mechanics 45 (2013) 357–378.
[20] S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, J. N. Kutz, Chaos as an intermittently forced linear system, Nature Communications 8 (2017).
[21] J. Bongard, H. Lipson, Automated reverse engineering of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 104 (2007) 9943–9948.
[22] M. Schmidt, H. Lipson, Distilling free-form natural laws from experi- mental data, science 324 (2009) 81–85.
[23] R. Tibshirani, Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society. Series B (Methodological) (1996) 267– 288.
[24] S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equa- tions from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113 (2016) 3932–3937.
[25] N. M. Mangan, S. L. Brunton, J. L. Proctor, J. N. Kutz, Inferring biological networks by sparse identification of nonlinear dynamics, IEEE Transactions on Molecular, Biological and Multi-Scale Communications 2 (2016) 52–63.
[26] W.-X. Wang, R. Yang, Y.-C. Lai, V. Kovanis, C. Grebogi, Predicting catastrophes in nonlinear dynamical systems by compressive sensing, Physical Review Letters 106 (2011) 154101.
[27] H. Schaeffer, R. Caflisch, C. D. Hauck, S. Osher, Sparse dynamics for partial differential equations, Proceedings of the National Academy of Sciences 110 (2013) 6634–6639.
[28] V. Ozoli¸nˇs, R. Lai, R. Caflisch, S. Osher, Compressed modes for varia- tional problems in mathematics and physics, Proceedings of the National Academy of Sciences 110 (2013) 18368–18373.
[29] A. Mackey, H. Schaeffer, S. Osher, On the compressive spectral method, Multiscale Modeling & Simulation 12 (2014) 1800–1827.
[30] S. L. Brunton, J. H. Tu, I. Bright, J. N. Kutz, Compressive sensing and low-rank libraries for classification of bifurcation regimes in nonlinear dynamical systems, SIAM Journal on Applied Dynamical Systems 13 (2014) 1716–1732.
[31] J. L. Proctor, S. L. Brunton, B. W. Brunton, J. Kutz, Exploiting spar- sity and equation-free architectures in complex systems, The European Physical Journal Special Topics 223 (2014) 2665–2684.
[32] Z. Bai, T. Wimalajeewa, Z. Berger, G. Wang, M. Glauser, P. K. Varsh- ney, Low-dimensional approach for reconstruction of airfoil data via compressive sensing, AIAA Journal (2014).
[33] G. Tran, R. Ward, Exact recovery of chaotic systems from highly cor- rupted data, arXiv preprint arXiv:1607.01067 (2016).
[34] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics informed deep learning (part ii): Data-driven discovery of nonlinear partial differential equations, arXiv preprint arXiv:1711.10566 (2017).
[35] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations, arXiv preprint arXiv:1711.10561 (2017).
[36] C. Basdevant, M. Deville, P. Haldenwang, J. Lacroix, J. Ouazzani, R. Peyret, P. Orlandi, A. Patera, Spectral and finite difference solutions of the Burgers equation, Computers & fluids 14 (1986) 23–41.
[37] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, J. M. Siskind, Automatic differentiation in machine learning: a survey, arXiv preprint arXiv:1502.05767 (2015).
[38] M. Raissi, P. Perdikaris, G. E. Karniadakis, Machine learning of linear differential equations using Gaussian processes, Journal of Computational Physics 348 (2017) 683 – 693.
[39] M. Raissi, G. E. Karniadakis, Hidden physics models: Machine learning of nonlinear partial differential equations, arXiv preprint arXiv:1708.00588 (2017).
[40] M. Raissi, P. Perdikaris, G. E. Karniadakis, Inferring solutions of dif- ferential equations using noisy multi-fidelity data, Journal of Computational Physics 335 (2017) 736–746.
[41] M. Raissi, P. Perdikaris, G. E. Karniadakis, Numerical Gaussian pro- cesses for time-dependent and non-linear partial differential equations, arXiv preprint arXiv:1703.10230 (2017).
[42] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, et al., Tensorflow: Large-scale machine learning on heterogeneous distributed systems, arXiv preprint arXiv:1603.04467 (2016).
[43] T. A. Driscoll, N. Hale, L. N. Trefethen, Chebfun guide, 2014.
[44] I. Goodfellow, Y. Bengio, A. Courville, Deep learning, MIT press, 2016.
[45] J. N. Kutz, S. L. Brunton, B. W. Brunton, J. L. Proctor, Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems, volume 149, SIAM, 2016.
[46] K. Taira, T. Colonius, The immersed boundary method: a projection approach, Journal of Computational Physics 225 (2007) 2118–2137.
[47] T. Colonius, K. Taira, A fast immersed boundary method using a nullspace approach and multi-domain far-field boundary conditions, Computer Methods in Applied Mechanics and Engineering 197 (2008) 2131–2146.