Data-driven PDE discovery with evolutionary approach

TL;DR

This paper introduces an evolutionary symbolic regression approach for discovering partial differential equations from data, overcoming limitations of sparse regression and tested on canonical PDEs with noise robustness analysis.

Contribution

It presents the EPDE method, a novel evolutionary approach for PDE discovery that offers greater flexibility in model form compared to traditional sparse regression techniques.

Findings

EPDE successfully discovers canonical PDEs from data.

The method demonstrates robustness to noisy data.

Fewer restrictions on PDE form compared to existing methods.

Abstract

The data-driven models allow one to define the model structure in cases when a priori information is not sufficient to build other types of models. The possible way to obtain physical interpretation is the data-driven differential equation discovery techniques. The existing methods of PDE (partial derivative equations) discovery are bound with the sparse regression. However, sparse regression is restricting the resulting model form, since the terms for PDE are defined before regression. The evolutionary approach described in the article has a symbolic regression as the background instead and thus has fewer restrictions on the PDE form. The evolutionary method of PDE discovery (EPDE) is described and tested on several canonical PDEs. The question of robustness is examined on a noised data example.

Figures

Tables

Table 1. Table 1: Equations used for algorithm validation.

Name	Equation
Burgers equation	$\frac{\partial u}{\partial t}=-u\frac{\partial u}{\partial x}+\mu \frac{{\partial }^{2}u}{\partial x^2}$
The wave equation	$\frac{{\partial }^{2}u}{\partial t^2}=\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial x^2}$
Korteweg-de Vries equation	$\frac{\partial u}{\partial t}+6u\frac{\partial u}{\partial x}+\frac{\partial u^3}{\partial x^3}=0$

Table 2. Table 2: Resulting coefficients for the equations

Term	1	u	$\frac{\partial u}{\partial t}$	$\frac{{\partial }^{2}u}{\partial t^2}$	$\frac{\partial u}{\partial x}$	$\frac{{\partial }^{2}u}{\partial x^2}$	$\frac{{\partial }^{2}u}{\partial t^2}\frac{{\partial }^{2}u}{\partial x^2}$	$u\frac{{\partial }^{2}u}{\partial t^2}$	$u\frac{\partial u}{\partial x}$	$\frac{{\partial }^{3}u}{\partial x^3}$
Wave Equation	0	0	0	1	0	-0.25006	0	0	0	0
Burgers	0	0	1	0	0	0.09985	0	0	-0.99986	0
KdV	0	0	1	0	0	0	0	0	-5.9992	-0.9997

Table 3. Table 3: Discovered structure of Burger‘s and KdV equations for different input matrix section.

Data part	Burger‘s correct	Burger‘s wrong	KdV correct	KdV wrong
0.9	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{2}u}{\partial x^2}$, $u\frac{{\partial }^{2}u}{\partial x}$	-	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{3}u}{\partial x^3}$, $u\frac{{\partial }^{2}u}{\partial x}$	-
0.8	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{2}u}{\partial x^2}$, $u\frac{{\partial }^{2}u}{\partial x}$	-	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{3}u}{\partial x^3}$, $u\frac{{\partial }^{2}u}{\partial x}$	-
0.7	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{2}u}{\partial x^2}$, $u\frac{{\partial }^{2}u}{\partial x}$	-	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{3}u}{\partial x^3}$, $u\frac{{\partial }^{2}u}{\partial x}$	-
0.6	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{2}u}{\partial x^2}$, $u\frac{{\partial }^{2}u}{\partial x}$	-	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{3}u}{\partial x^3}$, $u\frac{{\partial }^{2}u}{\partial x}$	-
0.5	$\frac{{\partial }^{2}u}{\partial x^2}$, $\frac{\partial u}{\partial t}$	$\frac{{\partial }^{2}u}{\partial x^2}\frac{\partial u}{\partial t}$	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{3}u}{\partial x^3}$, $u\frac{{\partial }^{2}u}{\partial x}$	-
0.4	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{2}u}{\partial x^2}$, $u\frac{{\partial }^{2}u}{\partial x}$	-	$\frac{\partial u}{\partial t}$, $\frac{{\partial }^{3}u}{\partial x^3}$, $u\frac{{\partial }^{2}u}{\partial x}$	-
0.3	$\frac{{\partial }^{2}u}{\partial x^2}$, $\frac{\partial u}{\partial t}$	$\frac{{\partial }^{2}u}{\partial x^2}\frac{\partial u}{\partial t}$	-	$\frac{{\partial }^{2}u}{\partial x^2}$, $\frac{{\partial }^{2}u}{\partial t^2}$
0.2	$\frac{{\partial }^{2}u}{\partial x^2}$, $\frac{\partial u}{\partial t}$	$\frac{{\partial }^{2}u}{\partial x^2}\frac{\partial u}{\partial t}$	-	$\frac{{\partial }^{2}u}{\partial x^2}$, $\frac{{\partial }^{2}u}{\partial t^2}$
0.1	$\frac{{\partial }^{2}u}{\partial x^2}$, $\frac{\partial u}{\partial t}$	$\frac{{\partial }^{2}u}{\partial x^2}\frac{\partial u}{\partial t}$	-	$\frac{{\partial }^{2}u}{\partial x^2}$, $\frac{{\partial }^{2}u}{\partial t^2}$

Equations

$$\begin{array}[]{cc}\frac{\partial u}{\partial t}=~{}u_{t}(x,t)=~{}\frac{u_{j}^{i+1}-u_{j}^{i-1}}{2\Delta t};~{}\frac{\partial u}{\partial x}=~{}u_{x}(x,t)=~{}\frac{u_{j+1}^{i}-u_{j-1}^{i}}{2\Delta x};\\ \frac{{{\partial}^{2}}u}{\partial{{t}^{2}}}=~{}\frac{u_{j}^{i+1}-2u_{j}^{i}+u_{j}^{i-1}}{\Delta{{t}^{2}}};~{}\frac{{{\partial}^{2}}u}{\partial{{x}^{2}}}=~{}\frac{u_{j+1}^{i}-2u_{j}^{i}+u_{j-1}^{i}}{\Delta{{x}^{2}}};\\ \frac{{{\partial}^{3}}u}{\partial{{x}^{3}}}=~{}\frac{u_{j+2}^{i}-2u_{j+1}^{i}+2u_{j-1}^{i}-u_{j-2}^{i}}{2\Delta{{x}^{3}}}\\ \end{array}$$

$${{f}_{1}}\left(t\right)=\left[\begin{array}[]{cc}1\\ \vdots\\ 1\\ \vdots\\ 1\\ \end{array}\right]\,,{{f}_{2}}\left(t\right)=\left[\begin{array}[]{cc}u\left(t,~{}{{x}_{0}}\right)\\ \vdots\\ u\left(t,~{}{{x}_{i}}\right)\\ \vdots\\ u\left(t,~{}{{x}_{n}}\right)\\ \end{array}\right]\,,{{f}_{3}}\left(t\right)=\left[\begin{array}[]{cc}{{u}_{x}}\left(t,~{}{{x}_{0}}\right)\\ \vdots\\ {{u}_{x}}\left(t,~{}{{x}_{i}}\right)\\ \vdots\\ {{u}_{x}}\left(t,~{}{{x}_{n}}\right)\\ \end{array}\right]\,,...$$

$$E=\left[\begin{array}[]{cc}{{f}_{1}^{N}}\left(t_{1}\right)\\ \vdots\\ {{f}_{1}^{N}}\left(t_{j}\right)\\ \vdots\\ {{f}_{1}^{N}}\left(t_{m}\right)\\ \end{array}\right]=\left[\begin{array}[]{cc}\frac{1}{\|{{f}_{1}}\left(t_{1}\right)\|_{2}}\\ \vdots\\ \frac{1}{\|{{f}_{1}}\left(t_{j}\right)\|_{2}}\\ \ \vdots\\ \frac{1}{\|{{f}_{1}}\left(t_{m}\right)\|_{2}}\\ \end{array}\right]\,,U=\left[\begin{array}[]{cc}{{f}_{2}^{N}}\left(t_{1}\right)\\ \vdots\\ {{f}_{2}^{N}}\left(t_{j}\right)\\ \vdots\\ {{f}_{2}^{N}}\left(t_{m}\right)\\ \end{array}\right]=\left[\begin{array}[]{cc}\frac{u(t_{1},x_{0})}{\|{{f}_{2}}\left(t_{1}\right)\|_{2}}\\ \vdots\\ \frac{u(t_{j},x_{0})}{\|{{f}_{2}}\left(t_{j}\right)\|_{2}}\\ \ \vdots\\ \frac{u(t_{m},x_{n})}{\|{{f}_{2}}\left(t_{m}\right)\|_{2}}\\ \end{array}\right]\,,...$$

$$F(j)=\left[\begin{array}[]{cc}(u^{\prime}(t_{1},x_{0})*u_{t}(t_{1},x_{0}))^{N}\\ \vdots\\ (u^{\prime}(t_{m},x_{n})*u_{t}(t_{m},x_{n}))^{N}\\ \end{array}\right]=U_{x}*U_{t}$$

$$Q=\sum _j\Vert w\cdot x_j-y_j{\Vert }_2^2+\lambda \Vert w{\Vert }_1\to \underset{w}{\min }$$

$$\underset{\alpha }{\min }\left(\sum _{k=0}^p\Vert F_k\alpha -F_{target,k}{\Vert }_2^2+\lambda \Vert \alpha {\Vert }_1\right)$$

$$S=\left[\begin{array}[]{cccc}|&|&|&\\ F(1)&F(2)&F(3)&\,...\\ |&|&|&\\ \end{array}\right]$$

$$f_{fitness}=\frac{1}{\Vert F\cdot \alpha -F_{target}{\Vert }_2}$$

$$\frac{{\partial }^{2}u}{\partial t^2}=\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial x^2}$$

$$\frac{\partial u}{\partial t}=-u\frac{\partial u}{\partial x}+\mu \frac{{\partial }^{2}u}{\partial x^2}$$

$$Q_{noise}=\frac{\Vert w_0-\tilde{w}{\Vert }_2}{\Vert w_0{\Vert }_2}\ast 100$$

$$E_{coeff}=\sqrt{\sum _{i=1}^N(w_i-w_{pred}{)}^2}$$

Full text

1 Introduction

Data-driven algorithms are usually considered as the source of models when the connection between the data samples is not known a priori. There are various data-driven models built on neural networks, regression, statistic, etc. As an example, deep neural networks models [1], recurrent-convolutional neural network models [9], statistical models [10], regression [6], combined evolutionary-based models [8] and other models and their combinations [5]. However, most of the existing models are unsuitable for interpretation. The expert usually is not able to determine, why the model made a particular decision.

On contrary, the models, where the connection between the data samples is based on the physical principles, are potentially interpreted. For example, one can built the physical law in form of function [13], ordinary differential equations system [3, 7]. However, most of the physical laws are written in form of the partial differential equation (PDE).

Partial differential equations (PDE) are able to represent a vast variety of processes, which occur in dynamic systems in nature and society. For example, the Navier-Stokes equation, that describes the flow of liquids or gasses, is often used to solve hydrodynamic problems. Maxwell’s equations, which represent electromagnetic fields, is comprised of four partial differential equations, characterizing interactions between changes of electric and magnetic fields and the effects of charges and electric currents.

However, derivation of these equations previously involved a primarily mathematical and physical approach to the problem, requiring a deep understanding of the phenomenon’s nature. This, consequently, reduces the ability to create models for systems that modern science has little to no knowledge about. For instance, there is no uniform theory that determines societal and economic dynamics. Furthermore, various unsolved problems can occur even is such explored studies as hydrodynamics and metocean science.

Data-driven algorithms are a solution for cases of systems, that we lack knowledge about. Nevertheless, in most cases raw observational data are available. The data-driven algorithms bring the ability to build the model for dynamical systems from time-series of data, received from in-field or laboratory observations. The development of the data-driven methodology of partial differential equations derivation, combined with recent advances in technologies of measurements and probing, brings new opportunities for studying of metocean dynamic systems.

Sparse regression is considered to be the main tool for selection of the leading terms of the differential equations [11]. The applied regularization is based on the addition of the L1 norm of the calculated weights to the least-square expression. One of the most popular methods used in PDE discovery is the least absolute shrinkage and selection operator (LASSO). The main feature of LASSO is the ability to mutate the loss function. Zero weights are chosen for terms, that poorly fit the input data, and, therefore, identify the structure of the PDE.

Previously, the problem of the discovery of the differential equation structure has been developed in a number of papers. From derivation of systems of equations, defining physical laws, by means of symbolic regression [4, 14] to study dynamic systems, that are represented by a system of partial differential equations [2, 12].

The methods of PDE derivation, used in previous papers, usually utilize regression over the set of the pre-determined terms, that are usually comprised of different polynomial combinations of derivatives and functions. This limitation provides only the discovery of equations, that have a corresponding structure. The method, presented in this paper is referred below as EPDE. It is based on a combination of sparse regression, that applies sparsity for the small set of potential terms. Sets, in turn, can have arbitrary form, obtained during the evolution process. Also, the proposed way of calculation of terms’ weights values includes the application of linear regression over the non-normalized data for selected terms.

The paper is organized as follows, Section 2 describes the problem of data-driven PDE discovery in details. Also, in Section 2 dataset for experiments is described. Section 3 describes the data-driven PDE discovery algorithm based on evolutionary optimization. Section 4 is dedicated to the analysis of algorithm precision, stability, and robustness. Section 5 concludes the paper.

2 Problem statement and data acquisition

The developed EPDE algorithm is aimed at the derivation of the dynamic systems governing equation by time series, containing information about the studied function (temperature, velocity, etc.). At first, the approach must be applied for test cases, including artificially created data, acquired from numerically solved equations. This simplification gives chances to exclude noise from data and, therefore, check the algorithm’s behavior independently from external conditions, such as faults of the measurement equipment. The ability to manually select parts of data and compare results gives opportunities for stability tests. Also, a noise of any magnitude can be added to data to investigate the reaction of the algorithm to it.

In this work, the algorithm was tested on the wave equation, Burgers and Korteweg-de Vries equations. They were solved numerically with the application of a finite-difference scheme to approximate time and spatial derivatives. For instance, the Crank-Nicolson method was utilized to solve the Burgers equation

From the acquired field of equation solution, its time and spatial derivatives are calculated in order to be utilized further in regression. These derivatives are calculated by the finite-difference method (Eq. 1) due to its simplicity and sufficient quality on noise-free data and small time and space steps.

[TABLE]

In cases, where the additional noise will be added or measurements are used, to numerically differentiate the solution, complex methods should be used. Tikhonov regularization and other types of regression are recommended due to their ability to remove noise so that the PDE discovery algorithm obtain cleaned data as input. After derivatives are obtained, it is possible to create vectors of spatial data for a specific time point:

[TABLE]

After that, the normalization of each of these time frames should be held. It can be done with the highest variable value for that time point, or by time frame’s L2-norm. Finally, data vectors are created by compositions of all time frames for the modeled period.

[TABLE]

Where $f_{j}(t)$ is the L2-normalized vector from Eq. 2. After normalized terms are found, the feature vectors F(j) are formed and written in terms Eq. 3 as, for example, the following product:

[TABLE]

On a balance, the data preparation step consists in representing data and their spatial and time derivatives in vectors in form Eq. 3. After these steps, features are collected in form Eq. 4 in order to perform the optimization procedure.

3 Algorithm description

The proposed algorithm includes two parts: the evolutionary algorithm that generates a small group of terms that are called individuals and sparse regression that allows choosing significant terms in the set of individuals. This process is shown schematically in Fig.1.

Such approach, in contrast to the existing algorithms, allows, on one hand, to generate more flexible space of terms for the sparse regression, and, on the other hand, to reduce the number of terms for the regression. Therefore, generation of the all possible terms is not required, as it is done within the evolutionary part. Sec. 3.1 contains the description of the sparse regression part of the algorithm. Sec. 3.2 contains the description of individual generation, crossover, and mutation.

3.1 Sparse regression

Sparse regression is one of the most common instruments for discovering partial differential equations due to its ability to create the form of terms that would compose the resulting PDEs. It is based least absolute shrinkage and selection operator that is commonly used in machine learning to prevent over-fitting by filtering insignificant and redundant features by penalizing coefficients before them. It is based on the addition of L1-norm of unknown coefficients into a function (Eq.5), used in the least-squares algorithm:

[TABLE]

In the sparse regression, w represents the searched sparse vector of weights between features, initialized by matrix x, and vector y is the learning target. $\lambda$ is a sparsity constant that is set before the learning process. To find values of the weights $\alpha$, that is representing the system’s PDE, it is possible to define the loss function (Eq. 6) in the following way, using the defined set of features and target vectors, created in the previous section:

[TABLE]

Where p is the number of features selected for the regression algorithm. This application of the regularized regression is not able to discover the true values of the weights due to the fact, that it uses normalized vectors of target and features. However, it is able to select leading ones with their sign. Due to the addition of L1-norm, the loss functions must be minimized, using optimization algorithms, that are able to work with non-differentiable functions, such as the subgradient method.

After the structure is found, the coefficients are defined with the non-normalized data, i.e. features are written back in the form (2) and regression is used to find the coefficients. Usually, in regression all possible combinations [12] of the feature vectors Eq. 4 are chosen for minimization problem Eq. 6. Thus, the optimization problem complexity grows exponentially as the maximal order of the derivative increases. With the evolutionary algorithm, described below, one can use multiple reduced optimization problems instead of full regression on a complete terms library. In existing algorithms as a target feature usually, the highest-order time-derivative is chosen

3.2 Evolutionary algorithm

The second element of the EPDE method is the evolutionary algorithm, that is aimed at the construction of the most complete set of terms. By its iterations, the evolutionary algorithm should be able to select and preserve the most appropriate elements of the resulting equation. Therefore, the sparse regression is done on every iteration of the evolutionary algorithm for every candidate in the population with a random selection of target among the set of terms.

To initiate the method, it is required to create a population of individuals, represented by chromosomes, where each gene represents a combination of functions and their derivatives. An evolutionary algorithm is able to vary the chromosomes in two ways: crossover, that represents the exchange of corresponding genes between two individuals, and mutation, which involves random alteration of chromosome’s genes. In the examined case, the mutation is held by the conversion of one term to the other randomly generated one.

Due to the specification of the task, every individual represents a specific case of the equation, having its own features matrix and the target vector. Vectors F(i), that compose the columns of the feature matrix S (Eq. 7), are created as a product of a randomly selected number of feature factors Eq. 4:

[TABLE]

It should be emphasized, that the number of feature vectors in Eq. 7 is the parameter of the evolutionary algorithm. The second remark is that, in contrast to the existing algorithms [11, 12, 7], the target feature is chosen randomly, whereas in the sparse-regression only cases time-derivative is used.

While mutation is usually applied to all individuals of the population, crossover occurs only between the most eligible of them. To select candidates for crossover, the fitness function should be introduced. For the task of partial differential equation derivation, it can be introduced by a norm of the difference between the target term and the expression with other ones i.e. regression error, calculated for all of the training data:

[TABLE]

A manner of the population’s participation in crossover should be defined before the initiation of an algorithm. In this research, the simplest way was adopted, executed by the selection of part the population part, that will breed. However, on some occasions that can result in stagnation of the evolution due to an occasional similarity between individuals with the highest fitness values, while the most eligible potential candidate requires gene from the fewer fit phenotypes. This limitation can be surpassed by the tournament selection. For that, a subset of the population shall be selected randomly and fitness tournaments must be held between them. The crossover procedure is schematically shown in Fig. 2.

In this case, the required diversity of the population will remain and, consequently, there are higher chances of better individual’s creation. The individuals, generated during the crossover, replace the least fit ones in order to keep the population quantity stable while leaving the phenotypes, that previously had good values of fitness function intact. The iterative process can be stopped at the moment when the growth of the fitness function decreases below the defined threshold. Finally, the individual with the highest fitness at the last iteration is considered to be the one, that represents the structure of the partial differential equation, describing the dynamic system.

After the sparse regression application, one more regression step is required. It is initialized over the set of terms, selected by non-zero weights in the previous step. In this step, non-normalized fields of variables are used as a feature and target vectors. This approach is uncommon in general machine learning due to its limitations on variables of different scale, where the algorithm is not able to properly generalize data and discover a contribution of each feature. However, in this particular case, the structure of the an equation, represented by weights of features, is already known, and these variables must be evaluated according to their scale.

The described algorithm allows one to reduce regression space. Additionally, it allows to theoretically find ordinary differential equation instead of the PDE since target feature is not restricted by the highest time-derivative. This is required for potential one-dimensional static problems ODE discovery.

4 Validation

To analyze the algorithm performance, it is necessary to make sure, that it has the following qualities: stability, approximation, and convergence. These qualities are dependent on each other, and to prove them, it is enough to check, if any two of them are fulfilled. Due to the reasons of convenience, in the research, stability, and convergence of the algorithm are studied. Convergence of the PDE deriving algorithm manifests in the improvement of the quality of the algorithm with the reduction of a step of the grid, from that is adopts data. Stability can be proved by addition of the noise to the input PDE solution and test, how this corruption affects the structure of the resulting equation.

The algorithm has proved to be capable of discovering partial differential equations structure and calculating the values of weights for the selected terms for all of the studied equations.

At first, the algorithm was tested on the wave equation (Eq. 9):

[TABLE]

where $c=2$ has the physical meaning of the sound velocity. It was solved on the spatio-temporal grid with 100x100 points. Time and spatial steps were set as 0.1. The equation was solved by the implicit method. The resulting coefficients obtained by the evolutionary PDE discovery algorithm are shown in Tab.2

To make sure that the resulting structure does not change with the reduction of studied space, the same algorithm was initiated over the parts of the matrix solution. For the case of the wave equation, due to its relative simplicity, the same results have been achieved on all studied sizes of the matrices: from the whole solution to 10$\%$ of it. These results show that such simple structures of dynamic systems can be easily detected by the utilized algorithm.

The Burgers’ equation (Eq. 10) presents an example of the more complex system to be studied:

[TABLE]

Where $\mu$ represents viscosity, which was set to a value of 0.1. The equation of solved on the grid of 256 x 256 spatial and time points correspondingly with steps of 16/256 and 10/256. To acquire data for the algorithm from the equation solution, the Crank-Nicolson method was used. the resulting coefficients are shown in Tab.2

The selected part of the solution matrix has influence over the results of regressions and, therefore, defines the equation’s structure. The results were tested on the parts of the matrix from 1.0 to 0.1 of its size. On the lesser sizes of the selected matrix part, especially for cases, when the selected part contains an only small part of the solution ridge, the algorithm can have difficulties, deriving wrong structures. Consequently, the incorrect set of terms prevents the calculation of their true weights during the second linear regression phase.

The results of the matrix division are presented in Table 3. Here it can be seen, that proposed algorithm fails to discover the structure of the dynamic system in cases of the low quantity of data: for sections, which are less than 0.5 of the solution and contain less than 130 points, the algorithm tends to create wrong terms of the equation. Therefore, to achieve the correct performance of the algorithm, data matrices of enough size shall be passed to it. This condition creates some limitations for the method application for cases of lack of data.

Similar results have been attained for the Korteweg-de Vries equation. It has been solved as shown in Table 2 - Table 2. We note that for Table 2 different number of points was taken in order to check the performance of the algorithm. For Burger’s equation, the 256x256 grid was taken whereas for the Korteweg-de Vries equation - 1024x1024 points.

The previously mention effect remained in this scenario: the algorithm only had issues in discovering the structure of the governing equation. For cases, when it succeeded, the true values of the weights were calculated correctly, even on minor parts of the equation’s solution matrix.

To check the evolutionary algorithm stability, the noise is added to the entire solution’s field. It is added from a normally distributed random variable with zero mean value and dispersion taken as the fraction of maximal value. As the invariant noise measure, Eq.11 is used.

[TABLE]

With $w_{0}$ in Eq.11 the initial (clean) solution field is designated, $\widetilde{w}$ is the field with noise added, $\|\cdot\|_{2}$ is the matrix’s Frobenius norm.

For comparison, we take the latest supplementary code for the article [11] from GitHub repository. Same Burger‘s equation solution field and same noise procedure implementation were taken. It should be noted, that we compare ”basic” versions of the algorithms. For the sparse regression more sophisticated derivative procedure and meta-parameter optimization for the regression algorithm could be implemented, which, definitely, increases the quality of both algorithms.

As the discovery precision metric coefficient root mean square error is taken as it is shown in Eq.12.

[TABLE]

In Eq. 12, $N$ is the number of terms taken for sparse regression, $w_{i}$ are the coefficients of clean (without noise) equation terms, $w_{pred}$ are the corresponding predicted coefficients.

Polynomial derivatives procedure was utilized, also for the sparse regression improved ridge regression with $\alpha=10^{-6}$ was taken. Comparison results for the Bruger‘s equation are shown in Fig.3.

After certain noise level limit $Q_{noise}\approx 0.11$ (it is seen in Fig.3 as the function jump discontinuity), the classical algorithm loses ability to discover the term $\frac{\partial^{2}u}{\partial x^{2}}$ without an additional regression tuning. However, it is still able to catch the leading term. In Fig.3 shown the maximal noise level range where the evolutionary approach is able to determine the equation structure.

We note that two other equations described in the article are not taken for the comparison due to the hyperbolic nature of the differential operators. Numerical solution of such equations requires specific methods that should be implemented for both algorithms. So, the comparison results, in this case, will show the quality of the implemented numerical schemes rather than the performance of the PDE discovery algorithms.

As seen the evolutionary approach allows one to extend the noise level which is allowed for all terms of the initial equation discovery. The term coefficients discovery precision is increased, which leads to more stable equation discovery and allows one to discover the equations in a more robust way.

5 Conclusions and discussion

In the paper evolutionary approach for PDE discovery is described. In contrast to the existing algorithms based on the regression on a complete terms library it has the following advantages:

• •

Regression is done on a reduced space, i.e. only a small amount of features is taken for the regression;

• •

More flexible features choice allows to obtain wider space of possible differential operators;

• •

No restriction on the target function is allowing to obtain more sophisticated forms of differential operators including ODEs;

The possible disadvantages could be:

• •

Possible extended computation time due to the stochastic process of the initial population initialization, population crossover and mutation;

• •

Additional procedures are required in order to maintain the robustness of the algorithm, i.e. in order to obtain the same model for the data of the same origin;

The proposed method can be considered as a base point for the data-driven PDE discovery with an evolutionary approach. In the article, the main stages of the methods are shown. Every stage could be improved, for example, a more sophisticated grid function differentiation method could be taken to increase precision and stability. Also, more advanced evolution methods could be used in order to increase computation efficiency and stability.

The questions of the stability concerning the different boundary conditions, more difficult differential operators are left out of the scope of the paper. As well as the question of real observational data model discovery. However, the main goal was to show the proof-of-concept algorithm and validate it on synthetic examples and partially compare it with the existing algorithms.