CC BY
5Russian Journal of Nonlinear Dynamics, 2020, vol. 16, no. 3, pp. 437-451. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd200303
NONLINEAR PHYSICS AND MECHANICS
MSC 2010: 37G15, 37G35
Analysis of Special Cases in the Study of Bifurcations of Periodic Solutions of the Ikeda Equation
E. P. Kubyshkin, A. R. Moriakova
This paper deals with bifurcations from the equilibrium states of periodic solutions of the Ikeda equation, which is well known in nonlinear optics as an equation with a delayed argument, in two special cases that have not been considered previously. Written in a characteristic time scale, the equation contains a small parameter with a derivative, which makes it singular. Both cases share a single mechanism of the loss of stability of equilibrium states under changes of the parameters of the equation associated with the passage of a countable number of roots of the characteristic equation through the imaginary axis of the complex plane, which are in this case in certain resonant relations. It is shown that the behavior of solutions of the equation with initial conditions from fixed neighborhoods of the studied equilibrium states in the phase space of the equation is described by countable systems of nonlinear ordinary differential equations that have a minimal structure and are called the normal form of the equation in the vicinity of the studied equilibrium state. An algorithm for constructing such systems of equations is developed. These systems of equations allow us to single out one "fast" variable and a countable number of "slow" variables, which makes it possible to apply the averaging method to the systems of equations obtained. Equilibrium states of the averaged system of equations of "slow" variables in the original equation correspond to periodic solutions of the same nature of sustainability. In the special cases under consideration, the possibility of simultaneous bifurcation from equilibrium states of a large number of stable periodic solutions (multistability bifurcation) and evolution of these periodic solutions to chaotic attractors with changing bifurcation parameters is shown.
Received November 17, 2019 Accepted May 18, 2020
The work was carried out with the financial support of the P. G. Demidov Yaroslavl State University development Program for the period 2017-2021 (task 0P-2G-01-2019).
Evgenii P. Kubyshkin
kubysh.e@yandex.ru
Alena R. Moriakova
alyona_moryakova@mail. ru
P. G. Demidov Yaroslavl State University
ul. Sovietskaya 14, Yaroslavl, 150000 Russia
One of the special cases is associated with the formation of paired equilibrium states (a stable and an unstable one). An analysis of bifurcations in this case provides an explanation of the formation of the "boiling points of trajectories", when a periodic solution arises "out of nothing" at some point in the phase space under changes of the parameters of the equation and quickly becomes chaotic.
Keywords: Ikeda equation, periodic solutions, bifurcation of multistability, chaotic multi-stability
1. Introduction
In this paper we consider the delay differential equation
x = j sin(x(t — h) — c) — x, (1.1)
proposed in [1] to describe the dynamics of a passive optical resonator, where the variable x(t) determines the phase shift of the electric field in the nonlinear absorbing medium of the ring resonator, h is the time of propagation of light in the ring resonator, 0 ^ c < 2n is constant phase shift, and the coefficient j > 0 characterizes the intensity of laser radiation. Equation (1.1), called the Ikeda equation, has rich nonlinear dynamics, depending on the values of the parameters c, h and j, can have a large number of different stable and unstable self-oscillating solutions and is widely used in nonlinear optics to explain optical effects. In this regard, it is of great interest to study the mechanisms of the emergence of self-oscillatory solutions of Eq. (1.1) and their development depending on changes in the parameters of the equation. The authors of [2] investigated the equilibrium states of Eq. (1.1), their dynamics and stability depending on the parameters of the equation, as well as the nature of the loss of stability. Bifurcations of periodic solutions from equilibrium states for the most common case of stability loss are considered. The possibility of simultaneous bifurcation of a large number of periodic solutions (multistability bifurcation) and their development by changing the parameters of the equation to chaotic at-tractors (chaotic multistability) is shown. There is also a brief overview of the work of other authors. In [2], two special cases of loss of stability of equilibrium states, which are important for applications, have not been considered. The first is the case of the loss of stability of the zero equilibrium state when c = 0 at the point j = 1, the second is the case of formation of a pair of equilibrium equations — a simultaneous formation of stable and unstable equilibrium states "out of thin air". These two cases share a single mechanism of loss of stability, different from that considered in [2]. In both cases, the possibility of bifurcation of a large number of periodic solutions and their development into chaotic multistability will be shown. The second case allows one to explain the formation of peculiar points of "boiling of trajectories" when under a change of the parameters of the equation a periodic solution arises at some point of phase space "out of nothing" and quickly becomes a chaotic solution. The possibility of such behavior of solutions based on numerical analysis is pointed out in [3]. This work complements Ref. [2] and allows one to create a fairly complete picture of possible bifurcations of periodic solutions from equilibrium states of Eq. (1.1) and of the development of these solutions when changing the parameters of the equation to chaotic attractors, and to understand the structure of the phase space of the equation.
2. Mathematical statement of the problem
Let us proceed in Eq. (1.1) to the dimensionless time t' = t/h (which we will omit later). In optical systems, the delay time h is significantly longer than the characteristic time of oscillations of the electric field, so e1 = h-1 ^ 1. As a result, we have the equation
£\x(t) + x(t) = H sin(x(t — 1) — c). (2.1)
Equilibrium states x*(i, c) of equation (2.1) are determined by the solutions of the equation
x = i sin(x — c), (2.2)
and their stability is determined by the position of the roots of the characteristic equation
P(A; e1) = e1X + 1 — i cos(x* (i, c) — c) exp(—A) = 0. (2.3)
Accordingly, when |icos(x*(i,c) — c)| < 1 (> 1) and 0 < e1 ^ e0, where e0 is small, all the roots of Eq. (2.3) lie in the left open complex half-plane (Eq. (2.3) has roots belonging to the right open complex half-plane). In the former case, the equilibrium state x*(i,c) is asymptotically stable; and in the latter, unstable [4]. The boundary case is hcos(x*(i,c) — c)| = 1. This equality defines in the plane ¡, c the set of points of loss of stability of the equilibrium state x*(i,c) and bifurcation of periodic solutions from it. Reference [2] studies bifurcations of periodic solutions of Eq. (2.1) from equilibrium states x*(i,c) under the condition of stability loss i cos(x* (¡,c) — c) = —1. This type of loss of stability is the most typical one. However, two special cases remain unexplored. One of them is the case with c = 0, when the equilibrium state x*(i, 0) = 0 at the point = 1 loses stability under the condition cos(x*(¡*, 0)) = 1. Two stable equilibrium states x*+(¡, 0) > 0 and x*-(¡, 0) < 0 (|icos(x*± (l, 0))| < 1) bifurcate from x*(i, 0) = 0 at i > 1; at i = & 2.21 they simultaneously lose stability under the condition cos(x*±(i*, 0)) = —1. As will be shown below, a certain number (depending on the values of e1) of periodic solutions bifurcate in addition to these equilibrium states of x*(i, 0) = 0 with i > 1. The second is the case of the emergence of paired equilibrium states of Eq. (2.1). For every c, there is a sequence of values 0 < i1 (c) ^ i2(c) ^ ..., at which multiple roots x-(pk(c),c) = x+(ik(c),c), ik(c)cos(x±(ik(c),c) — c) = 1 appear in equation (2.2); when i > (c), they meet paired (stable and unstable) equilibrium states x-(i, c) and x+(i, c) of Eq. (2.1). As the parameter i increases further, the equilibrium state x-(i,c) loses stability at some i* and at the same time i* cos(x-(i*,c) — c) = —1.
Figure 1 shows diagrams of the equilibrium states of equation (2.1) for the parameter values 0 ^ i ^ 10, c = 0 and c = n/3 (for example). In the diagrams, stable equilibrium states are marked with a thick line and unstable ones with a thin line. The studied equilibrium states x*(i*(c),c), losing stability under the condition i*(c) cos(x*(i*(c),c) — c) = 1, are marked with dark circles, and under the condition i*(c) cos(x*(i*(c),c) — c) = —1 (for comparison) they are shown as light circles. As the parameter i increases further, the picture of equilibrium states of Eq. (2.1) does not change qualitatively — new paired (stable and unstable) equilibrium states appear, the nature of the loss of their stability being similar.
Next, we consider bifurcations of periodic solutions from equilibrium states in these critical cases using the ideology and results of [2]. The main emphasis is placed on identifying the features characteristic of the cases considered.
c = 0
10
Z = 7r/3
H 0
-5
-10
10
h 0
-5
0 2 4 6 8 10 6
-10
0 2 4 6 8 10
e
Fig. 1. Picture of the equilibrium states of the Ikeda equation for the cases c = 0 and c = n/3.
3. Bifurcation analysis of the stability loss of the equilibrium state x* 0) = 0 in the case c = 0
In the case c = 0, the loss of stability of the equilibrium state x*(i, 0) = 0 occurs at the point = 1. Set i = 1 + e2, \e2\ ^ 1 and write Eq. (2.1) in the neighborhood of x*(i, 0) = 0 as
eiy(t) + y(t) - (1 + e2)y(t - 1) + f (y(t - 1); e) = 0,
where f (y) = (1 + e2)y3/6 + o(y3) is an analytic function. Consider the linear part of Eq. (3.1)
eiy(t)+ y(t) - (1 + e2)y(t - 1) =0
and study the location of the roots of its characteristic equation
P(A; ei) = eiA + 1 - (1 + e2) exp(-A) = 0, A = 7 + ia.
(3.1)
(3.2)
(3.3)
Equation (3.3) has one real root Ao(e) = 70(e), which is the solution of the equation 7 + ln(1 + e1Y) = ln(1 + e2) in the domain -e-1 < 7, and is defined by the formula
Yo(e) = ln(1 + e2) - ln (1 + ei(ln(1 + e2) - ln(1 + ei(ln(1 + e2) - ...)))).
(3.4)
In the domain Im A > 0 we write Eq. (3.3) as eA(eiA + 1) = (1 +e2), which is equivalent to the sequence of equations
eA+in(i+ei A) = eina+^+mfc, ln £ = ln\z\ + i arg z - < arg z<n), k = 2,4,....
To determine the roots of Eq. (3.3), it is enough to consider the sequence of equations
A + ln (1 + eiA) = ln(1 + e2) + ink, k = 2,4,.... (3.5)
In the domain {(w,ti): w G WXQ = {w: —xo ^ Rew ^ .To, Imw ^ 7r, xq > 0 being some small fixed number}, \ei\ ^ e0}, we introduce the function
Al(w; ei) = - ln (1 + ei(w - ln (1 + ei(w - ln (1 + ei(w - ...)))))), "j^_RUSSIAN JOURNAL OF NONLINEAR DYNAMICS, 2020, 16(3), 437-451_
which will be continuous over a set of variables analytic over w £ Wx0 for every fixed e1 (|e11 ^ e0) and analytic over ti for every fixed w £ Wxo.
Each of Eqs. (3.5) has a unique solution defined by the formula
Ak (e) = ink + ln(1+ e2) + A1 (ink + ln(1 + e2); e1), k = 2,4,.... (3.6)
The proof of this repeats almost verbatim the proof of the corresponding statement [2] for odd k. Combining (3.4) and (3.6), we formulate the following statement.
Theorem 1. There is an eo > 0 such that when |e| < eo the entire set of roots of equation (3.3) is defined by the formula,
Afc(t) = ink + In (1 + t2) + X1(ink + In (1 + ¿2); ti), A_fc(t) = Afc(t), k = 0,2,4,.... (3.7)
In this case, the equality is performed uniformly with respect to k = 0,2,4,...
Ak (e) = Yk (e) + i(nk + ak (e)) = ink + ln (1 + e2) — ln (1 + e^ink + ln (1 + e2))) + 0(|e|2), Yk (e) = ln (1 + e2) — ln ((1 + e1 ln (1 + e2))2 + ej n2k2)/2 + 0(|e|2),
ak(e) = — arccos ((1 + e1 ln (1 + e2))/((1 + e1 ln (1 + e2))2 + e2xn2k2)1/2) + 0(|e|2). (3.8)
The stability of the solutions of Eq. (3.2) is determined by the sign of the functions Yk(e), which are analytic at the point e = 0 and have the radius of convergence of the corresponding series rk = 0(k-1), k > 0. At the same time
Yk(e) = e2 — e1(nk)2/2 — e1 e2 — e2/2 + o(|e|3), k = 0,2,4,...,
that is, for small e and the inequality e2 > e2(nk)2/2, the kth root of the characteristic equation (3.3) has a positive real part.
We now describe the structure of solutions of Eq. (3.1) in the neighborhood of the zero equilibrium state. The phase space of Eq. (3.1) is the space of continuous real-valued functions C[—1;0], the norm in which we define in the standard way and denote ||v(s)||c. By solving Eq. (3.1) (Eqs. (2.1)) defined at t ^ t0, with the initial condition y0(s) £ C[—1;0] we will understand the function y(t + s; e) (—1 ^ s ^ 0), which inverts Eq. (3.1) (Eq. (2.1)) into the identity and satisfies the initial condition y(t0 + s; e) = y0(s).
Denote by l2 the space of complex sequences of the form 2 = (z0,z2,z-2,z4,z-4,...), 2o e R, Zk £ C, k = 2,4,..., Z-k = Zk, \\z\\'l = Efc=o,2,4,... \zk\2 < By l\ we denote the subspace l2 of complex sequences z = (z0,z2,z_2,z4,z-4,...) for which ||z||=
= k=0,2,4,... ^k (e)|2 ^k |2 < TO and || £ k=0)±2)±4)... zk Ak (e)ek (s; e)||c < to, ek (s; e) = eXk (e)s/(1 + + e1 + Ak(e)). Further, s(r0) = {z £ I2, Mh2 < r0}, s1^) = {z £ l^, ||z||zi < r0}, S(R0) = = {y(s) £ C[—1;0], ||y(s)||c <R0}.
We introduce a system of differential equations in the space l2
zk = Ak(e)zk + dkikik3 (e)zkLzk2zk:i + Z*(z; e), k = 0, ±2, ±4,..., (3.9)
(ki ,k2,k3)enk
with the definition area of the right side of s1(r0), in which Q/k = {(k1 ,k2,k3): k1 ,k2,k3 = = 0, ±2, ±4,...,k1 + k2 + k3 = k,k1 < k2 < k3},
dkikik3 (e) = (e1 — (1+ e2)(e-AW3(e) — e-Xk(e))/(Ak1k2k3 (e) — Ak (e)))-1 (1 + e1 + Ak (e))fki k2k3 (e). _RUSSIAN JOURNAL OF NONLINEAR DYNAMICS, 2020, 16(3), 437-451_"j^
Afc1k2fc3(e) = Afc1 (e)+Afc2 (s)+\k3(e), fk^ks(e) = —Pki k^fcs(1+e2) eki(-1; e)efc2 (-1; e)efc3 (-1; e)/6, where pkkks = 1, if ki = k2 = fe?; Pkik2ks = 3, if ki = k2 = k3, either ki = k3 = k2, or k2 = k3 = k1; pk1k2k3 = 6, if k1 = k2 = k3; the smooth operator Z*(z; e) = (Z*(z; e), Z-i(z; e),...)(Z * (0; e) = 0): s^) ^ s(ro) and ||z* (z; e)||Z2 = o(||z||f:i).
The following statement defines the behavior of the solutions of Eq. (3 . 1) in the neighborhood of the zero equilibrium position.
Theorem 2. There are R0 ,r0,e0 > 0 such that at 0 < |e| < e0 any solution y(t + s; e) of Eq. (3.1) defined at t > 0 with the initial condition y0(s; e) G S(R0), at t ^ t0 (t0 ^ 3) and those t at which y(t + s; e) G s(r0) can be represented as
y(t + s; e) = (zk(t; e) + gk(z(t; e); e))ek(s; e), (3.10)
k=0,±2,±4,...
where the operator g(z; e) = (g0 (z; e), g2(z; e), g_2(z; e),...) (g(0; e) = 0): s1(r0) ^ s1(r0) and ||g(z;e)||zi = o(||z||zi), z(t;e) = (z0(t;e),z2(t;e),z_2(t;e),...) G s1(r0), and zk(t;e) (k = = 0, ±2, ±4,...) are solutions of the system of equations (3.9) with initial conditions zk(t0; e) uniquely determined by y(t0 + s; e). Conversely, if z(t; e) G s1(r0) at t ^ t0 is the solution of the system of equations (3.9), then the function y(t + s; e) defined according to (3.10) is the solution of Eq. (3.1). The operators Z*(z; e) and g(z; e) in (3.9), (3.10) are uniquely defined.
The proof of the theorem is analogous to the proof of the corresponding theorem of [2], where the case of odd k is considered.
The system of equations (3.9) will be called the normal form of Eq. (3.1) in the neighborhood of the zero equilibrium state, normalized to cubic terms inclusive. In [2], an algorithm is proposed to construct the normal form (3.9) to any order of expansion by zk.
Consider the analysis of the system of equations (3.9). Go to (3.9) to the polar coordinates in the planes (zk,z_k), k = 2,4,..., putting zk = pke%Tk (pk ^ 0, ^k=24 k2p| < to, —to < Tk < to) and denoting for uniformity z0 = p0. This yields the system of equations
Pk = Yk (e)pk + R(k) (p,t ; e) + R^ (p,T; e), k = 0,2,4,..., (3.11)
Tk = nk + Gk(e) + (Tk3)(p,T;e) + t^(p,t;e))/pk, k = 2, 4,..., (3.12)
where Yk(e) and Gk(e) are determined in (3.8), p = (p0,p2,p4,...), pk ^ 0, k = 2,4,..., t = (t2,t4,...), the functionals Rk3(•), R^O, Tk3(■), Tk*\-) are smoothly dependent on input
(3) (3)
variables and parameters 2n-periodic in Tj, with Rk)(■) and Tk )(■) being cubic forms on pk and R(-)|, TO = o(||p||fi).
The structure of the system (3.11)-(3.12) allows one to enter one "fast" variable and a countable number of "slow" variables. This cannot be done uniquely, but they are all connected by linear relations. We use the following method of introducing such variables. As the "fast" variable, we choose the variable t2 and consider "truncated" finite-dimensional systems, successively assuming in (3.9) zk = 0, k = ±4, ±6, ±8,..., then k = ±8, ±10,..., etc. In the first case, the "slow" variables are p0, p2. In the second case, the "slow" variables are p0, p2, p4, 62 = t4 — t2. In the next step, two new equations are added for the variables z6 and z_6, and the resonance monomials z_6z4z2 and z6z_4z_2 appear on the right-hand side of the equation for p0. As new "slow" variables we will choose p6, d4 = t6 — t4 — t2. In general, the next step adds two new "slow" variables pko, dka_2 = Tko — Tko_2 — t2. Continuing this process, we will carry out
a transition to the "slow" variables p = (po, p2, p4, •••), 0 = (02, 04,...) and "fast" t2, and the corresponding system of differential equations for their determination will have the form
pk = Yk(e)pk + 43) (p, 0; e) + R*k(p, 0,T2; e), k = 0, 2, 4,..., (3.13)
0k = ¿k(e) + ©k3)(p,0; e) + ©k(p,0,T2; e), k = 2,4,..., (3.14)
T2 = 2n + 02(e) + T(3) (p, 0; e)+ T k(p,0,T2; e), (3.15)
in which the functionals R^(•), R*k(•), ©k3)(•), t(3)(0, T*(-) are 2n-periodic in 0k, the rest of their properties and the properties of the function 5k(e) (5k(0) = 0) are defined by the properties of functions and the functionals included in (3.11)-(3.12). The phase space of the system of equations (3.13)-(3.15) is the product of the spaces l2 ® c = E ® R?, where l2 = = {p = (po, p2, p4,...)}, pk > 0, k = 2,4,..., ||p||| = £r=o pk < to, c = {0 = (02,04,...)}, \\0\\c = supk ^k| < to. The domain of definition of the right-hand side of (3.13)-(3.15) is the product of spaces ¿2 ® c0 = E0 ® R, where t? = {p G l2, ||p|||i = po + Sfc=2 k2pi < to},
co = {0 = (02, 04,...) G c, 0 < 0k < 2n}.
Let us introduce variables Z ^ 0 and —n/2 < 0 < n/2, putting
Z = (e? + |e2|)1^2, ei = Z cos 0, e2 = Z2 sin2 0 sign 0, (3.16)
in the domain {(e1,e2),e1 > 0, |e| < e0}. Substitute (3.16) into the system of equations (3.13)-(3.15) and select the equation for t2
T2 = 2n - Z cos 0 + Z2 cos 02 + O(Z3). (3.17)
Consider the main part (omitting R*k(•), ©k(•), T*(•)) of the equations of "slow" variables while normalizing pk ^ Zpk, t ^ t/Z2. This yields the system of equations
pk = Yk(0,Z)pk + R3 (p, 0; 0,Z) (Rii)(p, 0; 0, 0) = R^p 0)), k = 0,2,..., (3.18) 0k = 5t(^,Z) + ©k3)(p,0; ^,Z)(©k3)(p, 0; 0, 0) = ©k3)(p, 0)), k = 2, 4,..., (3.19)
where Yk(0,Z) = Yk(0,Z)/Z2, 5k(0,Z) = 5k(0,Z)/Z2 are, according to (3.7)-(3.8), continuous functions at n/2 ^ 0 ^ n/2, 0 ^ Z ^ Z0. At the same time
Yl(0, 0) = Yk(0) = sin2 0 sign 0 - n2k2 cos2 0/2, 5k (0,0) = 0. (3.20)
In (3.18)-(3.19), the properties of functions and functionals are determined by the properties of functions and functionals of the system (3.13)-(3.15). Consider the system of equations
pk = Yk(0)pk + R(i3)(p, 0), k = 0,2,..., (3.21)
0k = ©k3)(p,0), k = 2,4,.... (3.22)
Consider in l2 ® c the system of nonlinear equations
Yk (0)pk + Rf(p,0) = 0, k = 0,2,..., ©k3)(p,0) = 0, k = 2,4,..., (3.23)
assuming (p,0) G E0 = t? ® c0. In (3.23), the operator y(0)p = (YO(0)pO,Y2(0)p2,...) is considered based on (3.20) symmetrically extended by t?.
Let (p*(/),d*(/)) e T ® co solve the system of equations (3.23), i.e., the equilibrium state of the system of equations (3.21)-(3.22).
We introduce the infinite matrix
B(/) (Yk (№j + ^fjdPn dRH/dj \ (k . n2 k . 24 ) (324) B(/) = (3^ 1 „n( 3),^ )(ki, ji =0 ,2 , ...,k2, j2 = 2 ,4 ...) , (3.24)
V d @k2/dP31 d0k2;/d°h)
calculated at p*(/),d*(/), where 5kj is the Kronecker symbol. The matrix (3.24) defines the linear operator
B(/)v (v = (p,d)): El = Tl ® co — E = ® c, (3.25)
\\v\\E = \\p\\t2 + \\e\\c, \\v\\Ei = \\p\\j + \\e\\c.
In [2] it is shown that the operator (3.25) has a countable number of eigenvalues ¡ik(/) (B(/)vk (/) = nk(/)v(/)), which can be numbered in ascending order of their modules and any bounded region of the complex plane can belong only to a finite number of eigenvalues of finite multiplicity. The limit can only be the infinity point, with lim^.^ Re pk(/) = —to. The real parts ¡ik(/) determine the characteristic parameters of the system of equations (3.21)-(3.22) linearized at equilibrium (p*(/), 9* (/)) and determine the stability of the equilibrium state. If all characteristic indicators are negative, then the equilibrium state is exponentially stable; if there are positive indicators, it is unstable. The following theorem is valid, the proof of which repeats the proof of a similar theorem in [2], where another case of loss of stability of the equilibrium state is considered.
Theorem 3. Suppose that for some / the system of equations (3.21)-(3.22) has a solution (p*(/),&*(/)) e 12 ® Co (pk(/) > 0, k = 2,4,...), and the matrix constructed by this solution (3.24) defines the operator (3.25), which has no eigenvalues lying on the imaginary axis of the complex plane. Then there exists (0 > 0 such that for 0 < Z ^ (0 Eq. (3.1), in which el and e2 are defined according to (3.16), has aperiodic solution y*(t + s; e) = y*(t + s;/,Z) admitting the representation
m 2k-2
y*(7-2;/,Z) = Z(p0(/) + 2^2 P2k(/) cos (kr + ^ 9j(/)) + O(Z3), (3.26)
k=i j=i
where t2 = T2(t) is defined from Eq. (3.17). The solution y*(t + s;/,Z) will be asymptotically orbitally stable if all eigenvalues of the operator (3.25) lie in the left complex half-plane and unstable if m eigenvalues (counting multiplicities) belongs to the right complex half-plane. In the latter case, the periodic solution has an m-dimensional unstable invariant periodic manifold.
We present some results of the numerical analysis of stable equilibrium states (p*(/), d*(/)) of the system of differential equations (3.21)-(3.22), which are solutions of the system of nonlinear algebraic equations (3.23) and determine, according to (3.17)-(3.26), stable periodic solutions of Eq. (3.1) bifurcating from the zero equilibrium state. Calculations were carried out according to the following scheme. We considered a sequence of truncated finite-dimensional systems (3.21)-(3.22) in which k = 0,2,...,N. A finite-dimensional system of differential equations was numerically integrated under different initial conditions. For numerical integration, we used the software package [5], which allowed us to simultaneously build a large number of integral trajectories (with an automatic selection of the integration step) for randomly "thrown" (for example, according to a uniform law) initial conditions. There was a steady state of equilibrium (pN(/), 0*N(/)). This equilibrium state was refined in the system of equations (3.21)-(3.22)
0.02
£ 0.01
D i
0
0.01
0.02
0.03
0.04
0.05
Fig. 2. The D-partition pattern of the parameter domain e\,e2 of Eq. (3.1).
at k = 0.2,...,N + 2, etc. as long as Hp*N(0) - pN+2(0)||fi/WpN(0)H^ < A (A = 10 2 was
selected). Figure 2 gives a picture of D-partitions ([6]) of the parameter plane (e?,e2). It shows curves on which equation (3.3) has roots lying on the imaginary axis. These curves are the boundaries of the regions of D-partitions. At values e2 < 0 all roots of the characteristic equation (3.3) are in the left complex half-plane, at e G D? one root (3.3) is in the right complex half-plane, and the others are in the left complex half-plane, at e G D3 three roots (3.3) (one the real and two pairs of complex-conjugate roots) are in the right complex half-plane, and the others are in the left one, etc. The dashed line denotes the boundaries of the subsequent regions of D-partitions, similar to those described above. Let us fix the value of the parameter 0 = 1.0. According to (3.16), this value 0 corresponds to the curve in the region D?. Note first that the equation for p0 in (3.21) (for k = 0) has two more stable equilibrium states (±p0 (0),pk (0) > 0,pk (0) = 0,k = 2, 4,...), which are generated by stable equilibrium states xk+ 0) > 0 and xk-(^, 0) < 0 of equations (2.1) at n> 1. For the specified value 0 it was also possible to identify six stable equilibrium states (pk(0), 0k(0),pk(0) > 0,k = 2, 4,...). For each such equilibrium state at Z = 0.1, a periodic solution was constructed according to (3.17)-(3.26). All solutions have nonzero average values, i.e., the coordinate of pk(0) = 0. The periods of all solutions are close to one. The solutions have mutual antisymmetry — for each periodic solution ^*(t2; 0, Z) there is a solution — ^*(t2 ; 0, Z). Graphs of half of these solutions (for pk (0) > 0) are shown in Fig. 3. The bold line shows the solutions constructed by formulas (3.17)-(3.26), the thin line shows the graphs of exact solutions obtained by numerical integration of Eq. (3.1). The integration was carried out by the Euler method with constant step selection, at the same parameter values. The integration step was reduced to complete stabilization of the calculation results. The functions obtained according to (3.17)-(3.26) were chosen as initial values for integrating Eq. (3.1). As the parameter 0 decreases, the number of stable equilibrium states of the system of differential equations (3.21)-(3.22) decreases. This can happen in different ways: stable equilibrium states can merge with unstable ones and disappear, or two stable equilibrium states can merge with an unstable one to form one stable equilibrium state. As the parameter 0 increases, the number of stable equilibrium states of the system of differential equations (3.21)-(3.22) increases. This also happens in different ways — stable equilibrium states can bifurcate to form their own domains of attraction, separated by an unstable equilib-
rium state, or can arise in pairs "out of thin air" (a stable and an unstable equilibrium state). In addition, as the boundary of the domain of D-partition is passed through, an unstable periodic solution having a zero mean bifurcates from the zero state of equation (3.1). As the parameter ^ increases, this unstable periodic solution splits into three solutions: a stable one with zero mean and two unstable ones with nonzero mean. The solution with zero mean value corresponds to a stable equilibrium state (p*(^), d*(^), p0(^) = 0, pk> 0, k = 2,4,...) of the system of equations (3.21)-(3.22), from which, as ^ increases, two stable equilibrium states bifurcate simultaneously; these have ±p0(0), p0(^) > 0 and generate periodic solutions of Eq. (3.1) having mutual antisymmetry. The values ^ = 1.4, ^ = 1.49 and ^ = 1.51 correspond to the curves in the areas D3, D5 and D7, respectively. 8, 10 and 12 periodic solutions were found on the respective curves; the periods of these solutions are close to one, and the structure is identical to those shown in Fig. 3. Let us choose ^ = 1.51 and consider how the found periodic solutions will change with increasing Z. As Z decreases, the type of solutions does not change, only the amplitude of oscillations decreases, which tends to zero as Z ^ 0. As Z increases, each periodic solution passes through a series of period-doubling bifurcations and turns into irregular (chaotic) oscillations. For chaotic oscillations, the maximum characteristic index was calculated, and it turned out to be positive. In this case, there is a range of values of the parameter Z at which there are simultaneously twelve different chaotic oscillations (attractors). This phenomenon is called chaotic multistability. Numerical integration of Eq. (2.1) and calculation of the maximum characteristic index were performed using the software package [5]. Figure 3 shows one of 12 periodic solutions, its projection on the plane (x(t),x(t — 1)), as well as projections of phase portraits of solutions on the plane (x(t),x(t — 1)), which correspond to a transition of the selected periodic solution from regular oscillations to chaotic ones for the values Z = 0.2, 0.4, 0.6, 0.8.
Fig. 3. Periodic solutions of Eq. (3.1) for ^ = 1.0 and Z = 0.1.
4. Bifurcation analysis of the behavior of solutions
of the Ikeda equation at the birth of paired equilibrium states x- and
Consider the case where paired equilibrium states are born. In this case, for ¡j(c) = j = 1,2,..., Eq. (2.2) has two-fold roots x_ = x-(i*,c) = x+(i*,c) and thus cos(x±(i*,c) — — c) = 1. These roots at i > ¡ik correspond to the pair equilibrium states x-(i, c) and x+(i, c) of Eq. (2.1). In the case c = n/3, shown in Fig. 1, these will be points shown as dark circles with coordinates ^(c) & 2.42, x_ & 2.21, ¡2(c) & 6.73, x_ & —6.61, ¡3(c) & 8.91, x_ & 9.09.
We study the behavior of solutions of Eq. (2.1) in the neighborhood of equilibrium states x-(i,c) and x+(i, c). Replace i ^ ¡(1 + ¡). From (2.2) we have the representation x±(i*,c) as convergent series in powers of ¡1/2
xf(p, c) = x* ± fj,1/2 + tj,/(Sx*) =F p3/2V2(1 + 1/(9x2))/2 + Op2)
and the following equation:
¿t*(l + n) cos (xf(p, c) - c) = 1 =f avV^f1/2 +----
Denote t2 =x* \/2p1/2 + ... and write Eq. (2.1) in the neighborhood of the state x~(p,c), putting in (2.1)
x(t) = x-(i,c)+ y(t). (4.1)
As a result, we get the equation
eiy(t) + y(t) - (1 + e2)y(t - 1) + x*y2(t - 1)/2 + (1 + e2)y3(t - 1)/6 + ... = 0, (4.2)
where the dots denote terms having a higher order of smallness in y(t - 1). The linear part of Eq. (4.2) is of the form (3.2), and the corresponding characteristic equation is of the form (3.3), whose roots \k(e), k = 0, ±2, ±4,... are determined by Theorem 1.
We describe the behavior of solutions of Eq. (4.2) with initial conditions of y0(s; e) G S(R0). In l2 (defined above), consider the system of differential equations
¿k = Ak(e)zk + ^ dk!k2 (e)zkiZk2 + Z_(z; e), k = 0, ±2, ±4,..., (4.3)
2 k
(fcifc2)en2
with the domain of definition of the right side of s?(r0), in which Q = {(k1,k2): k1,k2 = = 0, ±2, ±4,..., k = ki + k2},
dkik2(e) = (ei - (1 + e2)(e-Akik2(e) - e-Ak(e)/(-Ak^ (e) + Ak(e)))-1 x (4 4)
(1 + e1 + eAk(e)) * fklk2(e),
Akik2(e) = Aki(e) + Ak2 (e), fkk (e) = -Pkk(1 + e2)eki (-1; e)ek2 (-1; e))x*/2, where pkik2 = 1 if k1 = k2; pkik2 = 2, if k1 = k2, and the operator, which is smooth in the set of variables, Z*(z; e) = (Zk(z; e),Z*_ 1 (z; e),...)(Zk(0; e) = 0): s1^) ^ s(r0) and ||Z_(z; e)^ = o(||z||2).
The behavior of solutions of Eq. (4.2) with initial conditions of y0(s; e) £ S(R0) is determined by the conditions of Theorem 2, in the formulation of which it is necessary to replace the expressions "Eq. (3.1)" with "Eq. (4.2)" and "system of equations (3.9)" with "system of equations (4.3)". The system of equations (4.3) will be called the normal form of Eq. (4.2) in the vicinity of the pair equilibrium states x_(p,c) and x+(p,c).
Transforming (4.3) to polar coordinates by setting zk = pkeiTk (pk ^ 0, ^k=2 4 k2p2k < to, —to < Tk < to, k = 2,4,...) and denoting zo = po, we obtain the following system of equations:
pk = Yk(e)pk + R(2)(p, t; e) + R*(p, t; e), k = 0, 2, 4,..., (4.5)
Tk = nk + ak(e) + (Tk2) (p,T; e) + T^(p,T; e))/pk, k = 2, 4,..., (4.6)
where Yk(e) and ak(e) are determined in (3.8), p = (p0,p2,p4,...), pk ^ 0, k = 2,4,...,
T = (T2,T4, . . .), the functionals Rf (•), R(k*)(-), Tf (•), T{k*)(-) are smoothly dependent on input
(2) (2)
variables and parameters 2n-periodic in Tj, with Rk '(•) and T( '(•) being quadratic forms in pk and \Rk(•)!, T (•)! = 0(||p||2j).
Instead of the variables p and t, we introduce one "fast" variable and a countable number of "slow" variables. As a "fast" variable, we choose t2. We will introduce "slow" variables by considering "truncated" finite-dimensional systems and sequentially assuming in (4.3) zk = 0, k = ±6, ±8,..., then k = ±8, ±10,..., etc. In the former case, p0 and p2, are chosen as "slow" variables. In the latter case, the "slow" variables are p0, p2, p4, 62 = t4 — t2. The next step adds two new equations for the variables z6 and z_6. In this case, no resonance monomials appear on the right-hand side of the equation for p0, a resonance monomial z_4z6 appears on the right-hand side of the equation for z2, and a monomial z2z4 appears on the right-hand side of the equation for z6. Changing to polar coordinates results in terms that depend on the expressions t6 — t4 — t2, —t6 + T4 + T2. As a new "slow" variable d4 we take the expression t6 — t4 — t2. As a result, we have two additional "slow" variables p6 ,d4. When we consider the next "truncated" system, resonant monomials z_6z8, z_4z8, z_2z8 will appear on the right-hand sides of the equations for z2, z4, z6, respectively, and when we change to polar coordinates, terms depending on the expressions t8 — t6 — t2 = d6 will appear. We have two new "slow" variables p8, d6. In general, the next step adds two new "slow" variables pko, dka_2 = Tko — Tko_2 — t2. Continuing this process, we will carry out a transition to the "slow" variables p = (p0, p2,...), d = (d2, d4,...) and to the "fast" variable t2, and the corresponding system of differential equations for their determination will have the form
pk = Yk(e)pk + Rik)(p, 0; e) + Rk(p, 9,T2; e), k = 0,2,4,..., (4.7)
9k = ¿k(e) + ©k2)(p,0; e) + ©k(p,0,T2; e), k = 2,4,..., (4.8)
T2 = 2n + a2(e)+ T(2)(p, 0; e) + T*(p,0,T2; e), (4.9)
in which the functionals R(2)(•), Rk(•), ©k2)(•), i(2)(-), T*(•) are 2n-periodic in 0k, and the rest of their properties and the properties of the function 5k(e) (5k(0) = 0) are defined by the properties of functions and functionals included in (4.5)-(4.6).
The phase space of the system of equations (4.7)-(4.9) is defined similarly to the way the phase space of the system of equations (3.13)-(3.15) is defined.
In the parameter space {(ei,e2),ei > 0, |e| < e0}, we introduce the variables ( ^ 0 and n/2 < ^ < n/2 according to (3.16) and substitute (4.7) into the system of equations (4.9). Consider its "main part" (omitting Rk(•), ©k(•), T*(•)) of the equations of "slow" variables,
normalizing pk ^ Z2pk, t ^ t/Z2. This yields the system of equations
pk = Yk (4,Z)pk + R(k) (p, 0; 4,Z) (R{kX.p, 0; 4, 0) = Rk)(p, 0)), k = 0,2,..., (4.10)
0k = 5* (^,C) + ©k2)(p,0; 4,Z), (©k2)(p, 0; 4,0) = ©k2) (p, 0)), k = 2,4,..., (4.11)
where Yk(4,Z) = Yk(4,Z)/Z2, 5k(4,Z) = 5k(4,Z)/Z2 (5k(4,0) = 0) are, according to (3.7)-(3.8), continuous functions at n/2 ^ 4 ^ n/2, 0 ^ Z ^ Z0. For Yk(4,0) = Yk(4), the expression (3.20) is valid. In (4.10)-(4.11), the properties of functions and functionals are determined by the properties of functions and functionals of the system (4.7)-(4.9).
Consider the system of equations
pk = Yk(4)pk + R(k)(p, 0), k = 0, 2, 4,..., (4.12)
0k = ©k2) (p,0), k = 2, 4,.... (4.13)
Theorem 4. Suppose that for some 4 the system of equations (4.12)-(4.13) has an exponentially stable or unstable equilibrium state (p*(4),0*(4)) £ E0 (pk(4) > 0, k = 2,4,...). In the latter case, m characteristic exponents (counting multiplicities) of the system (4.12)-(4.13) linearized on (p*(4),0*(4)) are positive. Then there exists (0 > 0 such that for 0 < Z ^ (0 equation (4.2), in which e1 and e2 are defined according to (3.16), has a periodic solution y* (t + s; e) = y* (t + s; 4, Z) with the same stability pattern. The dimension of the unstable manifold of the periodic solution is m. For the periodic solution the following formula holds:
m 2k_2
y*(T2; 4, Z) = Z2(p0(4) + 2 ^ p2k(4) cos (kT + ^ 0j(4)) + O(Z3), (4.14)
k=i j=i
where t2 = T2(t) is defined from Eq. (3.17).
Here we present some results of numerical analysis of the stable equilibrium states (p *(4),0*(4)) of the system of differential equations (4.12)-(4.13) for the cases c = n/3 and n* = ^1(c), and describe the dynamics of the periodic solutions (constructed according to (3.17)-(4.14)) of Eq. (4.2) with increasing parameter Z. The picture of ^-partitions of the parameter plane (e1,e2) of equation (4.2) is similar to that of ^-partitions of Eq. (3.1) and has the form shown in Fig. 2. Choose 4 = 1.51. For this value of 4, 6 stable equilibrium states of the system of equations (4.12)-(4.13) were revealed. Calculations were carried out according to the scheme described in the previous section. Note that in this case equilibrium states arise (and disappear) always in pairs — a stable and an unstable one. For these 6 equilibrium states, periodic solutions of Eq. (2.1) were constructed according to (4.1) and (3.17), (4.14) for the case Z = 0.1. Figure 5 shows graphs of these periodic solutions. The periods of all periodic solutions are close to one. The thin line shows graphs of exact solutions obtained according to (4.1) and of numerical integration of Eq. (4.2), in which e1 and e2 were determined by the expressions (3.16), and the approximate solutions obtained according to (3.17) and (4.14) were chosen as initial conditions. As the parameter Z increases further, each periodic solution passes through a series of bifurcations of doubling of the oscillation period and turns into a chaotic attractor. In this case, there is a range of values of the parameter Z at which 6 chaotic attractors of Eq. (4.2) (Eq. (2.1)) coexist, i.e., chaotic multistability takes place.
In conclusion, we note that the method of constructing self-oscillating solutions proposed in [2] and in this paper can be used to analyze the vibrational solutions of more complex mathematical models of optical systems [7, 8], taking into account the spatial structure of light waves, as well as various distributed models of reaction-diffusion systems [9].
Fig. 4. Periodic solution and projections of self-oscillating solutions of Eq. (3.1) onto the plane (x(t),x(t - 1)) for 0 = 1.51 and Z = 0.1,0.2,0.4,0.6,0.8.
Fig. 5. Periodic solutions of Eq. (2.1) in the neighborhood of the pair equilibrium state x* & 2.21 at ¡1 (n/3) & 2.42.
References
[1] Ikeda, K., Multiple-Valued Stationary State and Its Instability of the Transmitted Light a Ring Cavity System, Opt. Commun., 1979, vol.30, no. 2, pp. 257-261.
[2] Kubyshkin, E.P. and Moriakova, A. R., Features of Bifurcations of Periodic Solutions of the Ikeda Equation, Russian J. Nonlinear Dyn., 2018, vol. 14, no. 3, pp. 301-324.
[3] Ikeda, K. and Matsumoto, K., High-Dimensional Chaotic Behavior in System with Time-Delayed Feedback, Phys. D, 1987, vol. 29, nos. 1-2, pp. 223-235.
[4] Bellman, R. and Cooke, K.L., Differential-Difference Equations, New York: Acad. Press, 1963.
[5] Glyzin, D.S., The Program Package for the Analysis of Dynamic Systems "Tracer", Appl. no. 2008610548 dated Feb 14, 2008, Certificate of State Registration of the Computer Program no. 2008611464, Registered in the Register of Computer Programs 24.03.2008 (Russian).
[6] Neimark, Yu. I., D-Partitions of the Space of Quasipolynomials (to the Stability of Linearized Distributed Systems), Prikl. Mat. Mekh, 1949, vol.13, no. 4, pp. 349-380 (Russian).
[7] Chesnokov, S. S., Rybak, A. A., Stadnichuk, V. I., Optical Turbulence Modes in a Nonlinear Optical System with Time-Delayed Distributed Feedback, Atmos. Oceanic Opt., 2002, vol. 15, pp. 515-520.
[8] Budzinskiy, S. S., Razgulin, A.V., Rotating and Standing Waves in a Diffractive Nonlinear Optical System with Delayed Feedback Under O(2) Hopf Bifurcation, Commun. Nonlinear Sci. Numer. Simul, 2017, vol.49, pp. 17-29.
[9] Prizzi, M., Rybakowski, K.P., The Effect of Domain Squeezing Upon the Dynamics of Reaction-Diffusion Equations, J. Differential Equations, 2001, vol. 173, pp. 271-320.
