Self-oscillations in a certain Belousov–Zhabotinsky model

We consider the dynamic properties of a system of three differential equations known as the oreganator model. This model depends on four external parameters and describes one of the periodic Belousov–Zhabotinsky reactions. We obtain broad conditions for the parameters that ensure the existence of nonstationary steady-state regimes in oregonator model. With classical values of the parameters, the localization of the limit (at a long time) dynamics in the phase space has been improved. In fact, using numerical analysis, we significantly narrow the bounded region of the phase space containing the trajectories of the system. An iterative procedure is proposed for the approximate localization of closed trajectories (cycles) of the system on algebraic surfaces in R. A promising problem of theoretical substantiation of the numerical convergence of this procedure is posed.


Introduction
The problem of the existence of periodic regimes in problems of mechanics, chemistry, biology, economics, etc. is highly relevant. Let us consider in this connection one of the models of the Belousov-Zhabotinsky reaction [1][2][3][4], known in the literature as the "oregonator". This model corresponds to the system of ODEs with positive constants s, w, f , q. According to [1][2][3], system (1) describe oscillations in a real chemical reaction at s = 77.27, q = 8.375 · 10 −6 and appropriate constraints on w, f . Various aspects of this model were considered in [5][6][7]. An interesting work [8] demonstrates the possibility of applying the oregonator model to problems not related to chemistry. For q < 1, system (1) admits [2] a bounded positively invariant domain D ⊂ R 3 + containing a single stationary pointp. We find broad (different from the well-known [2]) conditions on the parameters s, w, q and f ≤ 1, under which a given point is unstable. Under these conditions, almost all solutions of system (1) go to nonstationary regimes when t → +∞. An * e-mail: liudmila.kondratieva@inbox.ru unstable invariant manifold M u (p) of a unstable pointp is two-dimensional, and a stable manifold M s (p) is one-dimensional. System (1) generates in D a smooth phase semiflow {Φ t } t≥0 and for each point p ∈ D\M s (p) the limit set ω(p) does not containp, and therefore has a nontrivial geometric structure. Acting as in [2], the existence of a closed trajectory (cycle) Γ ⊂ D of system (1) can be established using the so-called torus principle [9, § 18].
We analyze the instability of the system at q < 1, f ≤ 1. It is convenient to consider these two parameters as control ones. In principle, reducing (1) to a two-dimensional manifold M u (p) and applying the Poincaré-Bendixson theory can provide asymptotic orbital stability of the cycle Γ. A similar technique was applied in [10][11][12] to the problem of satellite motion, as well as to one biochemical model of dynamics in cellular processes.
In the case s = 77.27, q = 8.375·10 −6 , w = 41.86, f = 1, we significantly narrow the positively invariant domain of system (1). The new domain still contains a stationary pointpand an unstable manifold M u (p). Thus, we improve the localization of the limit (as t → +∞) phase dynamics of system (1). Numerical analysis using the Maple package demonstrates that the system enters the self-oscillation mode (M u (p) contains a cycle Γ) and this cycle "well placed" on some algebraic surface of the fifth order in R 3 .

Stationary point instability
As seen from (1), the positive octant R 3 + is positively invariant. It is easy to show (see [2]) that the parallelepiped has the same property. We find the only stationary pointp = (x,ȳ,z) ∈ R 3 + ,p =p( f, q), of system (1) from the relations The convex bounded positively invariant domain D ⊂ R 3 + contains at least one stationary point; therefore,p ∈ D. The identity qx ( Consider the Jacobi matrix Since γ > 0, α > 0, then according to the Routh-Hurwitz criterion the instability of a stationary point provides a condition γ > αβ, which therefore the condition β ≤ 0. The only stationary pointp ∈ D is non-degenerate, its Poincare index is 1, and the Jacobi matrix J ′ (p) with all permissible parameter values has in Re λ > 0 either two (with multiplicity) of eigenvalues or none. So, with β ≤ 0 an unstable manifold of the pointp two-dimensional.
The value β is the sum of the main 2 × 2 minors of the matrix J ′ (p) and .
In particular, the condition β ≤ 0 is fair if (1 − f ). Taking into account estimates (2), we obtain sufficient conditions for the instability of the stationary pointp ∈ D as It is convenient to record estimates (3.1), (3.2) in the form of q < 2 µ 2 0 at f = 1 and (3), the stationary pointp is a saddle with a onedimensional stable and two-dimensional unstable manifold.

Localization and visualization of limit dynamics
After normalization of variables x → s −1 x, y → sy, z → s −1 z, the system (1) acquires the form with a positively invariant domain We preserve the designations of the coordinates x, y, z and the stationary pointp when moving from the system (1) to the system (4). We assume f = 1 and try to reduce the domain Ω. Let us show that on the surface the vector fieldF of system (4) is directed to the domain. Indeed, if we take the field of normals on S 1 in the formn = 2s −1 (x + y), 2s −1 (x + y), −1 , then the scalar product Numerical analysis shows that on the surface S 2 the vector field of system (4) is directed to the domain Ω 2 : z < g 2 (x, y), i.e.F ·n| S 2 < 0. The plane z = z 0 , z 0 > s/4, intersects the surfaces S 1 and S 2 in the straight line x + y = √ s(z 0 − s/4) and in the hyperbola (x + 48)(y − 263) = s(333 − z 0 ), respectively. We accept z 0 = 251, then (also, numerically) it can be shown that on a part of the plane z = z 0 , (x, y, z 0 ) ∈ Ω 1 ∩ Ω 2 , the vector fieldF is directed to the domain Ω 3 : z < z 0 , i.e. F ·n = w(x − z) < 0 for the normaln = (0, 0, 1).
As we can see, the domain is positively invariant for system (4) and homeomorphic to a ball. The domain Ω 0 contains a single stationary pointp, its two-dimensional unstable manifold M u (p), and a stable (as numerical analysis shows) limit cycle Γ ⊂ M u (p). At the same time, domain Ω 0 is significantly less than Ω. Let the y + , y 0 and z + , z 0 be the maximum coordinate values of y and z in the domains Ω and Ω 0 , respectively. Then we have y + ≈ 4.6 · 10 6 , y 0 ≈ 789.9 and z + ≈ 1543.3, z 0 = 251. Let us pose the problem of approximate localization of the cycle Γ on some algebraic set (algebraic surface). For an arbitrary polynomial Q const in R 3 , we put If S (Q) ⊃ Γ and gradient ∇Q 0 on Γ then the vectorF lies in the tangent plane of the surface S (Q) for all points of the cycle, i.e. ∇Q ·F = 0 on Γ. Let us introduce a linear differential operator L : Q → ∇Q ·F. We have for an arbitrary solution x(t), y(t), z(t) of system (4) in Ω 0 , i.e. LQ is the derivative of the function Q(x, y, z) by virtue of (4). For the polynomial Q 0 , we set Q 1 = LQ 0 and, in general, Q k = LQ k−1 for k ≥ 1. At each iteration, the degrees of the polynomials are increased by one. The sets S (Q k ), k ≥ 1, contain stationary points of system (4) in R 3 . Let us try to choose an initial approximation Q 0 such that, with increasing k, the distances (deviations) from the cycle Γ to S (Q k ) decrease. Figure 1 shows the limit cycle Γ and the trajectory approaching it.
At the second iteration, we find the polynomial Q 2 = LQ 1 = ∇Q 1 ·F of degree 4 and the surface S (Q 2 ). At the third iteration, we obtain the polynomial Q 3 = LQ 2 = ∇Q 2 ·F of One can see (visually) that the deviations ρ 0 , ρ 1 , ρ 2 , ρ 3 are sequentially decrease and the cycle Γ "well placed" on the surface of the fifth order S (Q 3 ) ⊂ R 3 . It is interesting to find out under what conditions for the considered iterative procedure the deviations ρ k → 0 at k → ∞. This task looks very promising and may become the subject of further research.

Conclusion
The area of variation of parameters s, w, f, q > 0 is described, at which non-stationary limiting modes arise in system (1) and in its equivalent system (4). In the case s = 77.27, q = 8.375 · 10 −6 , w = 41.86, f = 1, the known positively invariant domain is substantially reduced to one containing a single stationary point and a stable limit cycle Γ. An iterative procedure for the approximate localization of the cycle Γ on a suitable two-dimensional algebraic surface in three-dimensional phase space is described.