Approximations in Mean Square Analysis of Stochastically Forced Equilibria for Nonlinear Dynamical Systems

1. Introduction

Currently, nonlinear stochastic phenomena such as noise-induced transitions [1,2], stochastic excitement [3,4,5], noise-induced crisis [6,7], stochastic bifurcations [8,9], noise-induced chaos [10,11], stochastic and coherence resonances [12,13,14] are being actively studied in various fields of natural sciences. One of the key mechanisms of such effects is associated with the transition of random trajectories through separatrices detaching basins of coexisting attractors. For the parametric analysis of the conditions causing such transitions and the estimation of the corresponding critical values of the intensity of random disturbances that generate them, a new fairly general geometric approach has been used, based on the method of confidence domains [15,16,17].

The main idea of this approach is as follows. As the noise intensity increases, the size of the confidence domain increases. The critical value of the noise intensity is found from the condition of the intersection of the confidence domain and the separatrix. The key parameter that determines the configuration of the confidence domain is the stochastic sensitivity of the attractor [18,19]. Initially, the stochastic sensitivity technique was introduced in connection with the approximation [20] of a quasipotential [21] in the vicinity of an attractor. Currently, this technique has been developed for regular and chaotic attractors of both continuous and discrete systems (see, e.g. [22,23,24]). The stochastic sensitivity technique and the associated confidence domains method are actively used in the analysis of nonlinear stochastic phenomena [25,26,27,28] and control problems [29,30].

Using the stochastic sensitivity technique, an approximation of the mean square deviations of random states of a stochastic system from the deterministic attractor is constructed. Although formally this technique is applicable to both the case of additive and parametric noise, however, in the case of parametric noise, the error in the corresponding approximations may be such that the prediction made on the basis of this approximation may turn out to be incorrect.

This paper is devoted to the problem of approximation of probabilistic distributions of random states around stable equilibria of stochastic differential Ito’s equations with general multiplicative noise. Using first approximation systems, we construct an approximation of mean square deviations that explicitly takes into account the presence of parametric noise. This more accurate approximation is compared with an approximation based on the stochastic sensitivity technique. General mathematical results are illustrated with examples.

2. Mean Square Analysis of Stochastic Equilibria

Consider a nonlinear autonomous system of ordinary differential equations

$$dx=f\left(x\right)dt,$$

(1)

where x is an n-dimensional vector and $f\left(x\right)$ is a sufficiently smooth n-vector function. It is assumed that the system (1) has an exponentially stable equilibrium $\overline{x}$.

Definition 1. The equilibrium $\overline{x}$ is called exponentially stable in system (1) if for some neighborhood $\mathrm{U}$ of $\overline{x}$ there exist constants $K>0,$$l>0$ such that for all $t⩾0$ it holds that

$$\parallel x\left(t\right)-\overline{x}\parallel ⩽K{e}^{-lt}\parallel x_0-\overline{x}\parallel ,$$

where $x\left(t\right)$ is a solution of the system (1) with the initial condition $x\left(0\right)=x_0\in \mathrm{U}.$ Here, $\parallel ·\parallel$ is the Euclidean norm.

Along with the deterministic system (1), let us consider the stochastic Ito system

$$dx=f\left(x\right)dt+\sum _{i=1}^k{\sigma }_i\left(x\right)d{w}_i\left(t\right),$$

(2)

where ${\sigma }_i\left(x\right)$ are sufficiently smooth n-vector functions, $w_i\left(t\right)$ are scalar standard independent Wiener processes. The functions ${\sigma }_i\left(x\right)$ model the dependence of multiplicative disturbances on the system state.

Solutions of the stochastic system (2), leaving the deterministic equilibrium $\overline{x}$ under the influence of random disturbances, form some probability distribution. It is assumed that the probabilistic distribution of the states of system (2) stabilizes. The corresponding stable stationary distribution density satisfies the Fokker-Planck-Kolmogorov equation [21,31]. It is known that in the general case it is very difficult to directly use this equation to describe probability distributions, even for two-dimensional systems. Here, the apparatus of first approximation systems is useful.

2.1. First Approximation System and Its Mean Square Analysis

Let us consider the deviation $z=x-\overline{x}$ of the random state x of the system (2) from the exponentially stable equilibrium $\overline{x}$ of the system (1). Dynamics of the variable z is governed by the following first approximation linear system:

$$dz=Fzdt+\sum _{i=1}^k(S_{0,i}+S_{1,i}z)d{w}_i,$$

(3)

where

$${\displaystyle F=\frac{\partial f}{\partial x}\left(\overline{x}\right),S_{0,i}={\sigma }_i\left(\overline{x}\right),S_{1,i}=\frac{\partial {\sigma }_i}{\partial x}\left(\overline{x}\right).}$$

In our mean square analysis of the system (3) solutions, we will use first (m) and second (M) moments: $m=\mathrm{E}\left(z\right),$$M=\mathrm{E}\left(z{z}^{\top }\right)$. Dynamics of these deterministic characteristics is described by the following equations:

$$\dot{m}=Fm,$$

(4)

$$\dot{M}=FM+M{F}^{\top }+\sum _{i=1}^k\left(S_{0,i}{S}_{0,i}^{\top }+S_{0,i}{m}^{\top }{S}_{1,i}^{\top }+S_{1,i}m{S}_{0,i}^{\top }+S_{1,i}M{S}_{1,i}^{\top }\right).$$

(5)

To find an approximation of the dispersion of stationary distributed random states of the nonlinear stochastic system (2) around the deterministic equilibrium $\overline{x}$, we will use stationary solutions of the system (4), (5).

Due to exponential stability of the equilibrium $\overline{x}$, it holds that $\mathrm{Re}{\lambda }_i\left(F\right)<0$, where ${\lambda }_i\left(F\right),i=1,\dots ,n,$ are eigenvalues of the matrix F. In these circumstances, the system (4) has a unique stationary stable solution $\overline{m}=0$. Substituting $\overline{m}=0$ into (5), we get

$$\dot{M}=FM+M{F}^{\top }+\sum _{i=1}^k\left(S_{0,i}{S}_{0,i}^{\top }+S_{1,i}M{S}_{1,i}^{\top }\right).$$

(6)

So, the matrix $\overline{M}$ of the stationary solution of the equation (6) satisfies the following algebraic equation:

$$F\overline{M}+\overline{M}{F}^{\top }+\sum _{i=1}^k\left(S_{0,i}{S}_{0,i}^{\top }+S_{1,i}\overline{M}{S}_{1,i}^{\top }\right)=0.$$

(7)

Let us consider a deviation $\Delta \left(t\right)=M\left(t\right)-\overline{M}$ where $M\left(t\right)$ is a solution of the equation (6). For the function $\Delta \left(t\right)$, one get the homogeneous equation

$$\dot{\Delta }=F\Delta +\Delta F^{\top }+\sum _{i=1}^k{S}_{1,i}\Delta S_{1,i}^{\top }.$$

(8)

The matrix $\Delta \left(t\right)$ is the matrix of second moments $\Delta \left(t\right)=\mathrm{E}y\left(t\right)y{\left(t\right)}^{\top }$ for solutions $y\left(t\right)$ of linear homogeneous stochastic equation

$$dy=Fydt+\sum _{i=1}^k{S}_{1,i}yd{w}_i.$$

(9)

Thus, the question about stability of the stationary solution $\overline{M}$ of the equation (6) is reduced to the equivalent question about the exponential mean square stability of the trivial solution $y=0$ of the stochastic system (9).

Definition 2. Solution $y\equiv 0$ of the stochastic system (9) is called exponentially stable in mean square, if there exist constants $K>0,$$l>0$ such that for all $t⩾0$ it holds that

$${\mathrm{E}\parallel y\left(t\right)\parallel }^2⩽K{e}^{-lt}\mathrm{E}{\parallel y_0\parallel }^2,$$

where $y\left(t\right)$ is a solution of the system (9) with the initial condition $y\left(0\right)=y_0$.

Let us consider the matrix $S_0={\sum }_{i=1}^k{S}_{0,i}{S}_{0,i}^{\top }$ and operators

$$\mathcal{A}\left[M\right]=FM+M{F}^{\top },\mathcal{S}\left[M\right]=\sum _{i=1}^k{S}_{1,i}M{S}_{1,i}^{\top },\mathcal{P}=-{\mathcal{A}}^{-1}\mathcal{S}.$$

Rewrite equations (6), (7), and (8) as follows:

$$\dot{M}=\mathcal{A}\left[M\right]+\mathcal{S}\left(M\right)+S_0,$$

(10)

$$\mathcal{A}\left[\overline{M}\right]+\mathcal{S}\left[\overline{M}\right]+S_0=0,$$

(11)

$$\dot{\Delta }=\mathcal{A}[\Delta ]+\mathcal{S}[\Delta ].$$

(12)

Note that the existence of the operator ${\mathcal{A}}^{-1}$ follows from the condition $\mathrm{Re}{\lambda }_i\left(F\right)<0$.

Basic theoretical connections are presented in the following theorem.

Theorem 1. 

The following statements are equivalent:

(a) System (10) has a stationary exponentially stable solution $\overline{M}$ satisfying (11);

(b) The solution $\Delta \equiv 0$ of the system (12) is exponentially stable;

(c) The solution $y\equiv 0$ of the stochastic system (9) is exponentially stable in mean square;

(d) It holds that $\mathrm{Re}{\lambda }_i\left(F\right)<0(i=1,...,n)$ and $\rho \left(\mathcal{P}\right)<1$, where $\rho \left(\mathcal{P}\right)$ is the spectral radius of the operator $\mathcal{P}$.

The statements of this theorem were proven or follow from more general results presented in [32,33,34,35].

Remark. 

In the one-dimensional case ($n=1$), we have

$$S_0=\sum _{i=1}^k{\sigma }_i^2\left(\overline{x}\right),\mathcal{A}\left[M\right]=2f^{\prime }\left(\overline{x}\right)M,\mathcal{S}\left[M\right]=\sum _{i=1}^k{\left({\sigma }_i^{\prime }\left(\overline{x}\right)\right)}^{2}M,$$

and the condition $\rho \left(\mathcal{P}\right)<1$ has an explicit parametric representation:

$$\rho \left(\mathcal{P}\right)=-\frac{1}{2f^{\prime }\left(\overline{x}\right)}\sum _{i=1}^k{\left({\sigma }_i^{\prime }\left(\overline{x}\right)\right)}^2<1.$$

In this case, for the mean square variance of random states around the equilibrium $\overline{x}$ the following estimation can be written:

$$M=-\frac{{\sum }_{i=1}^k{\sigma }_i^2\left(\overline{x}\right)}{2f^{\prime }\left(\overline{x}\right)+{\sum }_{i=1}^k{\left({\sigma }_i^{\prime }\left(\overline{x}\right)\right)}^2}.$$

(13)

2.2. Asymptotics for the Case of Weak Noise, Stochastic Sensitivity of the Equilibrium

Consider the stochastic system

$$dx=f\left(x\right)dt+\epsilon \sum _{i=1}^k{\sigma }_i\left(x\right)d{w}_i\left(t\right),$$

(14)

where $\epsilon$ is a scalar small parameter of the disturbance intensity. For this system, the equation (11) for the covariance matrix M of the equilibrium $\overline{x}$ has the form

$$\mathcal{A}\left[M\right]+{\epsilon }^2\mathcal{S}\left[M\right]+{\epsilon }^2{S}_0=0.$$

Let us study the dependence of the solution $M\left(\epsilon \right)$ of this equation on the parameter $\epsilon$. Let $W\left(\epsilon \right)$ be the solution to the equation

$$\mathcal{A}\left[W\right]+{\epsilon }^2\mathcal{S}\left[W\right]+S_0=0.$$

(15)

Then $M\left(\epsilon \right)={\epsilon }^{2}W\left(\epsilon \right).$ For $W\left(\epsilon \right)$, one can write the following decomposition:

$$W\left(\epsilon \right)=-{(\mathcal{A}+{\epsilon }^2\mathcal{S})}^{-1}\left[S_0\right]=-{\left(I+{\epsilon }^2{\mathcal{A}}^{-1}\mathcal{S}\right)}^{-1}{\mathcal{A}}^{-1}\left[S_0\right]=$$

$$={\left(I-{\epsilon }^2\mathcal{P}\right)}^{-1}\left[W\left(0\right)\right].$$

For small $\epsilon$, it holds that

$$W\left(\epsilon \right)=\sum _{m=0}^{\infty }{\epsilon }^{2m}{\mathcal{P}}^m\left[W\left(0\right)\right]=W\left(0\right)+{\epsilon }^2\mathcal{P}\left[W\left(0\right)\right]+{\epsilon }^4{\mathcal{P}}^2\left[W\left(0\right)\right]+\dots$$

As a result, for the matrix function $M\left(\epsilon \right)$ we get the expansion in powers of the small parameter:

$$M\left(\epsilon \right)=\sum _{m=0}^{\infty }{\epsilon }^{2m+2}{\mathcal{P}}^m\left[W\left(0\right)\right]={\epsilon }^{2}W\left(0\right)+{\epsilon }^4\mathcal{P}\left[W\left(0\right)\right]+{\epsilon }^6{\mathcal{P}}^2\left[W\left(0\right)\right]\dots$$

(16)

In this series, the matrix $W\left(0\right)$ plays an important role in the asymptotic analysis of the spread of random states around the equilibrium. Because of $W\left(0\right)={lim}_{\epsilon \to 0}\frac{1}{{\epsilon }^2}M\left(\epsilon \right)$, this matrix characterizes the stochastic sensitivity of the equilibrium to the impact of weak noise. Thus, in the first approximation, we have

$$M\left(\epsilon \right)\approx M^{\left(1\right)}\left(\epsilon \right)={\epsilon }^{2}W,$$

(17)

where W is a solution of the following equation

$$FW+W{F}^{\top }+S_0=0.$$

(18)

If the noise in the system (14) does not depend on the state, then $\mathcal{P}=0$ and the first approximation coincides with the exact value: $M\left(\epsilon \right)={\epsilon }^{2}W.$ In general, using W as an approximation for $M\left(\epsilon \right)$ one get an underestimation of the covariance of random states. Indeed, since the operator $\mathcal{P}$ is positive, the inequality $M\left(\epsilon \right)⪰{\epsilon }^{2}W$ is valid.

In one-dimensional case, the stochastic sensitivity of the equilibrium $\overline{x}$ for the system (14) is given by the formula

$$W=-\frac{{\sum }_{i=1}^k{\sigma }_i^2\left(\overline{x}\right)}{2f^{\prime }\left(\overline{x}\right)}.$$

(19)

3. Examples

Let us consider how these theoretical results can be applied to the approximation of mean square deviation of random states from the equilibrium in some stochastic systems.

Example 1. 

Consider a simple one-dimensional stochastic system

$$dx=-axdt+\epsilon ({\sigma }_{1}d{w}_1\left(t\right)+{\sigma }_{2}xd{w}_2\left(t\right)),$$

where $a,{\sigma }_1,$ and ${\sigma }_2$ are non-negative parameters, $\epsilon$ is the intensity of random disturbances, $w_i\left(t\right)$ are uncorrelated scalar Wiener processes. The parameters ${\sigma }_1$ and ${\sigma }_2$ specify the weights of additive and multiplicative disturbances, accordingly.

For $a>0$, the corresponding deterministic system (with $\epsilon =0$) has an exponentially stable equilibrium $\overline{x}=0.$ Second moments $M\left(t\right)=\mathrm{E}{(x\left(t\right)-\overline{x})}^2$ of deviations of solutions $x\left(t\right)$ from the equilibrium satisfy the equation

$$\dot{M}=-2aM+{\epsilon }^2({\sigma }_1^2+{\sigma }_2^{2}M).$$

This equation has a stationary solution

$$M\left(\epsilon \right)=\frac{{\epsilon }^2{\sigma }_1^2}{2a-{\epsilon }^2{\sigma }_2^2}.$$

Following the decomposition (16), for weak noise $M\left(\epsilon \right)$ has the following asymptotics

$$M\left(\epsilon \right)={\epsilon }^{2}W+O\left({\epsilon }^4\right),$$

where W characterizes the stochastic sensitivity of the equilibrium $\overline{x}$. Here, Wy satisfies (see (18)) the equation:

$$-2aW+{\sigma }_1^2=0.$$

Using W one can write the first approximation for the function $M\left(\epsilon \right)$:

$$M^{\left(1\right)}\left(\epsilon \right)={\epsilon }^{2}W=\frac{{\epsilon }^2{\sigma }_1^2}{2a}.$$

Formally, the approximation $M^{\left(1\right)}\left(\epsilon \right)$ is defined for any $a>0$ while the approximated function $M\left(\epsilon \right)$ is defined only for $a>\frac{{\epsilon }^2{\sigma }_2^2}{2}.$ In absence of multiplicative noise $({\sigma }_2=0$), the values $M^{\left(1\right)}$ and M are identical. At ${\sigma }_2\ne 0$, they can essentially differ.

This difference is clearly seen in Figure 1 where plots of the functions M (solid line) and $M^{\left(1\right)}$ (dashed line) are shown versus parameter a. Note that the approximation $M^{\left(1\right)}$ is always less than M (this fact was shown above for the general case). Moreover, in the interval $0<a<\frac{{\epsilon }^2{\sigma }_2^2}{2},$ where the approximation gives finite values, the original function is not defined at all: the second moments $M\left(t\right)$ tend to infinity. In the interval $a>\frac{{\epsilon }^2{\sigma }_2^2}{2}$, the approximation error monotonically increases and tends to infinity as it approaches the bifurcation value $a_{*}=\frac{{\epsilon }^2{\sigma }_2^2}{2}$. For the relative error, an explicit representation can be written:

$$\left|\frac{M-M^{\left(1\right)}}{M}\right|=\frac{{\epsilon }^2{\sigma }_2^2}{2a}.$$

Let us continue the comparison of these two methods for estimating the dispersion of random states around the equilibrium using the two-dimensional system as an example.

Example 2. 

Consider the van der Pol model with hard excitation of self-oscillations:

$$\begin{array}{c}dx=ydt\hfill \\ dy=(-x+(a+b{x}^2-x^4)y)dt+{\sigma }_{1}d{w}_1+{\sigma }_{2}yd{w}_2.\hfill \end{array}$$

(20)

Here, ${\sigma }_1$ is the intensity of the additive noise, ${\sigma }_2$ is the intensity of the multiplicative noise, and $w_1\left(t\right),w_2\left(t\right)$ are scalar standard independent Wiener processes.

Let us fix $a=-0.1,b=2$. For this set of parameters, the deterministic system (20) with ${\sigma }_1={\sigma }_2=0$ is bistable and exhibits the coexisting attractors: the stable equilibrium $\overline{x}=0,$$\overline{y}=0$ and stable limit cycle. Basins of these attractors are separated by the orbit of the unstable limit cycle. In Figure 2a and Figure 3a, the equilibrium is shown by black filled circle, the stable cycle is plotted by blue curve, and the unstable cycle (the separatrix) is shown by red curve.

Let us consider the behavior of trajectories of the stochastic system (20) solutions starting at the equilibrium $(0,0)$. Under the influence of weak random disturbances, trajectories leave the stable equilibrium and form a stationary probability distribution concentrated in a small neighborhood of the origin. This type of dynamics corresponds to the unexcited mode of the oscillator (see Figure 2 for ${\sigma }_1=0.02,{\sigma }_2=0.2$).

As the noise intensity increases, random trajectories cross the separatrix (unstable limit cycle) and continue to oscillate near the stable cycle. This means a transition to the excitation mode (see Figure 2 for ${\sigma }_1=0.05,{\sigma }_2=0.4$).

For the analytical approximation of the dispersion of random states, we will use theory presented above.

For the system (20), parameters of the equation (7) are following:

$$F=\left[\begin{array}{cc}0& 1\\ -1& a\end{array}\right],M=\left[\begin{array}{cc}{m}_{11}& m_{12}\\ m_{21}& m_{22}\end{array}\right],$$

$$S_{0,1}={\sigma }_1\left[\begin{array}{c}0\\ 1\end{array}\right],S_{0,2}=\left[\begin{array}{c}0\\ 0\end{array}\right],S_{1,1}=\left[\begin{array}{cc}0& 0\\ 0& 0\end{array}\right],S_{1,2}={\sigma }_2\left[\begin{array}{cc}0& 0\\ 0& 1\end{array}\right].$$

Now, we can write the matrix equation (7) as the following system for the elements $m_{ij}$ of the symmetric matrix M:

$$-m_{11}+a{m}_{12}+m_{22}=0,m_{12}=m_{21}=0,-2m_{12}+2a{m}_{22}+{\sigma }_1^2+{\sigma }_2^2{m}_{22}=0.$$

From this system, we have solution

$$m_{11}=m_{22}=-\frac{{\sigma }_1^2}{2a+{\sigma }_2^2},m_{12}=m_{21}=0.$$

Thus, the matrix M that defines mean square deviation of random states from the equilibrium $(0,0)$ is

$$M=-\frac{{\sigma }_1^2}{2a+{\sigma }_2^2}\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].$$

Note that the asymptotic method of stochastic sensitivity gives for mean square deviation another approximation:

$$M^{\left(1\right)}=-\frac{{\sigma }_1^2}{2a}\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].$$

The difference in these approximations can lead to the qualitative differences in the prediction of the results of the noise influence. Let us consider how these two estimations work in the context of the confidence domains method. For diagonal $2\times 2$-matrices, the equation of the confidence ellipse is written as

$$\frac{x^2}{m_{11}}+\frac{y^2}{m_{22}}=-2ln(1-P).$$

Here, P is fiducial probability.

Confidence ellipses are effectively used in predicting noise-induced transitions through the separatrix. Figure 4 shows two confidence ellipses constructed using the matrices M (larger ellipse) and $M^{\left(1\right)}$ (smaller ellipse) for the stochastic system (20) with $a=-1,{\sigma }_1=0.05,{\sigma }_2=0.4.$ The larger ellipse captures the basin of attraction of the limit cycle, which allows us to make a prediction about the generation of large-amplitude oscillations (excitation mode). The smaller ellipse is entirely contained in the basin of attraction of the stable equilibrium and therefore predicts the unexcited mode of the oscillator. As we see, an error in estimating the second moments can lead to qualitative errors in solving important prediction problems.

4. Conclusions

The paper is devoted to the problem of approximation of probability distributions of random states near the equilibrium of the stochastic system with multiplicative noise. The system of nonlinear stochastic differential Ito’s equations is used as a basic mathematical model. For the linear first approximation system, we present conditions of existence of the mean square stable stationary distribution. For mean square deviation, an expansion in powers of the small parameter of noise intensity is derived. These approximations are compared with the widely used approximation based on the stochastic sensitivity technique. The general mathematical results are applied to the analysis of noise-induced generation of large-amplitude oscillations in the van der Pol oscillator with hard excitement.