Stochastic modelling of dynamical systems with several attractors

TL;DR

This paper explores stochastic models for systems with multiple attractors, comparing mean-field and Langevin approaches through numerical simulations to better understand their long-term behavior.

Contribution

It introduces a mean-field approach for finite systems with multiple attractors and compares it with traditional Langevin models using numerical simulations.

Findings

Mean-field hypothesis captures bifurcation behavior in finite systems.

Numerical simulations show agreement between mean-field and Langevin approaches.

Long-term dynamics are better understood with the proposed stochastic modelling.

Abstract

The usual Langevin approach to describe systems driven by noise fails to describe the long time behavior of systems with multiple attractors. The solution of the associated linear Fokker-Planck equation is always unique, even though it might show one or more maxima. In this context, it is customary to call transitions to changes in the shape of the equilibrium distribution function and relate the maxima to the attractors. Some years ago, a theory was developed for a system with interacting elements or subunits that, starting from the Langevin description of all the variables, leads to bifurcating \textquotedblleft one-particle\textquotedblright\, distribution functions when the number of elements tends to infinity. In this paper, a mean-field hypothesis has been used to deal with systems with a finite number of elements. We carry out numerical simulations yielding bifurcation solutions for the probability density of a collective variable. We also compare the results of the mean-field hypothesis with those obtained with the Langevin approach for finite systems.

Equations

$${\dot{x}}_i(t)=x_i(t)-x_i^3(t)+\frac{\theta }{N}\sum _{j=1}^N(x_j(t)-x_i(t)),$$

$${\dot{x}}_i(t)=(1-\theta )x_i(t)-x_i^3(t)+\theta s(t),(i=1,2,\dots ,N).$$

$$s(t)=\frac{1}{N}\sum _{j=1}^N{x}_j(t).$$

$${\dot{X}}_i(t)=(1-\theta )X_i(t)-X_i^3(t)+\theta S(t)+{\xi }_i(t),$$

$$S(t)=\frac{1}{N}\sum _{j=1}^N{X}_j(t),$$

$$⟨{\xi }_i(t)⟩=0;⟨{\xi }_i(t){\xi }_j(s)⟩=2D{\delta }_{ij}\delta (t-s).$$

$$\frac{\partial f_N}{\partial t}$$

$$U=\sum _{i=1}^{N}V(x_i)+\frac{\theta }{4N}\sum _{i=1}^N\sum _{j=1}^N(x_i-x_j{)}^2,$$

$$V(x)=-\frac{x^2}{2}+\frac{x^4}{4}.$$

$$\frac{\partial f_1(x,t)}{\partial t}$$

$$⟨x(t)⟩={\int }_{-\infty }^{\infty }x{f}_1(x,t)dx.$$

$${\dot{X}}_i(t)=(1-\theta )X_i(t)-X_i^3(t)+\theta ⟨S(t)⟩+{\xi }_i(t),$$

$$P(s,t)=\int dx_{1}\cdots dx_{N}\;\delta\Big{(}s-\frac{1}{N}\sum x_{i}\Big{)}f_{N}(x_{1},\cdots,x_{N},t)$$

Full text

I Introduction

The addition of noise to nonlinear dynamical systems has unveiled a number of new phenomena with both fundamental and practical interest in many scientific areas CITAS1 . For deterministic systems with a single attractor, one can turn noise-free variables into fluctuating ones by introducing a “fluctuacting force”  to the equations of the model. This approach was pioneered by Langevin many years ago in his study of Brownian motion and, since then, it has been applied successfully to a great variety of problems in physics, chemistry and biology ARNOLD1981 .

Nevertheless, the aforementioned program fails when applied to deterministic systems having a number of coexisting attractors, as it might be the case when the dynamics of the system is nonlinear. If we introduce the thermal effects (noise effects) by means of the usual Langevin construction, the corresponding Fokker-Planck equation for the probability density turns out to be a linear equation. Then, the possibility of having several distribution functions does not exist GANG1990 . It is possible that, as the value of the noise term is varied, the shape of the distribution function changes, so that for some range of noise values it can be multimodal, while it becomes monomodal for other noise values. When the distribution is multimodal one can think of associating each of its maxima to the corresponding deterministic attractors. But, we should keep in mind that, for a finite system, the multimodality of the distribution function indicates that the noise is able to connect the maxima, and consequently the attractors. This is in sharp contrast with the fact that the deterministic attractors are not connected.

The behavior of the average value provided by the Langevin description can not even qualitatively coincide with that of the deterministic variables for all values of the parameters. For example, if the system symmetry leads to a multimodal stationary distribution with zero average value, the random trajectories will not remain indefinitely around one of the attractors, even for very small noise strengths. Instead, they will explore the entire space, yielding a zero average value. Furthermore, even for very small noise intensities, the fluctuations about the average value will be very large, in consonance with the fact that the noise forces the random trajectories to explore all the attractors. On the other hand, as the noise strength is increased, the random trajectories will explore more easily the entire space, thus connecting the attractors much more frequently and leading, for a sufficiently large noise value, to a monomodal stationary distribution. Consequently the fluctuations will tend to decrease in intensity as the noise strength is increased. This is bizarre as one is normally used to the idea that thermal fluctuations tend to increase as the temperature increases.

Consider a finite dynamical system characterized by, say, the non random variables $x_{i},\,i=1,\ldots,N$ satisfying a set of nonlinear equations $\dot{x}_{i}=F_{i}(x_{1},\ldots,x_{j};\lambda)$ where $\lambda$ refers to some parameters. Imagine that the observation of the system leads to the conclusion that for some parameter values the system —as characterized by a global variable $s(t)$— has a single attractor, while for other set of parameter values, several attractors might be possible. Each of those attractors has its own basin of attraction. Then, the behavior of the system with time and its asymptotic state depends on its initial preparation.

Imagine now that temperature plays a role, so that one observes that the transition between a single attractor to several attractors might depend not only on the system parameters, but also on the temperature. The presence of temperature will imply that the variables $x_{i}$ will no longer be deterministic, but they will have fluctuations. Let’s denote by $X_{i}(t)$ the corresponding random variables. The question that arises is: what is the correct set of stochastic evolution equations compatible with the observations?

An example that comes to mind is the following: consider a permanent magnet which has been magnetized with its magnetic moment along a certain axis. If one does not apply any other external field, the magnetization will remain along the initial axis for as long as we observe it. Had we magnetized the sample along a different direction, it would stay like that for later times. Spontaneous transition of the magnetization among those directions are not observed as the temperature remains fixed. Now, suppose that we heat up the sample and the magnetization dissapears after a certain temperaure has been reached. Furthermore, assume that there are small fluctuations of the system magnetizations around each attractor value. Those fluctuations are not large enough to connect the two attractors during any reasonable amount of time. The nonstochastic dynamics can not describe these fluctuations. If the stochastic dynamics is modelled by the addition of a unbounded noise term representing the thermal effects, then we have fluctuations. But, if the Fokker-Planck description of the probability distribution is a linear one, the stationary distribution is unique. Thus, there must be fluctuations as large as one needs to connect the two attractors and the average magnetization will, sooner or later, be zero. If this is consistent with the experiments, fine. But if the system at the given temperature does not lose its magnetization in a spontaneous way, the Langevin (or linear Fokker-Planck) description is not correct.

II The model

Let us consider a model system consisting of $N$ “coordinates”  $x_{i}$, whose dynamics are given by the set of coupled differential equations (in dimensionless form)

[TABLE]

where $\theta$ is the strength of the coupling. We can rewrite this as

[TABLE]

where $s(t)$ is a collective variable defined by

[TABLE]

The numerical solution of Eqs.(2) leads to an $s(t)$ that always ends in one of the attractors of the individual dynamics ($1$ or $-1$). However, the temporal evolution of the individual variables is very different to that of the $s(t)$ and depends on initial conditions. The collective variable is representative of the global behavior of the system even for the deterministic dynamics.

A stochastic version of the model was introduced some years ago by Kometani and Shimizu KOMETANI1975 in a biophysical context and it was later analyzed by Desai and ZwanzigDESAI1978 , DawsonDAWSON1983 and ShiinoSHIINO1987 from a more general statistical mechanical perspective. The Langevin description amounts to add a stochastic term to each one of the deterministic equation, thus giving

[TABLE]

where

[TABLE]

is the collective stochastic variable and $\xi_{i}(t)$ are Gaussian white noises with

[TABLE]

An alternative formulation can be casted in terms of the linear Fokker-Planck equation for the joint probability distribution $f_{N}(x_{1},x_{2},\ldots,x_{N},t)$:

[TABLE]

where $U$ is the potential energy relief,

[TABLE]

with the single particle potential

[TABLE]

A problem with the above stochastic description is that, for any finite system, the joint probability distribution satisfies a linear Fokker-Planck equation. Consequently, the multiplicity of solutions in the deterministic limit does not manifest in the stochastic description as the joint probability distribution is unique. For some values of the parameters, the unique stationary distribution for the collective variable might show two peaks, reminiscent of the two deterministic attractors. Regardless of the initial conditions, the stochastic trajectories will explore the two attractors, at variance with the deterministic limit where the two attractors are independent.

As pointed out by Desai and Zwanzig, one can propose a cumulant moment expansion treatment of the linear Fokker-Planck equation. Truncation of the cumulant expansion at some reasonable level (say, the Gaussian level) leads to a nonlinear evolution equation for the averages. For some values of the parameters, this truncation provides an adequate qualitative description ($\theta>1)$ for the model at hand, but for other values of the parameter, any truncation leads to wrong qualitative results Drozmor1996 .

In the limit $N\to\infty$, Desai and Zwanzig DESAI1978 were able to write a Nonlinear Fokker-Planck Equation (NLFPE) (see also FRANK ; KURSTEN2016 ) describing the behavior of the “one-particle”  probability distribution,

[TABLE]

where

[TABLE]

In this approach, the infinite system undergoes an order-disorder phase transition for some values of parameters $\theta$ and $D$. This phase transition is signaled by a stochastic bifurcation leading to multiple equilibrium distributions $f_{eq}(x)$. In Fig.1 the bifurcation line separating the zones of one and three stationary solutions in the parameter space $\theta$ and $|z|=|\theta-1|/\sqrt{2D}$ has been depicted DESAI1978 .

Let us introduce an alternative approach to the description of a system with a finite number of elements. We will substitute the stochastic variable $S(t)$ in the finite set of equations (4) by its mean field version in the Weiss spirit, thus writing

[TABLE]

where $\langle\ldots\rangle$ represents an average over the noise. The interactions between the $X_{i}(t)$ variables are then replaced by the average value of the collective one. Note that within this description, the joint probability distribution $f_{N}(x_{1},x_{2},\ldots,x_{N},t)$ factorizes as a product of single particle distributions, each of them satisfying a nonlinear Fokker-Planck, Eq. (10). Also note that $\langle S(t)\rangle=\langle x(t)\rangle$. The probability distribution for the collective variable is defined by

[TABLE]

A purpose of this work is to study the distribution function $P_{eq}(s)$ for the collective variable of a finite system in the range of parameters where there are multiple attractors and to compare the results with those obtained with the finite set of Langevin equations, Eq.  (6). We are also interested in analyzing the relaxation of the average collective behavior towards a stationary situation, starting from different initial conditions.

III Numerical results

We have not been able to find a suitable analytical approximation for the process $S(t)$. Thus, we will rely on numerical simulations to get the results presented here. In what follows, we have chosen $\theta=0.5$, thus allowing the noise strength $D$ to be the only bifurcation parameter of the system. Other values of the parameter $\theta$ lead to qualitatively analogous results.

We have solved numerically the systems of equations, Eqs. (4) and (12) for a small number of particles $(N=11)$ using different realizations of the noise. After a convenient relaxation time so that we are sure a stationary probability distribution has been reached, we record during a long time interval and for all the random trajectories, the times that the collective variable is within the different bins of the suitably discretized collective variable. In this form we have been able to construct histograms that approximate the equilibrium probability distributions $P_{eq}(s)$ for the Langevin and the mean-field descriptions.

In Fig.2 we have represented the equilibrium distribution $P_{eq}(s)$ for the collective variable in the mean-field description for some values of the noise strength $D$. In all cases the initial conditions of the individual elements were chosen so that $S(0)<0$. As we can observe, there is a range in $D$ where $P_{eq}(s)$ is a zero-centered, monomodal distribution while, as the value of $D$ decreases, the distribution peaks around a finite (negative) value of $s$. For initial conditions with $S(0)>0$, and the same small values of the noise strength, a symmetric $P_{eq}(s)$ (not depicted) is reached that peaks around a positive value of $s$.

The above behavior is in contrast with that obtained by using the Langevin description. In this case, as shown in Fig. 3, the equilibrium distribution $P_{eq}(s)$ is always unique regardless of whether the initial preparation yields $S(0)<0$ or $S(0)>0$. The equilibrium distribution is always zero centered, although its shape changes from a single peak distribution at high noise intensities to a bimodal one as the noise strength is decreased. This is to be expected as the strong noise blends the two symmetric deterministic attractors into a single equilibrium point with the highest probability. For low noise values, the random trajectories are essentially fluctuating around the deterministic attractors, except for some sporadic transitions between them. As the system size increases, the two maxima of the bimodal distribution become narrower while the local minima around $s=0$ becomes much smaller. It is not surprising, then, that in the limit $N\to\infty$, the two peaks are completely separated and one finds two equilibrium monomodal distributions, either one of them being reached depending on the initial conditions.

A further comparison has been made between these two descriptions by calculating the dependence of the average value $\langle S\rangle_{eq}$ and dispersion $\sigma_{s}=\sqrt{\langle S^{2}\rangle_{eq}-\langle S\rangle_{eq}^{2}}$ on $D$ for both cases. We have numerically evaluated the integrals leading to $\langle S\rangle_{eq}$ and $\langle S^{2}\rangle_{eq}$ using the histograms. As we can observe in Fig. 4, the average value is always zero for the Langevin description whether or not the distribution is monomodal o bimodal, due tho the symmetric character of the distribution. For the mean-field description, however, the average value bifurcates from zero to finite (positive or negative) values of $s$ as $D$ is decreased.

The behavior of the dispersion as $D$ changes is also very different in both descriptions. As we can see in Fig.5, as the noise strength is increased, the value of $\sigma_{s}$ increases in the case of mean-field dynamics and decreases in the Langevin description. The fact that the dispersion decreases as the noise increases in the Langevin description is tied to the fact that the corresponding equilibrium distribution $P_{eq}(s)$ changes from having a single peak for large noise values to be bimodal at low noise intensities.

Another aspect of our analysis is that of relaxation towards equilibrium. It is useful to compare the relaxation time of the average behavior with that of the deterministic relaxation and also to compare the relaxation times in both stochastic descriptions for the same values of the noise strength. This last comparison can be made by observing Figs. 6 and 7. For small noise values, the differences in the relaxation times are noteworthy.

In Fig. 8 we depict an example of the relaxation of the average collective variable as given by the two stochastic descriptions analyzed above and by the deterministic dynamics. As we can see, the Langevin dynamics shows an average behavior towards a zero value, even though the deterministic dynamics shows a relaxation towards $s=-1$, the attractor reached by the initial condition $S(0)=-0.3$. The mean-field dynamics, on the other hand, indicates a relaxation towards a negative value, which differs from the deterministic attractor due to the renormalization of the deterministic nonlinear dynamics by the coupling of noise and nonlinearity, a well known effect in stochastic nonlinear systems NORDHOLM1974 .

IV Conclusions

We have discussed in this paper a possible way to analyze noise effects in the description of finite size nonlinear systems which, in the deterministic limit, would present two disconnected stationary attractors. The variable of interest is a collective variable characterizing the whole system. Due to the nonlinear nature of the problem, an analytical solution is out of the question and we have relied on numerical simulations.

We have compared the characteristics of the equilibrium distribution function of the collective variable for two different approaches to describe the fluctuating dynamics of a finite system with interacting subunits. A possible description is to add a noise term to each of the individual deterministic degrees of freedom as it is typically done in a Langevin description. As an alternative, we propose a stochastic description leading to nonlinear Fokker-Planck equations. In this approach, we make use of an extension of the Weiss procedure to explain the paramagnetic-ferromagnetic transition, and incorporate the collective variable average value into the stochastic dynamics of each variable.

Our numerical results show that the equilibrium probability distribution for the collective variable $P_{eq}(s)$ in the mean-field approach presents a bifurcating character signaling the occurrence of a order-disorder transition. For high values of the noise strength, $P_{eq}(s)$ is an almost Gaussian distribution centered around zero, regardless of the initial condition. As the noise strength is decreased, two stable stationary distributions centered around positive or negative values are possible. Which one of them is reached depends on the initial preparation.

The symmetrical mean values obtained from the bifurcating solutions are not identical to the stationary values in the deterministic description. This is not surprising, due to the renormalization effects associated to the noise, nonlinearity and cooperativity.