# The identification of continuous, spatiotemporal systems

thanks: published in: Phys. Rev. E 57 (1998) 2820  

## Abstract

We present a method for the identification of continuous, spatiotemporal dynamics from experimental data. We use a model in the form of a partial differential equation and formulate an optimization problem for its estimation from data. The solution is found as a multivariate nonlinear regression problem using the ACE–algorithm. The procedure is successfully applied to data, obtained by simulation of the Swift–Hohenberg equation. There are no restrictions on the dimensionality of the investigated system, allowing for the analysis of high-dimensional chaotic as well as transient dynamics. The demands on the experimental data are discussed as well as the sensitivity of the method towards noise.

The unstable dynamics observed in spatially extended systems attracted huge experimental and theoretical research activity in the last decades (see  and references therein). Progress has been achieved by describing the dynamics in the vicinity of bifurcations with the help of universal amplitude equations, vastly reducing the complexity of the involved models. Additionally, the research has concentrated on the classification of the observed instabilities and the resulting patterns, and the investigation of scaling laws and intermittency effects . For most investigations the models for spatiotemporal systems arise from mainly theoretical considerations and their validity is affirmed by the comparison with experimental measurements. Here we address the problem of finding a model, which describes the dynamics of an observed system, directly from experimental data. For systems which exhibit temporal low–dimensional chaotic motion this was accomplished with the help of nonlinear maps . Other general methods rely on some sort of mode–expansion and were successfully applied . A different approach consists in a nonparametric model identification as proposed for systems with a time delay feedback .

In this letter we extend that approach to the identification of the underlying evolution equation of spatially extended systems. At first we formulate an optimization problem for finding a model equation from the data. Then, we rewrite the equations in the form of a multivariate nonlinear regression problem. As a last step, we use a novel kind of numerical algorithm for solving the problem. Our approach does not include any parameter dependencies, rather those are delivered as a by-product. We discuss the identification of homogeneous and autonomous partial differential equations, but emphasize that the ideas are quite general and can be applied to other problems like finite-difference equations, coupled-map lattices or integro–differential equations.

We assume the dynamics of the system under consideration to be governed by a PDE of the form

$$[\stackrel{}{},_t,\stackrel{}{u}(\stackrel{}{x},t)]=0,$$
(1)   
where $`\stackrel{}{u}`$ is the field variable with $`N`$ components, $``$ is an $`N`$–dimensional function of $`\stackrel{}{u}(\stackrel{}{x},t)`$ and its spatial and temporal derivatives, with $`t`$ and $`\stackrel{}{x}`$ the temporal and spatial variables, respectively.

To ease the analysis and the representation, we only consider systems with a single component $`u`$ and at most two spatial dimensions. We note that the following considerations are valid in principle also for multi–component systems and higher dimensions in space.

In the following, we discuss the procedure to estimate a PDE of the form (1) from experimental data. Therefore, we distinguish between the solution $`u(\stackrel{}{x},t)`$ and the data $`v(\stackrel{}{x},t)`$. For the sake of simplicity, we denote both the continuous space–time variables of the model field and the discrete space–time variables of the data by the same symbols $`\stackrel{}{x}`$ and $`t`$. Considering the data $`v`$ as random variables, we act on a probability space and denote all entities estimated from the data by a hat $`\widehat{}`$ .

Since one can get only an estimate for the true function $``$ from the data $`v`$, all one can achieve is to estimate the analogue of Eq. (1)

$$\widehat{}[\widehat{},\widehat{_t},v(\stackrel{}{x},t)]=0,$$
(2)   
where $`\widehat{}`$ is the estimate of $``$ and the derivatives have been substituted by estimates computed from the data $`v`$.

To obtain $`\widehat{}`$, we formulate from (1) the following optimization problem :

$$\underset{}{inf}=e^2.$$
(3)   
The optimization lies in varying $``$ until $`e^2`$ converges to the infimum.

The function $``$ is defined as a function of operators, e.g. $`^2,_t,\mathrm{id},\mathrm{}`$. We denote the set of constituting operators by $`\{𝒪_i\}_{i=1,\mathrm{},K}`$. Note that in our definition of the differential operators $`𝒪_i`$ also any product terms like $`u^2_{xx},_t_x`$ are included. Since we consider a nonlinear problem, it is useful to split the function $``$ into sub–functions $`\mathrm{\Phi }_i`$, which have as arguments $`𝒪_iu`$, e.g. $`\mathrm{\Phi }(_{xx}u)=(_{xx}u)^\alpha `$ or $`\mathrm{\Phi }(\mathrm{id}u)=u^\beta `$. These $`\mathrm{\Phi }_i`$ are elements of some class of functions $`S`$, which has to be specified according to the problem. Thus, the final representation of Eq. (1) reads

$$=\underset{i=0}{\overset{K}{}}\mathrm{\Phi }_i(𝒪_iu)=0.$$
(4)   
We want to determine the constituents $`𝒪_iu`$ and the functions $`\mathrm{\Phi }_i`$ of $``$ from a data set $`v`$. Therefore, we estimate at first $`𝒪_iu`$ from the data by finite differencing or alternative schemes. Especially in the presence of noise, Fourier methods or kernel estimation could be helpful. The result consists of $`K`$ random variables $`\widehat{𝒪}_i\widehat{𝒪_iu}`$. If the underlying PDE is linear, i.e. the function $``$ a multilinear one, we could solve the problem with linear regression methods. Secondly, to obtain a solution for the nonlinear problem, we solve

$$\underset{\mathrm{\Phi }S}{inf}\underset{i=1}{\overset{K}{}}\mathrm{\Phi }_i(\widehat{𝒪}_i)=e^2$$
(5)   
in varying the functions $`\mathrm{\Phi }_i`$. The result are the estimates $`\widehat{\mathrm{\Phi }}_i`$. Up to now, we anticipated to know the number and type of operators $`𝒪_i`$. For unknown systems, which are the ones we treat, it would be necessary to extend the number of operators $`K`$ to infinity. In practice one has to select a finite number. Redundant terms then deliver $`\widehat{\mathrm{\Phi }}_i0`$ as result. It is important to find a reasonable initial guess for the operators $`𝒪_i`$ of Eq. (4). If there exists already some description for the system in a special state, one starts with the operators which appear in the known equations and tries to determine additional terms which may appear when leaving that state. This situation appears e.g. near some critical points, where one can start with already derived amplitude equations.

In the case of data analysis, the optimization problem (5) becomes a multiple nonlinear regression problem which can be solved using the alternating conditional expectation algorithm (ACE) . It has already been successfully applied in related fields of data analysis . In this case the norm of Eq. (5) is the $`L_2`$ norm. Using this algorithm, the functions $`\mathrm{\Phi }_i`$ of Eq. (4) are seen as so-called optimal transformations , which are defined to solve regression problems of the form

$$E[(\mathrm{\Phi }_0(\widehat{𝒪}_0)\underset{i=1}{\overset{K}{}}\mathrm{\Phi }_i(\widehat{𝒪}_i))^2]\stackrel{!}{=}\mathrm{min}.$$
(6)   
The $`\mathrm{\Phi }_i`$ are varied in the space $`S`$ of the Borel–measurable functions with the additional requirement of zero expectation, $`E[\mathrm{\Phi }_i]=0(i=0,\mathrm{},K)`$, and $`E[\mathrm{\Phi }_0^2]=1`$. The convergence of the ACE algorithm in the case of discrete samples rather than random variables has also been proved in . The minimization of the expectation value (6) is equivalent to the calculation of the maximal correlation $`\mathrm{\Psi }`$ . Instead of $`e^2`$ we use the maximal correlation as a measure for the quality of the result: $`\mathrm{\Psi }[0,1]`$ equals 1 for perfect estimation, the smaller it is the worse is the estimation. The ACE–algorithm achieves the maximization by iteratively transforming each $`\widehat{𝒪}_i`$ by suitable, generally nonlinear, transformations such as to yield a linear relationship between the new random variables $`\mathrm{\Phi }_i(\widehat{𝒪}_i)`$.

As a first example, we analyze data $`v(\stackrel{}{x},t)`$ from the numerical integration of the Swift–Hohenberg equation . The model is of the form

$`_tu`$ $`=`$ $`\left[r(^2+k^2)^2\right]uu^3`$ (7)   
$`=`$ $`(rk^4)uu^32k^2(_{xx}+_{yy})u`$ (9)   
$`(_{xxxx}+_{yyyy}+2_{xxyy})u.`$   
The global dynamics of the model can be derived from a potential, such that the asymptotic time dependence is trivial . Therefore, we analyze a transient state to have a sufficient variation in the time derivative. For the identification procedure we use data produced by an explicit Euler integration scheme with a time step of $`10^4`$, a spatial discretization of $`\mathrm{\Delta }x=\mathrm{\Delta }y=0.25`$, and periodic boundary conditions. As initial conditions we choose uniformly distributed independent random numbers from the interval $`[10,10]`$. The parameters are $`r=0.1`$ and $`k=1`$.

The differential operators are estimated by symmetric differencing schemes, e.g. $`\widehat{_t}v(\stackrel{}{x},t)=(v(\stackrel{}{x},t+\mathrm{\Delta }t)v(\stackrel{}{x},t\mathrm{\Delta }t))/2\mathrm{\Delta }t`$. Thus, to estimate the time derivatives of first order in each spatial data point, we need three consecutive ”pictures” of data. The field size is $`100\times 100`$ points, i.e., the data set $`v(\stackrel{}{x},t)`$ contains $`310^4`$ values. The data for the central time point are shown in Fig. 1. In the following we can drop without ambiguity the hat $`\widehat{.}`$ .

To identify the unknown system, we use an ansatz with non–mixed terms (like $`_xv`$) up to fourth order in the spatial derivatives. To show how to handle mixed terms we additionally include the terms $`v_xv`$, $`v_yv`$, $`_xv_yv`$.

$`\mathrm{\Phi }_0(_tv)`$ $`=`$ $`\mathrm{\Phi }_1(v)+\mathrm{\Phi }_2(_xv)+\mathrm{\Phi }_3(_yv)`$ (15)   
$`+\mathrm{\Phi }_4(_{xx}v)+\mathrm{\Phi }_5(_{xy}v)+\mathrm{\Phi }_6(_{yy}v)`$   
$`+\mathrm{\Phi }_7(_{xxx}v)+\mathrm{}+\mathrm{\Phi }_{10}(_{yyy}v)`$   
$`+\mathrm{\Phi }_{11}(_{xxxx}v)+\mathrm{}+\mathrm{\Phi }_{15}(_{yyyy}v)`$   
$`+\mathrm{\Phi }_{16}(v_xv)+\mathrm{\Phi }_{17}(v_yv)`$   
$`+\mathrm{\Phi }_{18}(_xv_yv)`$   
with the nine redundant terms $`\mathrm{\Phi }_2(_xv)`$, $`\mathrm{\Phi }_3(_yv)`$, $`\mathrm{\Phi }_5(_{xy}v)`$, $`\mathrm{\Phi }_7(_{xxx}v)`$, $`\mathrm{\Phi }_8(_{xxy}v)`$, $`\mathrm{\Phi }_9(_{xyy}v)`$, $`\mathrm{\Phi }_{10}(_{yyy}v)`$, $`\mathrm{\Phi }_{12}(_{xxxy}v)`$, $`\mathrm{\Phi }_{14}(_{xyyy}v)`$, $`\mathrm{\Phi }_{16}(v_xv)`$, $`\mathrm{\Phi }_{17}(v_yv)`$, $`\mathrm{\Phi }_{18}(_xv_yv)`$. We choose this ansatz as a compromise between generality and computational effort. Comparing Eq. (9) with Eq. (15), one expects in particular the following for the solution of Eq. (6): Up to an arbitrary common factor, $`\mathrm{\Phi }_0`$ should be the identity, $`\mathrm{\Phi }_1`$ should be a polynomial of third order, and for $`i=4,6,11,13,15`$ the $`\mathrm{\Phi }_i`$ should be linear functions $`\alpha _i𝒪_i`$. All other estimates should vanish. Furthermore, we expect for the slopes of the linear functions, after division by the slope of the l.h.s to remove the arbitrary common factor, $`\alpha _4=\alpha _6=\alpha _{13}=2`$, $`\alpha _{11}=\alpha _{15}=1`$.

Performing the ACE–algorithm, we find a maximal correlation of $`0.9993`$ and optimal transformations as shown in Fig. 2. All functions approximate the expected shape well, and the terms which were expected to vanish are indeed very small compared to the other ones. Note that this does not mean that the values for the redundant terms are vanishing themselves but that these are independent from all the other terms involved. The comparison of the slope of the linear functions yields the possibility to estimate parameters; we obtain $`\alpha _4=\alpha _6=1.9`$, $`\alpha _{13}=2.0`$, $`\alpha _{11}=\alpha _{15}=1.0`$. The slopes of the terms which are expected to vanish are all smaller than $`0.03`$ in absolute value. Since it is difficult to estimate the errors of these values, it is recommended to check the range of validity of the reconstructed model by integration and comparison with the dynamics. In Fig. 3 the nonlinearity $`0.9uu^3`$ of Eq. (9) is compared explicitly with the estimate, taking now an ansatz with only the non–vanishing terms of the above result. Inside a range of 98% of the data values this term can be estimated with high accuracy.

¿From a practical point of view, it is essential to discuss the stability of the identification method against noise. The effect of noise is to increase the errors of the estimates for the partial differential operators applied to the data, and therefore also the value of the error estimate $`e^2`$. The identification detects the remaining correlations of the vector field and its derivatives via the function $``$. In general, if the vector field and its derivatives are still strongly correlated for a reasonable low noise level of a noisy system, the minimum according to Eq. (5) should still be detectable.

To examine the stability against noise, we disturb each data point $`v(\stackrel{}{x},t)`$ by additive Gaussian white noise with a standard deviation of 0.5% of the data standard deviation. The estimates of the partial derivatives are disturbed severely due to error propagation. In spite of this, using again Eq. (15), the overall shape of all functions can still be recovered satisfactory, while being distorted (Fig. 3). The maximal correlation has decreased to 0.926.

Two other systems which were analyzed with similar success as the example above, are: 1. The Kuramoto–Sivashinsky equation (in the form as given in Ref. ) in the fully chaotic regime. 2. A reaction–diffusion system with two components , where the nonlinearity in the inhibitor dynamics is a non–analytic function. We reconstructed that function with high accuracy for different states like a single rotating spiral and spiral–defect chaos.

Summarizing our experiences, we find the following requirements on the data: 1) In order to identify a system with N components it is sufficient to measure N independent variables. 2) The sampling of the vector fields in space as well as in time has to be appropriate to estimate the partial derivatives with respect to time and space properly. 3) In the preceding example we analyzed data at one time point but for an extended region in space. In other situations it may be useful to perform the analysis at a single spatial point but for a longer interval in time, or even a combination of both procedures. However, as the example shows, even a small amount of data can be sufficient. 4) Since the reconstructed functions can only be determined at points which are attained by the data and the estimates for the differential operators, the data have to show a sufficient variability in space and time. 5) The inference from the data to a PDE is often not unique. But within the errors of the algorithm, our method reasonably reconstructs the dynamics on the trajectory given by the data. 6) We believe that due to the statistical nature of the method a higher noise level can typically be compensated by a larger data set.

We would like to stress the fact that the procedure presented above is essentially independent of the part of phase space the system moves on and the boundary and initial conditions. Thus we do not require the dynamics to be close to or away from an attractor. This kind of analysis is thus applicable in the case of high-dimensional chaotic motion, as exhibited by nonlinear spatially extended systems, as well as transient motion. We did not discuss non-homogeneous or non-autonomous evolution laws, but in principle the arguments apply equally well, with some minor modifications, for those systems.

In conclusion, we presented a method for the identification of spatiotemporal systems by numerical reconstruction of the PDE which describes the system. The identification procedure consists in solving an optimization problem which results in a nonlinear regression problem. Using the ACE–algorithm we showed that this task can indeed be solved in the case of numerical examples. We consider the method as a useful tool for the analysis of spatiotemporal systems and expect it to find a large variety of applications. The limits of this method and the requirements on the data have been discussed. Future work will concentrate on an extension of the ACE–algorithm to the solution of multi–component problems, and on the application to real data. Extensions to more general systems are envisaged.

We acknowledge helpful discussions with H. Kantz, J. Kurths, J. Parisi, A. Pikovsky, M. Zaks, and financial support by the Max–Planck–Gesellschaft.