Numerical Aspects for Approximating Governing Equations Using Data

2018·Arxiv

Abstract

Abstract

We present effective numerical algorithms for locally recovering unknown governing differential equations from measurement data. We employ a set of standard basis functions, e.g., polynomials, to approximate the governing equation with high accuracy. Upon recasting the problem into a function approximation problem, we discuss several important aspects for accurate approximation. Most notably, we discuss the importance of using a large number of short bursts of trajectory data, rather than using data from a single long trajectory. Several options for the numerical algorithms to perform accurate approximation are then presented, along with an error estimate of the final equation approximation. We then present an extensive set of numerical examples of both linear and nonlinear systems to demonstrate the properties and effectiveness of our equation recovery algorithms.

Keywords: Ordinary differential equation, differential-algebraic equation, measurement data, data-driven discovery, regression, sequential approximation.

1. Introduction

Recently there has been a growing interest in discovering governing equations of certain physical problems using observational data. It is a result of the wide recognition that all physical laws, e.g., Newton’s law, Kepler’s law, etc., were established based on empirical observational data. Since all models and physical laws are approximations to the true physics, a numerically constructed law or governing equation, which albeit may not possess a concise analytical form, could serve as a good model, so long as it is reasonably accurate.

Early work in this direction includes [3, 35], where symbolic regression was used to select physical laws and determine the structure of the underlying dynamical system. More recently, a sparsity-promoting method [7] was proposed to discover governing equations from noisy measurement data. It was based on a sparse regression method ([41]) and a Pareto analysis to construct parsimonious governing equations from a combinatorially large set of candidate models. This methodology was further applied to recovering partial differential equations [31, 34], as well as inferring structures of network models [25]. Conditions were provided in [42] for recovering the governing equations from possibly highly corrupted measurement data, along with a sparse convex optimization method to recover the coefficients of the underlying equations in the space of multivariate polynomials. A nonconvex sparse regression approach was proposed in [32] for selection and identification of a dynamical system directly from noisy data. A model selection approach was developed in [26] for dynamical systems via sparse regression and information criteria. Combining delay embedding and Koopman theory, a data-driven method was designed in [6] to decompose chaotic systems as an intermittently forced linear system. Very recently, three sampling strategies were developed in [33] for learning quadratic high-dimensional differential equations from under-sampled data. There are other techniques that address various aspects of dynamical system discovery problem. These include, but are not limited to, methods to reconstruct the equations of motion from a data series [12], artificial neural networks [19], nonlinear regression [43], equation-free modeling [21], normal form identification [24], empirical dynamic modeling [40, 48], nonlinear Laplacian spectral analysis [18], modeling emergent behavior [30], and automated inference of dynamics [36, 14, 15]. The readers are also referred to the introduction in [5].

The focus of this paper is on the study of local recovery/approximation of general first-order nonlinear ordinary differential/algebraic equations from the time-varying measurement data. The major new features of the proposed numerical framework include the following. First, our method seeks “accurate approximation” to the underlying governing equations. This is different from many of the existing studies, where “exact recovery” of the equations is the goal. The justification for “approximation”, rather than “recovery”, is that all governing equations will eventually be discretized and approximated by certain numerical schemes during implementation. As far as computations are concerned, all numerical simulations are conducted using “approximated” equations. Therefore, as long as the errors in the equation approximation/recovery step is small, comparable to the discretization errors, the approximated equations can be used as sufficient model for simulations. By relaxing the requirement from “exact recovery” to “accurate approximation”, we are able to utilize many standard approximation methods using standard basis functions such as polynomials, which are used in this paper. Consequently, one does not need to choose a (combinatorially) large dictionary in order to include all possible terms in the equations. Since most of the governing equations we encounter in practice have relatively simple and smooth forms, polynomial approximation can be highly accurate at moderately low degrees. Secondly, we propose the use of multiple, in fact a large number of, trajectory data of very short span, instead of using a single long trajectory data as in many existing studies. The length of each burst of the trajectory can be exceptionally short – it needs to be merely long enough to enable accurate estimate of the time derivatives. For noiseless trajectories, the length can as little as 2, for first-order estimation of the time derivative, or 3 for second-order estimation, and so on. We do require that the number of such short bursts of trajectories to be large so that they provide “good” samples in the domain of interest. Upon doing so, we convert the equation recovery problem into a function approximation problem with over-sampled data. A large number of well established approximation algorithms can then be employed. (The readers are referred to earlier books such as [11, 16, 27, 29] for the classical results and methods in function approximation theory.) In this paper, we discuss the results obtained by using least squares method, which is one of the most established methods; regression method, which is capable of removing data corruption ([9, 38]); and a matrix-free sequential algorithm ([39, 45, 46]), which is suitable for extraordinarily large data sets or stream data. We remark that the use of a very lagre number of trajectory bursts does not necessarily induce large simulation cost, as each burst is of very short length. In fact, in many cases the total amount of data in the large number of bursts is less than that in a single long trajectory. We also remark that the idea of using multiple bursts of trajectories has also be discussed in [33]. The major difference between [33] and this work is that the goal in [33] is for exact recovery of the equation of quadratic form. The bursts are random selected based on the theory of sparse recovery (compressed sensing). In this paper, the bursts are selected using approximation theory, as we only seek to approximate the equations and not recover them. And the equations can be of arbitrary form, beyond polynomials and including algebraic constraints. The choices of the bursts can be random or deterministic, as long as they enable one to obtain accurate approximation.

This paper is organized as follows. After the basis problem setup in Section 2, we present the main method in Section 3, which includes the general formulation, basis selections, trajectory selections, choices of approximation algorithms, and an error bound. We then present, in Section 4, an extensive set of numerical examples, covering a wide variety of linear/nonlinear differential/algebraic equations, to demonstrate the effectiveness of the proposed algorithms.

2. Problem Setup

We consider an autonomous system in the form of

where u = (1, are state variables, the initial time, and f defines the evolution of the states. Throughout this paper we call ) a trajectory, whose dependence on the initial state is explicitly shown. The right-hand-side f is the “true” and unknown dynamics. Our goal is to construct an accurate representation of f using measurement data of the state variables u.

2.1. Domain-of-Interest

While the state variables u may reside in the entire space , we restrict our equation recovery to a bounded domain . We refer to D as the domain-of-interest, which stands for the range of the state variables where one focuses on the underlying problem. Depending on the (unknown) dynamical system and one’s interest in studying it, the definition of the domain-of-interest varies.

be a closed interval and represent a range of the state variable , 1 , where one wishes to understand the behavior of the system. We then define D as the product of these intervals, i.e.,

Effectively this creates a hypercube in d dimensions. We then equip it with a measure , which is absolutely continuous and (without loss of generality) satisfies= 1. We then invoke the standard definition of inner product

and its induced norm .

Remark 1. Upon defining the domain-of-interest, we have restricted our equation recovery effort to a “local” manner. This is because recovering the equation in the entire space , albeit desirable, may be highly challenging, when the underlying dynamical system exhibits complex behaviors that are different in different parts of the space. The use of the domain-of-interest allows one to focus on a region of the domain where the problem can be well behaved and f possesses a regular/smooth form, which would be highly advantageous for the recovery effort. The size of D depends on the underlying problem and is not necessarily small.

2.2. Data

Let be a sequence of time instances. We assume that the true state variables u can be observed on these time instances and denote

the observed data for the state trajectory, where stands for observation errors. The trajectory data can be from different sequences of trajectories, which are originated by different initial states. Let 1, be a set of different initial states. We use the following shorthanded notation to denote the trajectory data from the m-th initial condition,

where is the number of intervals in the m-th trajectory data. The collection of all these data shall be used to estimate the true state equation.

2.3. Extension to Differential-Algebraic Equations

It is worth noting that the method discussed in this paper is applicable to semi-explicit differential algebraic equations (DAEs), i.e., DAEs with index 1, in the following form,

where u and v are the state variables and G represents algebraic constraints. Let us assume that the implicit function theorem holds such that the constraint (5b) induces the following explicit relation

Then, by letting f(u) = F(u, g(u)), the DAEs can be written in the same form as the autonomous system (1).

3. Main Numerical Algorithm

In this section we discuss numerical aspects of our equation approximation/recovery method. We first present the general formulation, which is rather similar to many of the existing work (e.g., [3, 7, 32]). The details of the numerical algorithm, which is different and new, are then presented, along with an error analysis.

3.1. General Formulation

Let {x} be a collection of trajectory data, in the form of (4). We assume that the time derivatives of the trajectory are also available at some time instances , either produced by the measurement procedure or approximated numerically using {x}. Let

be the derivative data, where ) stands for the errors or noises. This creates a set of data pairings at the time instances , where both the state variables and their derivatives are available. The collection of such pairings are denoted as

where the total number of data pairings is

We then choose a (large) set of dictionary Φ = to represent f in the governing equation (1). For the vector function f = (), its approximation is typically conducted component-wise. That is, for each = 1, . . . , d, we seek

Here the functions in the dictionary should contain all the possible functions one wishes to incorporate, for the accurate approximation of f. Some examples are

Let A be a () model matrix

where () represents the i-th data pairing in (7). The equation recovery problem is then cast into an approxi-

mation problem

where ˙is the vector of -th component of all available derivatives in the pairings (7), and = (is

as an approximation of the true system (1).

It should be noted that converting the equation recovery problem into an approximation problem in the form of (10) is a rather standard framework in the existing studies. In most of the existing studies, one typically chooses the dictionary (9) to be very large (often combinatorially large) and then use a certain sparsity-promoting method to construct a parsimonious solution, e.g., minimization for underdetermined system. By doing so, if the terms in the true equations are already included in the dictionary (9), it is then possible to exactly recover the form of the governing equations.

Hereafter, we will present the detail of our numerical approach, which is different from most of the existing approaches.

3.2. Equation Approximation using Standard Basis

Contrary to most of the existing studies, our method here does not utilize a combinatorially large dictionary to include all possible terms in the equations. Instead, we propose to use a “standard” set of basis functions, particularly polynomials in term of x, to approximate f directly. One certainly is not restricted to polynomials and can use other basis such as radial basis functions. Here we choose polynomials merely because it is one of the most widely used basis for function approximation. By doing so, we give up the idea of “exact recovery” of f, as exact recovery is not possible unless the true governing equations are strictly of polynomial form (which in fact is not uncommon). Instead, our approach focuses on obtaining an approximation of f with sufficiently high accuracy. We remark that this is satisfactory for practical computations, as the exact governing equations (even if known) are often discretized and approximated by a certain numerical scheme after all. So long as the errors in the equation recovery/approximation step is sufficiently small, compared to the discretization errors, the approximated governing equations can be as effective as the true model. Even though our method will not recover the true equations exactly in most cases, we shall continue to loosely use the term “equation recovery” hereafter, interchangeable with “equation approximation”.

To construct the approximation ˜f, we confine ourself to a finite dimensional polynomial space ) with dim 1. Without loss of generality, we take , and let the subspace V be Π, the linear subspace of polynomials of degree up to 1. That is,

where i = () is multi-index with . Immediately, we have

Let be a basis of V . We then seek an approximation ˜in a component-wise manner. That is, for the -th component, we seek

Upon adopting the following vector notations

Our goal is to compute the expansion coefficients in = 1, . . . , d. Later in Section 3.4 we will introduce several algorithms to compute the coefficients .

Unless the form of the true system (1) is of polynomial, we shall not obtain exact recovery of the equation. For most physical systems, their governing equations (1) consist of smooth functions, we then expect polynomials to produce highly accurate approximations with low or modest degrees. On the other hand, if one possesses sufficient prior knowledge of the underlying system and is certain that some particular terms should be in the governing equations, then those terms can be incorporated in the dictionary, in addition to the standard polynomial basis.

3.3. Selection of Trajectories

Upon setting the basis for the approximation, it is important to “select” the data. Since the problem is now casted into an approximation problem in , many approximation theories can be relied upon. One of the key facts from approximation theory is that, in order to construct accurate function approximation, data should fill up the underlying domain in a systematic manner. This implies that we should enforce the data pairings (7) to “spread” out the domain-of-interest D. Since the trajectories of many dynamical systems often evolve on certain low-dimensional manifolds, using a single trajectory, no matter how long it is, will unlikely to produce a set of data pairings for accurate approximation, except for some special underlying systems such as the chaotic systems [42]. This is the reason why we propose to use multiple number of short trajectories. Specifically, the trajectories are started from different initial conditions. The length of each trajectories may be very short — just long enough to allow one to estimate the time derivatives. More specifically, we conduct the following steps.

• Let 1 be the number of trajectories, and

We propose to use a large number of short bursts of trajectories, i.e.,

The length of each trajectory, (+1), is very small – it should be merely long enough to allow accurate estimation of the time derivatives. In case the time derivatives are available by other means, e.g., direct measurement, then can be set to 0.

There are several approaches to choose different initial states from the region D. They include but are not limited to

• Sampling the initial states from some random distribution defined on D, e.g., the uniform distribution.

• Sampling the initial states by quasi-Monte Carlo sequence, such as the Sobol or Halton sequence, or other lattice rules.

• If D is a hypercube, the initial state set can be taken as the set of quadrature points (e.g., the Gaussian points) or its random subset [49]. One can also sample the quadrature points according some discrete measure [45].

3.4. Numerical Strategies

We now discuss several key aspects of the numerical implementation of the method.

3.4.1. Computation of time derivatives

In most of the realistic situations, only the trajectory data are available, and the time derivatives ) need to be approximated numerically. Given a set of trajectory data , we seek to numerically estimate the approximate time derivative ) at certain time instances . This topic belongs to the problem of numerical differentiation of discrete data.

If the data samples are noiseless, the velocity can be computed by a high-order finite difference method, e.g., a second-order central approximation with = 2 and is

However, when data contain noise, estimation of derivatives can be challenging. In this case, if one uses the standard finite difference methods, then the resulting noise level in the first derivatives ) scales as . (See [44] for more analysis.) Several approaches were developed for numerical differentiation of noisy data. They include least squares approximation methods [22, 44], Tikhonov regularization [13], total variation regularization [10], Knowles and Wallace variational method [23], etc. For an overview, see [22].

In this paper, we use the least squares approach to denoise the derivative ). The approach consists of constructing a least squares approximation of the trajectory, followed by an analytical differentiation of the approximation to obtain the derivatives. With the m-th burst of trajectory data, m = 1, . . . , M, we first construct a polynomial p(t) = as an approximation to the trajectory ), by using the least squares method with 1 . We then estimate the derivative as ) = .

We remark that other methods for derivative computation with denoising effect can be certainly used. We found the least squares method to be fairly robust and easy to implement. When the noise level in the data become too large to obtain accurate derivatives, our method would become unreliable. This is, however, a universal challenge for all methods for equation recovery.

We also remark that even though the number of trajectories can be very large, i.e., 1. Since each trajectory is very short with (1), m = 1, . . . , M, the total amount of data for the state variables scales as O(M). In many cases this is fairly small.

3.4.2. Recovery algorithms: Regression methods

The regression methods seek the expansion coefficients by solving the following minimization problems

where the norm can be taken as the vector -norm , which corresponds to the well known least squares method and can handle noisy data. It can also be taken as the vector -norm, which is the least absolute deviation (LAD) method and can eliminate potential corruption errors in data (cf. [9, 38]). If one would like to seek the sparse coefficients , the LASSO (least absolute shrinkage and selection operator) method [41] may be used,

where is a parameter to be chosen. Since our method here seeks accurate approximation to the governing equations, rather than exact recovery, we focus on the use of least squares and LAD.

Remark 2. The methods discussed here are based on component-wise approximation of the state variable u. That is, each component of the states is approximated independently. It is possible to seek joint approximation of all components of u simultaneously. A joint approximation may offer computational advantages and yet introduce new challenges. We shall leave this to a separate study.

3.4.3. Recovery algorithms: Matrix-free method

When the data set is exceptionally large (1), or unlimited as in stream data, matrix-free method such as the sequential approximation (SA) method (cf. [39, 45, 46, 37]), can be more effective.

The SA method is of iterative nature. It aims at conducting approximation using only one data point at each step. By starting from arbitrary choices of the initial approximation, which corresponds to an arbitrary choice of the coefficient vectors = 0 step, the iteration of the coefficient vector at 1 step is

where 0 are parameters chosen according to the noise level in the k-th piece of data, and ˙denotes the -th component of ˙x. For noiseless data, one can set = 0.

The scheme is derived by minimizing a cost function that consists of data mismatch at the current arriving data point and the difference in the polynomial approximation between two consecutive steps. The method is particularly suitable for stream data arriving sequentially in time, as it uses only one data point at each step and does not require storage or knowledge of the entire data set. Moroever, the scheme uses only vector operations and does not involve any matrices. For more details of the method, including convergence and error analysis, see [39, 45, 46, 37].

3.5. Error analysis

Assuming the aforementioned numerical algorithms can construct a sufficiently accurate approximation of f with a error bound , we can estimate the error in solution of the recovered system (11).

Let x(t) be the exact solution of the recovered system (11) and u(t) the exact solutions of true system (1). Assume that ) = , we then have the following error bound for the difference between x(t) and u(t).

Theorem 1. Assume f is Lipschitz continuous on D,

If , for , then

where the operator denotes the component-wise projection on to V , and K(x) denotes the “diagonal” of the reproducing kernel of V . With any orthogonal basis of V , the kernel K can be expressed as K(x) = (x).

The proof can be found in Appendix A. We remark that the error bound here is rather crude, as it addresses a general system of equation. More refined error bounds can be derived when one has more knowledge of the underlying system of equations.

4. Numerical Examples

In this section we present numerical examples to verify the performance of the proposed methods for recovering the ordinary differential equations (ODEs) and semi-explicit differential algebraic equations (DAEs).

In all the test cases, we assume that there is no prior knowledge about the governing equations. Without loss of generality, we assume that . Unless otherwise stated, ∆t is taken as 0.005, the initial states are sampled by the uniform distribution over the region D, and the basis is taken as . We employ a variety of test cases of both linear and nonlinear systems of ODEs. The examples include relatively simple and well understood ODE systems from the textbook [4], as well as more complex nonlinear systems of coupled ODE and DAEs. Moreover, we also introduce corruption errors in the trajectory data, in addition to the standard random noises. We examine not only the errors in approximating the exact equations, but also the errors in the predictions made by the recovered equations.

Our primary focus is to demonstrate the flexibility of using polynomial approximation for accurate equation approximation, rather than exact recovery, as well as the effect of using a large number of short bursts of trajectories. For the actual computation of the approximation, we primarily use the standard least squares method. When data corruption exists, we also use LAD, i.e., regression, and employed the software package -Magic [8]. We also include an example with a very large data set and use the sequential method to demonstrate the flexibility of the matrix-free method.

4.1. Linear ODEs

We first present a set of numerical results for simple linear ODEs with d = 2, to examine the performance of the methods in recovering the structure of ODEs around different types of critical points. All examples are textbook problems from [4].

Example 1. Six linear ODE systems are considered here, and their structures cover all the elementary critical points. The description of the systems and the corresponding numerical results are displayed in Table 4.1. In all the tests, n = 1 is used. The first and second tests use M = 10 sets of noiseless trajectory data. The third and forth tests use M = 10 sets of noisy trajectory data, where the noise follows uniform distribution in [01].

The last two tests considers M = 20 sequences of trajectory data, with two uncertain (randomly chosen) sequences (10% of the total sequences) have large corruption errors following the normal distribution N(0.5, 1). It is seen from Table 4.1 that the proposed methods are able to recover the correct structure of the ODE systems from few measurement data. The values of approximate coefficients in the learned governing equations are very close to the true values. We now discuss the need for using a large number of short bursts of trajectory data. Let us consider the fourth linear ODE system in Table 4.1, which is the nodal sink system

Assume that the observed data are from a single trajectory

originated from the initial state . If belongs to the line Γ = {u = () : 1 = 0}, then all the data on the trajectory belong to Γ. It is straightforward to see that the solution of the following linear

system

can exactly fit the data on Γ, for arbitrary . This implies that recovering (27) becomes an ill-posed problem, as there are infinite number of linear systems (28) that correspond to such a single trajectory data. From the point of view of function approximation, this is not surprising, as a single trajectory of a given system often falls onto a certain low dimensional manifold (unless the system is chaotic). Approximating an unknown function using data from a low-dimensional manifold is highly undesirable and can be ill-posed. The use of a large number of short bursts of trajectory effectively eliminates this potential difficulty.

4.2. Nonlinear ODEs

We then test the proposed methods for recovering nonlinear ODE systems exhibiting more complex dynamics. In general, the nonlinear ODEs do not have explicit analytical solutions, we therefore generate the measurement data by numerically solving the ODEs using the forth-order explicit Runge-Kutta method with a small time step-size.

Example 2. We first consider the undamped Duffing equation

where is a small positive number and specified as 10. Taking and = ˙u gives

The region D is specified as [0, 2] 1]. From M = 30 bursts of noiseless trajectory data with = 2, we try to learn the governing system using different regression methods with n = 3. The errors of approximate expansion coefficients in the learned systems are displayed in Fig. 1. As we can see, the coefficients are recovered with high accuracy by and regressions.

Example 3. We consider

whose qualitative behavior was studied in [4]. We take D = [2]3], which contains the four critical points, an unstable node (0, 0), two asymptotically stable nodes (1, 0) and (0, 2), and a saddle point (0.5, 0.5). Suppose we have M = 30 sets of trajectory data with = 30, and all the data contain i.i.d. random noises following uniform

Table 1: Example 1: Local recovery of linear ODEs

Figure 1: Example 2: Results for recovering the system (29). The symbols “” represent the results of regression methods, respectively. The coefficient errors which are not shown in the figures are zero errors.

distribution in [01]. Learning from these data, the least squares regression method with n = 2 produced

the following approximate system

Although in this case the data contain noises at a relatively large level, the proposed method can still recover the coefficients with errors less than 4 10. In Figure 2, the phase portraits of the original and recovered systems are plotted to provide a qualitative comparison between the true and the learned dynamics. Good agreements are observed in the structures, and the four critical points are correctly detected. Figure 3 shows the solutions of the original and recovered systems at the initial state (1.25, 1.75), as well as the errors in the learned solution.

An important feature of this example is that we assume the state variables in the unit disk + are not accessible. Therefore we have data only from outside this region. This rather arbitrary choice of accessibility of data is to mimic the practical situation where data collection may be subject to certain restrictions. Also, we introduce relatively large corruption errors to the data. Therefore, we have a nonstandard region D = [2], where data are available and corrupted.

where 10% of the data contain relatively large corruption errors following i.i.d. normal distribution N(0.5, 1). (Note the relatively large bias in the mean value.) We then use -regression to learn the governing equations from these data, and obtain

We see that the -regression can effectively eliminate the sparse corruption errors, and the coefficients are recovered with high accuracy. The phase portraits in Figure 4 further exhibit the accurate dynamics of the learned system (33).

Figure 2: Example 3: Phase portraits for the true system (30) (left) and the learned system (31) (right).

Figure 3: Example 3: Validation of the recovered system (31) with initial state (1.25, 1.75). Left: solutions. Right: error . The horizontal axis denotes the time.

We also validate the solution of the recovered system (33) with an arbitrarily chosen initial state (325, 1.874), see Figure 5. The solution x(t) agrees well with the reference solution u(t) obtained from the exact equation. The increase of the error over time is expected, as in any numerical solvers for time dependent problems.

Figure 4: Example 4: Phase portraits for the true system (32) (left) and the learned system (33) (right).

Figure 5: Example 4: Validation of the recovered system (33) with initial state (874). Left: solutions. Right: error . The horizontal axis denotes the time.

We now consider data only with standard random noises and study the effect of noise level on the accuracy of recovered equations. Suppose we have M = 40 sets of trajectory data with = 10, and all the data contain i.i.d. random noises following uniform distribution in []. Here the parameter represents the noise level, and = 0 corresponds to noiseless data. The errors of the expansion coefficients in the learned systems are displayed in Fig. 6, for different noise levels at = 0, 0.005, 0.01, 0.05, 0.1, 0.2. The vector 2-norms of the coefficient errors are shown in Fig. 7. It is clearly seen that the errors become larger at larger noise level in the data, which is not surprising.

Figure 6: Example 4: Errors in each terms of the recovered equation at different noise levels.

Figure 7: Example 4: The norm of the errors in the recovered coefficients of the equations.

Example 5. We now consider a damped pendulum problem. Let l be the length of the pendulum, be the damping constant, and g be the gravitational constant. Then, for the unit mass, the governing equation is well

where is the angle from its equilibrium vertical position. By defining and = ˙, we rewrite the equation as

We set l = 1.1 and = 0.22, and use a large number of short trajectory data (M = 200= 2) to recover the equation. The initial state of each trajectory is randomly and uniformly chosen from the region D = [[]. The least squares -regression is used, with polynomial approximation at different degrees n. Figure 8 shows the phase portraits for the system (34) and the learned systems. As n increases, the recovered dynamics becomes more accurate and can capture the true structures of the solution. More detailed examination of the learned system is shown in Figure 9, with polynomial degree n = 6. We observe good agreements between the solution of the learned system x(t) and true solution u(t), at an arbitrarily chosen initial state (193876). The error evolution over time remains bounded at 10level.

4.3. Nonlinear DAEs

We now study semi-explicit DAEs (5) from measurement data. The learned approximating DAEs are denoted by

where x(t) and y(t) are approximations to u(t) and v(t), respectively. In the tests, the measurement data are first generated by numerically solving the true DAEs using the forth-order explicit Runge-Kutta method with a small time step-size. Possible noises or corruption errors are then added to the trajectories to generate our synthetic data.

Example 6. The first DAEs model we consider is a model for nonlinear electric network [28], defined as

where denotes the node voltage, and and are branch currents. Following [28], the physical parameters are specified as C = 10= 10= 1, 1 and = 0.25. In our test, D = [2] 2], M = 200 sequences of trajectory data with ∆t = 10and = 19 are generated by numerically solving (36) and then adding random noises. The added noises in and follow uniform distribution on [0001, 0.0001], while the noises in and follow normal distribution N(0, 0.001) and N(0, 0.005) respectively. We then use the -regression with n = 8 to learn the governing equations from these data. The phase portraits for the system ˙u = f(u) and the learned system ˙x = ˜f(x) are compared in Figure 10. It is observed that the proposed method accurately reproduces the dynamics and detects the correct features. Comparison the solutions between the recovered DAEs and the true system is shown Figure 11, using an arbitrary choice of the initial state (0, 0.1). The evolution of the numerical errors in the solution of the learned system are shown in Figure 12.

Example 7. We now consider a model for a genetic toggle switch in Escherichia coli, which was proposed in [17] and numerically studied in [47]. It is composed of two repressors and two constitutive promoters, where each promoter is inhibited by the repressor that is transcribed by the opposing promoter. Details of experimental measurement

Figure 8: Example 5: Phase portraits for the system (34) and its approximate systems learned from measurement data with different n.

Figure 9: Example 5: Validation of the recovered system (n = 6) with initial state (876). Left: solutions. Right: error . The horizontal axis denotes the time.

Figure 10: Example 6: Phase portraits for the system ˙u = f(u) (left) and the learned system ˙

Figure 11: Example 6: Validation of the recovered system with initial state (0, 0.1). The solid lines denote the solution of the exact system, and the symbols “” represent the solution of the learned system. The horizontal axis denotes the time.

Figure 12: Example 6: Same as Figure 11 except for the errors (left) and

can be found in [10], where the following model was also constructed,

In the last equation, we add a small perturbation term sin to increase the nonlinearity. Here, and are the concentration of the two repressors, respectively; and are the effective rates of synthesis of the repressors; and represent cooperativity of repression of the two promoters, respectively; and [IPTG] is the concentration of IPTG, the chemical compound that induces the switch. The values of these parameters are taken as = 156.25, = 15= 2= 1, = 2.0015. In the test, D = [0, 20]and we generate the measurement data by numerically solving (37) with [IPTG] = 10and = 0.01. Suppose we have M = 500 bursts of trajectory data with = 2, and five (5) arbitrarily chosen bursts of data contain corruption errors following normal distribution N(1, 1). We employ the -regression approach with n = 16 order of nominalized Legendre polynomial basis. The phase portraits for the system ˙u = f(u) and the learned system ˙x = ˜f(x) are plotted in Figure 13. As we can see, the learned structures agree well with the ideal structures everywhere in D. This demonstrates the ability of our method in finding the underlying governing equations. The accuracy of the learned model is further examined in Figures 14 and 15.

Figure 13: Example 7: Phase portraits for the system ˙u = f(u) (left) and the learned system ˙

Example 8. Our last example is a kinetic model of an isothermal batch reactor system, which was proposed in [2] and studied in [1, 20], among other literature. It was mentioned in [1] that the example in its original form was given by the Dow Chemical Company as a challenge test for parameter estimation methods. More details about this model can be found in [1]. The model consists of six differential mass balance equations, an electroneutrality

Figure 14: Example 7: Validation of the recovered system with initial state (19.5, 16.5). The horizontal axis denotes the time.

Figure 15: Example 7: Same as Figure 14 except for the errors (left) and

In our test, the values of the parameters are specified as = 25.1911, = 43.1042, = 1.1904 10, = = 2.575 10= 4.876, and = 1.7884 10, and [Q] = 0.0131. Assume that one can collect many sets of the time-series measurement data of u and v starting from any initial states D = [0[6[0[0[0[0, 0.02]. We collect 500, 000 bursts of trajectory data with ∆t = 10= 2 and sampled from the tensorized Chebyshev distribution on D. The data set is big as to enable possible accurate discovery of the complicated reaction process. Handling such a (relatively) large-scale problem with the standard regression algorithms can be challenging, primarily due to the large size of matrices. Without using more sophisticated regression methods to deal with large matrices (which is a topic outside our scope), we employ the sequential approximation (SA) algorithm to demonstrate its capability of handling large data sets. We employ the nominalized Legendre polynomial basis of up to order n = 12. This induces N = 18, 564 basis functions. The convergence result is examined by the errors in the function approximation ˜in the k-th step iteration compared with the exact functions f, g. To evaluate the errors we independently draw a set of 20,000 samples uniformly and compute the difference between ˜and f, g at these sampling points. The vector 2-norm is then used to report the difference. The relative errors in each iteration of SA are displayed in Figure 16 with respect to the number of iterations k. The SA method uses one piece of data at a time, and exhibits good effectiveness in this case. As the iteration continues, more data are used and more accurate approximations are obtained, yielding more accurate models. The exponential type convergence of the errors is observed. This is consistent with the theoretical analysis in [39, 46, 37]. We also validate the finial system with the initial state = (1.5776, 8.32, 0, 0, 0, 0.0142), and observe good agreements between the solutions of the recovered system and the exact system in Figure 17.

Figure 16: Example 8: Convergence history of the sequential approximation.

Figure 17: Example 8: Validation of the learned system. The solid lines denote the solution of the exact system, and the symbols “” represent the solution of the learned system. The horizontal axis denotes the time.

5. Conclusion

In this paper we studied several effective numerical algorithms for data-driven discovery of nonlinear differen-tial/algebraic equations. We proposed to use the standard basis functions such as polynomials to construct accurate approximation of the true governing equations, rather than the exact recovery. We also discussed the importance of using a (large) number of short bursts of trajectory data, rather than data from a single (long) trajectory. In conjunction with some standard approximation algorithms, e.g., least squares method, the overall method can produce highly accurate approximation of the true governing equations. Using an extensive set of test examples, ranging from textbook examples to more complicated practical examples, we demonstrated that the proposed method can be highly effective in many situations.

Appendix A. Proof of Theorem 1

We first introduce the following lemma.

Lemma 1. If ), we have

Proof. The triangular inequality gives

Let = (be any orthogonal basis of V , and

where the Cauchy-Schwarz inequality and the orthogonality of basis have been used. Taking -norm in the above inequality and using (A.2) give (A.1). Based on the above lemma, the proof of Theorem 1 is given as follows.

Proof. Integrating the systems (1) and (11) respectively gives

Subtracting (A.3) from (A.4) with ) = ), we obtain

where ) = is a convex function satisfying the Jensen’s inequality

This together with the triangular inequality yield

which implies (25). Applying Gr¨onwall’s inequality to (25) gives

Acknowledgment

This work was partially supported by AFOSR FA95501410022 and NSF DMS 1418771.

References References

[1] V. M. Becerra, P. D. Roberts, and G. W. Griffiths. Applying the extended Kalman filter to systems described by nonlinear differential-algebraic equations. Control Eng. Pract., 9(3):267–281, 2001.

[2] L. T. Biegler, J. J. Damiano, and G. E. Blau. Nonlinear parameter estimation: a case study comparison. AIChE Journal, 32(1):29–45, 1986.

[3] J. Bongard and H. Lipson. Automated reverse engineering of nonlinear dynamical systems. Proc. Natl. Acad. Sci. U.S.A., 104(24):9943–9948, 2007.

[5] S. Brunton, J. N. Kutz, and J. Proctor. Data-driven discovery of governing physical laws. SIAM News, 50:1, 2017.

[6] S. L. Brunton, B. W. Brunton, J. L. Proctor, Eurika Kaiser, and J. N. Kutz. Chaos as an intermittently forced linear system. Nature Communications, 8, 2017.

[7] S. L. Brunton, J. L. Proctor, and J. N. Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl. Acad. Sci. U.S.A., 113(15):3932–3937, 2016.

[8] E. Candes and J. Romberg. -MAGIC: Recovery of Sparse Signals via Convex Programming. www.acm.caltech.edu/l1magic/downloads/l1magic.pdf, 2005.

[9] E. J. Candes and T. Tao. Decoding by linear programming. IEEE Trans. Inform. Theory, 51(12):4203–4215, 2005.

[10] R. Chartrand. Numerical differentiation of noisy, nonsmooth data. ISRN Applied Mathematics, 2011, 2011.

[11] E.W. Cheney. Introduction to Approximation theory. McGraw-Hill, New York, 1966.

[12] J. P. Crutchfield and B. S. McNamara. Equations of motion from a data series. Complex Systems, 1(417-452):121, 1987.

[13] Jane Cullum. Numerical differentiation and regularization. SIAM J. Numer. Anal., 8(2):254–265, 1971.

[14] B. C. Daniels and I. Nemenman. Automated adaptive inference of phenomenological dynamical models. Nature Communications, 6, 2015.

[15] B. C. Daniels and I. Nemenman. Efficient inference of parsimonious phenomenological models of cellular dynamics using S-systems and alternating regression. PloS One, 10(3):e0119821, 2015.

[16] P.J. Davis. Interpolation and approximation. Dover, 1975.

[17] T. S. Gardner, C. R. Cantor, and J. J. Collins. Construction of a genetic toggle switch in escherichia coli. Nature, 403(6767):339, 2000.

[18] D. Giannakis and A. J. Majda. Nonlinear Laplacian spectral analysis for time series with intermittency and low-frequency variability. Proc. Natl. Acad. Sci. U.S.A., 109(7):2222–2227, 2012.

[19] R. Gonzalez-Garcia, R. Rico-Martinez, and I. G. Kevrekidis. Identification of distributed parameter systems: A neural net based approach. Comput. Chem. Eng., 22:S965–S968, 1998.

[20] S. C. Kadu, M. Bhushan, R. Gudi, and K. Roy. Modified unscented recursive nonlinear dynamic data recon- ciliation for constrained state estimation. J. Process Control, 20(4):525–537, 2010.

[21] I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidid, O. Runborg, C. Theodoropoulos, et al. Equation- free, coarse-grained multiscale computation: Enabling mocroscopic simulators to perform system-level analysis. Commun. Math. Sci., 1(4):715–762, 2003.

[22] I. Knowles and R. J. Renka. Methods for numerical differentiation of noisy data. Electronic Journal of Differential Equations, 21(2012):235–246, 2014.

[23] I. Knowles and R. Wallace. A variational method for numerical differentiation. Numer. Math., 70(1):91–110, 1995.

[24] A. J. Majda, C. Franzke, and D. Crommelin. Normal forms for reduced stochastic climate models. Proc. Natl. Acad. Sci. U.S.A., 106(10):3649–3653, 2009.

[25] N. M. Mangan, S. L. Brunton, J. L. Proctor, and J. N. Kutz. Inferring biological networks by sparse identifi- cation of nonlinear dynamics. IEEE Transactions on Molecular, Biological and Multi-Scale Communications, 2(1):52–63, 2016.

[26] N. M. Mangan, J. N. Kutz, S. L. Brunton, and J. L. Proctor. Model selection for dynamical systems via sparse regression and information criteria. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 473(2204), 2017.

[27] M.J.D. Powell. Approximation theory and methods. Cambridge University Press, 1981.

[28] R. Pulch. Polynomial chaos for semiexplicit differential algebraic equations of index 1. Int. J. Uncertain. Quantif., 3(1), 2013.

[29] T.J. Rivlin. An introduction to the approximation of functions. Dover publication Inc, 1969.

[30] A. J. Roberts. Model Emergent Dynamics in Complex Systems. SIAM, 2014.

[31] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz. Data-driven discovery of partial differential equations. Science Advances, 3(4):e1602614, 2017.

[32] H. Schaeffer and S. G. McCalla. Sparse model selection via integral terms. Phys. Rev. E, 96(2):023302, 2017. [33] H. Schaeffer, G. Tran, and R. Ward. Extracting sparse high-dimensional dynamics from limited data. arXiv preprint arXiv:1707.08528, 2017.

[34] Hayden Schaeffer. Learning partial differential equations via data discovery and sparse optimization. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 473(2197), 2017.

[35] M. Schmidt and H. Lipson. Distilling free-form natural laws from experimental data. Science, 324(5923):81–85, 2009.

[36] M. D. Schmidt, R. R. Vallabhajosyula, J. W. Jenkins, J. E. Hood, A. S. Soni, J. P. Wikswo, and H. Lipson. Automated refinement and inference of analytical models for metabolic networks. Physical Biology, 8(5):055011, 2011.

[37] Y. Shin, K. Wu, and D. Xiu. Sequential function approximation with noisy data. J. Comput. Phys., 371:363– 381, 2018.

[38] Y. Shin and D. Xiu. Correcting data corruption errors for multivariate function approximation. SIAM J. Sci. Comput., 38(4):A2492–A2511, 2016.

[39] Y. Shin and D. Xiu. A randomized algorithm for multivariate function approximation. SIAM J. Sci. Comput., 39(3):A983–A1002, 2017.

[40] G. Sugihara, R. May, H. Ye, C. Hsieh, E. Deyle, M. Fogarty, and S. Munch. Detecting causality in complex ecosystems. Science, 338(6106):496–500, 2012.

[41] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methodological), pages 267–288, 1996.

[42] G. Tran and R. Ward. Exact recovery of chaotic systems from highly corrupted data. Multiscale Model. Simul., 15(3):1108–1129, 2017.

[43] H. U. Voss, P. Kolodner, M. Abel, and J. Kurths. Amplitude equations from spatiotemporal binary-fluid convection data. Phys. Rev. Lett., 83(17):3422, 1999.

[44] J. Wagner, P. Mazurek, and R. Z. Morawski. Regularised differentiation of measurement data. In XXI IMEKO . Prague, Czech Republic, 2015.

[45] K. Wu, Y. Shin, and D. Xiu. A randomized tensor quadrature method for high dimensional polynomial approximation. SIAM J. Sci. Comput., 39:A1811–A1833, 2017.

[46] K. Wu and D. Xiu. Sequential function approximation on arbitrarily distributed point sets. J. Comput. Phys., 354:370–386, 2018.

[47] D. Xiu. Efficient collocational approach for parametric uncertainty analysis. Commun. Comput. Phys., 2(2):293–309, 2007.

[48] H. Ye, R. J. Beamish, S. M. Glaser, S. C. H. Grant, C. Hsieh, L. J. Richards, J. T. Schnute, and G. Sugihara. Equation-free mechanistic ecosystem forecasting using empirical dynamic modeling. Proc. Natl. Acad. Sci. U.S.A., 112(13):E1569–E1576, 2015.

[49] T. Zhou, A. Narayan, and D. Xiu. Weighted discrete least-squares polynomial approximation using randomized quadratures. J. Comput. Phys., 298:787–800, 2015.

designed for accessibility and to further open science