On the Proof of Existence of a Limit Cycle for the Prigogine Brusselator Model

This research was partially supported by the Research Grant of the Institute of Mathematics and IT of Academy of Sciences of Uzbekistan, No. FA-F1-F006 and the Research Grant (RUGS) of the University Putra Malaysia, No. 05-04-10-1005RU. Abstract Existence of a limit cycle for the dynamical system well-known as the Prigogine brusselator model is proved when parameters of the system take concrete values. The proof is conducted by the new method called discrete numerical (DN) tracking of trajectory combined with the Poincare-Bendickson theorem.


Introduction
Existence of a limit cycle for a differential equation is considerably nonlinear phenomena being of importance from the point of view of practice as well as a theory.Therefore, establishing the fact that the given dynamical system has a closed trajectory (cycle) is one of the essential problems of the Theory of Dynamical Systems.In spite of the rich arsenal of quantitative and qualitative tools of study for periodical solutions (Anosov, 1986, Hale, 1963, Cesari, 1959, Arnold, 1986, Lefever, 1971, Hassard, 1981, Andronov, 1966, Pliss, 1964, Krasnosel'skii, 1966, He, 2006), it should be admitted that in each concrete case the proof of the fact that the dynamical system has a cycle stays not a simple task.On the other hand, the majority of elaborations starting from the method of the first return function of Poincare up to recently ones are based on the assumption about existence of a cycle.However, even for the systems in the plane, though qualitative picture is thoroughly clear (Andronov, 1966), the problem of existence of a limit cycle is still a significant problem (Arnold, 1985, ch. 5, sec. 4).
In the present work, the method of discrete numerical tracking (for short, DN tracking) of trajectory, which is developed with the aim of establishing the existence of a closed trajectory when the right parts of a dynamical system are polynomials of low degree, is demonstrated on the example known as the Prigogine brusselator model (Lefever, 1971, Hassard, 1981): where a and b are positive parameters.The essence of the method was presented on the international conference "Pontryagin-2008" (Azamov, 2008).It was earlier applied to establish the existence of a closed trajectory for the system with the simplest quadratic nonlinearity (Azamov, 2009).
The bifurcation method enables to establish that as the point (a, b) goes through the parabola a 2 + 1 = b a limit cycle arises (Arnold, 1986, Hassard, 1981).Our purpose is to establish the existence of that for concrete values of a and b when the bifurcation method, generally, does not work.
Note that if then the system (1) has the unique equilibrium point P = (a, b/a) which is unstable focus.After shifting the origin to this point the system (1) takes the form where z = (x 1 , x 2 ) and f (z) = ( f 1 (z), f 2 (z)).Further constructions will be carried out for the values a = 1 and b = 2.01.
Remark.The reason for choice of these values for parameters is purely technical.This choice makes it possible to conduct the proof strictly taking into account the errors in performing arithmetical and comparison operations by means of computer.In fact, the method works for a certain range of values of a and b, for example, the method works for their values in the set of (2) cut by the condition 0.1 < a < 2. The more exact the approximately integrating method and results of the rounded arithmetical operations, the wider the set of values of a and b covered by the method.For the chosen values, it turns out, the Runge-Kutta method with the accuracy of second order is just enough (see section 2) allowing to reduce some details which are irrelevant for the purpose of the paper.
We shall construct a bounded region C which is invariant as t → +∞.Its boundary consists of arcs of the trajectories and transversal segments.

Main Results
Let z 0 (t) be the solution of the system (3) with the initial condition z 0 (0) = (0.32, 0), and z n = (x 1n , x 2n ) be the successive approximations obtained in applying the following numerical integrating scheme with the step h: Note that K 0 ⊂ K and dist(K 0 , ∂K) = 0.02.Clearly, the rectangle K contains the equilibrium point O = (0, 0).As the values a = 1, b = 2.01 satisfy the condition (2), then it is an unstable focus.Let N = 1055456, h = 10 −6 , and T = Nh = 1.055456.
Assumption 2. The solution z 0 (t) exists on the interval [0, T ] and satisfies the condition z 0 (t) ∈ K.
It should be noted that the existence theorems are not applicable to the system (1) on a priori given intervals, in particular, on [0, T ], which we are considering, not to mention the extendability of the solution z 0 (t) on [0, +∞).Therefore, it is not beforehand known that whether a solution z 0 (t) exists on the interval [0, T ], and if the condition z 0 (t) ∈ K, which is necessary to estimate the accuracy of the Runge-Kutta method, holds.
The estimate of accuracy of approximate integrating for the cubic dynamical systems being used here is based on the relation where are the Euclidian norms of a vector, a matrix, and a vector consisting of bilinear forms.
From now on, we consider z ∈ K.Note that validity of the inequality ( 5) rests upon the Assumptions 1 and 2.
Estimates of the type ( 5) are well-known (see, for example, Bakhvalov, (1973), p. 465, formula ( 9)).It should be noticed, that it was performed there only for one equation.Moreover, in definition of L (that plays the role of M * 1 ) the operation sup was taken only over the independent variable.As a result, L remains undefined when the solution is not known.The strict proof of the inequality ( 5) can be provided immediately for the system of differential equations by using Adamar's and Gronwall's lemmas.
We derive the relation (5).It is convenient for us to use the noninvariant norms (Kartan, 1971) If the dimensions of vectors z and f are equal to d, then Thus, in the case d = 2, the estimation ( 5) becomes worse a little The following inequalities can easily be verified for the system (3) As a result, in the case we deal with, the estimation (6) takes the following concrete form Let t ∈ [0, T ] and nh be the point of the net {nh| n = 0, 1, 2, ..., N} closest to the number t, i.e. n is the greatest integer that is less than or equal to t/h.The estimation M 0 < 2.88 implies Thus under Assumptions 1 and 2, the sequence z n allows to trace the real trajectory z 0 (t) to within the accuracy 4.2 However, in practice, it is impossible to deal with the exact values of z n because the formulas (4) are pure theoretical, while in calculating z n+1 based on z n one has to round results of arithmetical operations.As the expressions on the right part of the system (3) are polynomials of x 1 and x 2 , then there are no other kinds of errors of calculations.To take into account the rounding error we consider the sequence ζ n = (ξ n , η n ) along with z n .Consecutive terms of ζ n are calculated by the formula (4) with z n replaced by ζ n , of course, only by rounding, naturally, on computer.
We assume that there is possibility to realize arithmetical operations with rounding with an accuracy of at least 15 decimal places.More precisely, if s + , s − , and s × are values written in decimal form with 15 decimal places of the sum, difference and product, respectively, of two numbers a, |a| < 1, and b, |b| < 1, then each of the quantities does not exceed 1 2 • 10 −15 .Such a result is quite achievable on modern computers without any special algorithms to increase the accuracy.As it is obvious from (3) that to find ξ n+1 and η n+1 based on the values ξ n , η n for each 40 operations Proof.It can immediately be verified in process of calculations on computer that ζ n ∈ K 0 for all n = 0, 1, 2, ..., N. Therefore, dist(z n , K 0 ) < 0.01.As dist(∂K, K 0 ) = 0.02, then Assumption 1 holds.
Next, as z 0 ∈ K 0 , then the inclusion z 0 (t) ∈ IntK is true anyway for small positive values of t.We show that z 0 (t) can not reach the boundary ∂K of the rectangle K at any t ∈ (0, T ].Suppose, contrary to our claim, that z 0 (t) belongs to ∂K at some t * ∈ (0, T ], and n * h is the closest to t * point of the net {nh| n = 0, 1, 2, ..., N} from the left, i.e. n * = [t * /h].Then according to ( 8)-( 10) As z 0 (t * ) ∈ ∂K, then ζ n * K 0 , that contradicts the fact that ζ n ∈ K 0 for all n = 0, N.
holds for all t ∈ [0, T ], where t and n are connected with the relation n = [t/h].
Thus, we are able to follow the behavior of the trajectory z 0 (t) by means of the points ζ n .In particular, at N = 1055456 we have ζ N = (0.40909..., 0.56978...) where all the written digits are correct.
(We adopt the convention to write the rounded values of numbers with five decimal places being followed three points for unwritten digits.Absence of three points at the end of a decimal number means that its value is exact.) Set of the length 0.02 is parallel to the x 2 -axis and serves as the transversal for the system (3) (see Lemma 2 below).
Next, we let the trajectory z 1 (t) of the system (3) go out from the point z 1+0 .For this trajectory, we perform operations similar to those done in constructing the arc z 0 z 1−0 of the trajectory z 0 (t).Appropriate constructions and reasoning will be repeated 6 times and each time one will be convinced that they remain valid also for the trajectories going out from the points ζ i+0 , i = 0, 5 (ζ i+0 = ζ 0 for i = 0).The corresponding data are given in Table 1 and illustrated in Figure 1.
Duration of the initial five arcs has been exactly chosen equal to Nh = 1.055456 while the duration of the last arc equals the same value only approximately within ±10 −5 .The length of all segments S i , i = 1, 5 with the vertices ζ i−0 , ζ i+0 equals 0.02 sharply.But, the length of the first segment S 0 with the vertices ζ 0−0 , ζ 0+0 will be less than 0.02, where ζ 0−0 = ζ 6−0 .
Proof.It is enough to show that the segments S i , i = 0, 5 are transversals to the vector field (3).As an example, let us take the first segment 0.3 These estimations show that vector field (3) along S 0 is directed down, which means transversality of S 0 .The similar property can be verified for the segments S i , i = 1, 5.
Theorem.If a = 1, and b = 2.01, then the system (3) has at least one closed trajectory.
Indeed, a trajectory z * (t) which goes out from an inner point of the transversal [z 1−0 , z 1+0 ] at t = 0 does not leave the region C which contains the unique equilibrium point (0, 0) of the type of unstable focus.Therefore, by the Poincare-Bendickson theorem ω-limit set for the trajectory z * (t) must be closed trajectory.
Corollary 2. Given a = 1 and b = 2.01, the rectangle K contains at least one limit cycle of the system (3).
Indeed, it is known (Arnold, 1985, Ilyashenko, 1984) that an analytical dynamic system on the plane with regular equilibrium points may possess only finite number of closed trajectories.Particularly this property holds for the system (3).If it has the unique closed trajectory then it is obviously attractive.Now suppose that C 1 , C 2 , . . ., C k are closed trajectories of the system (3).They all must surround the focus O and situated one inside another.As C 1 is attractive from inside and C k is attractive from outside, then at least one of them must be attractive from both side.
In reality, the closed trajectory of the system (1) is, seemingly, unique.Unfortunately, the method of proof exposed here allows neither to determine nor estimate the number of closed trajectories).Numerical experiments show that the closed trajectory passes through the point x 0 ≈ 1.0829, y 0 = 2.01 and has the period τ ≈ 6.2866 (see Figure 1).
We have constructed "the bag of Bendicson" which border is sewed from 6 arcs of the system (3).If we take instead of T a value comparable with τ, for example, 6.2867, in the right part of the estimate ( 5) or ( 6) then the estimate becomes useless because of the multiplier e 2M 1 τ .