Neimark–sacker Bifurcation Analysis and 0–1 Chaos Test of an Interactions Model between Industrial Production and Environmental Quality in a Closed Area

A discrete-time model is presented to describe the complex interaction between industrial production and environmental quality in a closed area. Its Neimark–Sacker bifurcation and chaos are discussed based on Wen's explicit Neimark–Sacker bifurcation criterion, Kuznetsov's normal form method and center manifold theory and Gottwald and Melbourne's 0–1 test algorithm. Numerical simulations are employed to validate the main results of this work.


Introduction
In the handicraft economy, manufacturing was often done in people's homes or small and rural shops by using hand tools or basic machines, and life for the average person was difficult, with meager incomes and malnourishment.In that period, the environment was a utopia, almost free from human disturbance, so the pollution risk to the environment was really quite negligible, and thus, the ecological environment's capability of self-adjustment and restoration was quite strong.In contrast, after the Industrial Revolution, people produced the bulk of their own food, clothing, furniture and tools by using new chemical manufacturing and iron production processes, improving the efficiency of water power, increasing the use of steam power and developing machine tools.Unfortunately, industrial production can produce many forms of pollution as follows: air pollution, light pollution, littering, noise pollution, soil contamination, radioactive contamination, thermal pollution, visual pollution, water pollution, plastic pollution and others, increasing the risk of various occupational hazards, such as asbestosis and pneumoconiosis.Though the industrial production brought about such a greater volume and variety of industrial goods and improved the living quality for many people, particularly for the middle and upper classes, the poor and working classes were lowly paid and struggling to improve their dangerous and monotonous working conditions.What is more, they could enjoy less blessings poured out from industrial production, but suffer more from industrial pollutants than the middle and upper classes.From this point of view, studying the interactions between industrial production and environmental quality would be helpful to promoting their harmonious development, on the one hand, and improving social equity and preventing the gap between the rich and the poor from widening, on the other hand.
Based on the interaction between product and environment, Salomone et al. [1][2][3] proposed their pioneering work, the product-oriented environmental management system (POEMS), a sustainable management framework.To address the complex relationship adequately, many researchers employed various appropriate theoretical frameworks to represent the nexus between industrial production and environmental quality.For example, Zhao et al. [4] proposed a plant-level aggregation method to estimate the relationship between production's spatial distribution and regional water environmental carrying capacity in a a small region.Aşıcı [5] explored the relationship between economic growth and the pressure on nature from the environmental sustainability perspective.Dinda [6], Copeland and Taylor [7] investigated the relationship between environmental degradation and economic development.Paraschiv [8] discovered relationships between the textile industry and sustainable development and conducted a Holt-Winters forecasting investigation for the eastern European area.Empirical methods were adopted in the above research.By using differential equations, Aliehyaei et al. [9] reported the results of exergy, economic and environmental analyses of simple and combined heat and power internal combustion engines.In this paper, difference equations and numerical simulations will be employed to discover complex interactions between industrial production and environmental quality in a closed area.
As a kind of dynamical bifurcation, Neimark-Sacker bifurcation [10][11][12] is a crucial phenomenon that has drawn considerable attention in many discrete-time systems, such as financial systems [13], investment competition models [14] and Hénon systems [15].The Neimark-Sacker bifurcation emerges with a closed invariant curve from a fixed point in discrete-time dynamical systems, when the fixed point changes stability via a pair of complex eigenvalues with unit modulus [10].We will focus on the existence, stability and direction of Neimark-Sacker bifurcation in the interaction system.Gottwald and Melbourne [16] firstly proposed the 0-1 test algorithm, a reliable and efficient binary test method for chaos, which is one of the simplest and most effective ones, though there have been many methods for detecting chaotic attractors.The 0-1 test algorithm has already been successfully implemented in various discrete or continuous systems, such as in [17][18][19][20][21][22][23][24][25].
The remainder of this paper is organized as follows.In Section 2, we formulate the model of the interaction between production and environment in a closed area.In Section 3, we analyze the stabilities of the real-valued fixed point of the model.In Section 4, we study the existence of Neimark-Sacker bifurcation by using Wen's Neimark-Sacker bifurcation criterion [11], and prove the stability and direction of Neimark-Sacker bifurcation by means of Kuznetsov's normal form method and center manifold theory [10].In Section 5, we detect the chaos by using the 0-1 test algorithm [16].Finally, the conclusion in Section 6 closes the paper.

Model
In this section, we focus on depicting a complex interaction between production and environment in a closed area by using nonlinear dynamics.For the purpose of our description of the model, let us first give some notations and assumptions as follows.

Notations
The following notations in Table 1 are used throughout this paper: Table 1.Notations.

Variables/Parameters Descriptions
x n the environmental quality index at period n; indicates the change of environmental quality, either improvement or deterioration; y n the production amount at period n; indicates the change of production amount, which represents a production decision, either an increase or a decrease; a the low threshold of environmental quality; c the high threshold of environmental quality; b the payment on the pollution emission right; δ the step size of a decision or measurement; a general, small parameter.

Assumptions
The following assumptions in Table 2 are available throughout this paper: Table 2. Assumptions.

Assumptions Descriptions Expressions
Assumption 1 The environmental quality has a positive linear impact on the production decision. ( Assumption 2 The last production amount has a negative linear impact on the current production decision due to congestion in the product market.

Model Formulation
Based on Assumptions 1-5, if we let c = 1 and readjust the dimensions of the production amount and the environmental quality index, then we can get the following equations: where a, b, δ > 0 As we know, in some cases, x changes more quickly than y; in other cases, x changes more slowly than y.Therefore, we can employ a general small parameter to proportionally synchronize x and y as follows: where > 0.
Case 1. : = 0.This means that y n+1 = y n , i.e., the production amount is a constant, for instance, because production resources are rigidly constrained or regulated by their government.
Case 2. : 0 < < 1.This means that x is a fast variable and y is a slow variable, for instance, because the environmental quality can projectively change in accordance to the producers' production amounts, whereas producers cannot change instantaneously their production amounts in a centrally-planned economy.
Case 3. : = 1.This means that x and y are completely synchronized.More prosaically, x and y have the same step size of decision or measurement.
This means that x is a slow variable and y is a fast variable, for instance, because the environmental quality changes gradually under government control, whereas producers can instantaneously and independently adjust their production amounts in a market economy.

Stability of the Fixed Points
The fixed points of system Equation (2) satisfy the following equations: which have only a real-valued fixed point E 0 = (x 0 , y 0 ), where: The Jacobian matrix of system Equation ( 2) at E 0 is given by: where Its characteristic equation can be written as: where: According to the relations between roots and coefficients of a quadratic equation, one can get the following proposition.Proposition 1.Let p(λ) = λ 2 + p 1 λ + p 2 .Suppose that p(1) > 0, λ 1 and λ 2 are two roots of p(λ) = 0.Then, the fixed point E 0 is: (i) locally asymptotically-stable if and only if p(−1) > 0 and p 2 < 1; (ii) a saddle if and only if p(−1) < 0; (iii) locally unstable if and only if p(−1) > 0 and p 2 > 1; (iv) non-hyperbolic if either p(−1) = 0 and p 1 0, 2 or p 2  1 − 4p 2 < 0 and p 2 = 1.

Existence of Neimark-Sacker Bifurcation
In order to discuss the existence of Neimark-Sacker bifurcation, an explicit criterion of Neimark-Sacker bifurcation needs to be introduced as follows.
Lemma 1. [11] For an n-th order discrete-time dynamical system, assume first that at the fixed point x 0 , its characteristic polynomial of Jacobian matrix A = (a i j ) n×n takes the following form: where a j = a j (µ, k), j = 1, • • • , n, µ is the bifurcation parameter and k is the control parameter or the other to be determined.Consider the sequence of determinants According to Lemma 1, for n = 2, we can get the following equalities and inequalities: By using the Mathematics software to solve Equations ( 6)- (10), the critical value of Neimark-Sacker bifurcation of system Equation ( 2) can be obtained as the two following expressions: when Thus, it follows from Equation (1) that the eigenvalues satisfy the condition (H1) in Lemma 1, that is Neimark-Sacker bifurcation occurs at the fixed point E 0 = (x 0 , y 0 ).

Direction and Stability of the Neimark-Sacker Bifurcations
In this section, we will use Kuznetsov's normal form method and center manifold theory [10] to investigate the direction and stability of the Neimark-Sacker bifurcations in the system Equation (2).Since the fixed point E 0 = (x 0 , y 0 ) is not origin O(0, 0), the E 0 need to be transformed to the origin by the following change of variables: x = x 0 + u, y = y 0 + v This transforms system Equation ( 2) into the following equivalent system: This system can be written as: where X n = (u n , v n ) T is the vector of the transformed system and J is the Jacobin matrix of system Equation ( 13) evaluated at the origin O(0, 0) as follows.
Additionally, the multilinear functions B : R 2 × R 2 → R 2 and C : R 2 × R 2 × R 2 → R 2 are defined respectively by: which take on the planar vectors ξ For the system Equation ( 13), The eigenvalues of the matrix J(O) are: where: Let q ∈ C 2 be a complex eigenvector of the matrix J corresponding to λ 1 given by Equation ( 15) and satisfying: Jq = e iθ 0 q Let p ∈ C 2 be a complex eigenvector of the transposed matrix J corresponding to λ 2 given by Equation ( 15) and satisfying: J T p = e −iθ 0 p Then, we can obtain: T where: T , to normalize p, let: where: We have p, q = 1, where ., .means the standard scalar product in C 2 : p, q = p1 q 1 + p2 q 2 .Therefore, the coefficients of the normal of the system Equation ( 13) can be obtained by the following formulas: g 20 = p, B(q, q) g 11 = p, B(q, q) g 02 = p, B(q, q) g 21 =2 p, B q, (I n − J) −1 B(q, q) + p, B q, e i2θ 0 I n − J Then, the direction coefficient of bifurcation of a closed invariant curve can be obtained by the following formula: Thus, the following theorem holds.

Numerical Example
In this section, we will justify the above analytic results by means of a bifurcation diagram, phase portrait and evolution series diagram.
In what follows, we let a = 0.5, b = 0.4 and the initial state (x s , y s ) = (0.4,0.6); then, system Equation ( 2) can be rewritten as the following form: which have a unique real-valued fixed point E 0 = (0.37, −0.03) and two complex conjugate fixed points E 1 = (0.57+ i0.87, 0.17 + i0.87), E 2 = (0.57− i0.87, 0.17 − i0.87).The E 1 and E 2 are always complex conjugate and will not be considered here.Figure 1 is the bifurcation diagram of x of the system Equation ( 18) with two parameters δ and .The bifurcation diagram illustrates the possible long-term values of the system Equation ( 18) as parameters δ and are varied.To produce Figure 1, varies from zero to 1.2 with an increment of 0.0024; and δ is sampled 12 values from zero to 1.2, and for x is used 300 data points after skipping 1000 transient data.Figure 1 can be regarded as a kind of superposition of bifurcation diagrams for 12 different values of δ.We plot Figure 1 to show bifurcation slices of xwith two parameters δ and .Each slice exhibits a bifurcation of x with the varying and a fixed δ. Figure 2 is the critical value curve of Neimark-Sacker bifurcation of system (18), which can be obtained from Equations ( 11) and (12).What is more, for combinations of parameters and δ in Figure 2, once crossing over the blue critical curve of Neimark-Sacker bifurcation and falling into the region on the left, the system Equation ( 18) must begin to undergo the Neimark-Sacker bifurcation.It is obvious that Figure 1 agrees well with Figure 2.   The real-valued fixed point E 0 is shown in Figure 3, and the bifurcation diagram is shown in Figure 4.The above two figures indicate that the system Equation ( 18) is asymptotically stable with = 0.4 and δ < 0.62773.
When one gives a small perturbation δ = 0.001, a sufficiently small real number, i.e., δ = δ * + δ = 0.62733 + 0.0001 = 0.62833, the system Equation ( 18) has a stable, closed invariant curve around the the fixed point (quasi-periodic solution), as shown in Figure 5. Furthermore, Figure 6 represents that the solution in the system Equation ( 18) asymptotically approaches a unique invariant closed circle, i.e., the Neimark-Sacker bifurcation is supercritical.
Step 1: Choose a random number c ∈ ( π 5 , 4π 5 ), then define the new coordinates (p c (n), s c (n)) as follows. where: Step 2: Define the mean square displacement M c (n) as follows: Step 3: Define the modified mean square displacement D c (n) as follows: Step 4: Define the median value of correlation coefficient K as follows: where: ), and the covariance and variance are defined with vectors x, y of length q as follows: x( j) var(x) = cov(x, x) Step 5: Interpret the outputs as follows: (1) K ≈ 0 indicates that the underlying dynamics is regular (i.e., periodic or quasi-periodic), whereas K ≈ 1 indicates that the underlying dynamics is chaotic.
(2) Bounded trajectories in the (p, s)-plane imply that the underlying dynamics is regular, whereas for Brownian-like (unbounded) trajectories, the underlying dynamics is chaotic.
In periodic dynamics, there are isolated values of c for which K c is large due to resonances, so it is suggested to use the median of the computed values of K c as K = median(K c ).In practice, 100 choices of c is sufficient to get various values K c versus c [21,22].In this work, we let N = 9500 and use 100 random c values in [π/5, 4π/5].
As mentioned in Subsection 3.2 and shown in Figure 2, if = 0.4 is fixed, the point ( , δ) = (0.4,0.6273) is the critical point of the Neimark-Sacker bifurcation of system Equation (18).Therefore, the region with δ ∈ [0, 0.6273) is stable, and the other region with δ ∈ (0.6273, 1.37] is unstable, where chaos may occur.With δ varying from zero to 1.37 in increments of 0.001, one can get the diagram of the K value and Lyapunov exponents shown in Figure 7.There is a very good agreement between Figure 7 and the bifurcation diagram shown in Figure 4. From Figures 4 and 7, one can find that chaos occurs when δ ∈ [1.35, 1.37], where K ≈ 1, which means that the system has chaotic movement.When δ = 1.36, Figure 8 presents a chaotic attractor in the original state space (x, y), and Figure 9 denotes the plots with the transformed coordinates (p, s), where Brownian-like trajectories denote that the system Equation ( 18) is chaotic., where K ≈ 1, which means that the system is chaotic.Let = 0.32; Figure 12 presents a strange attractor in the original state space (x, y), and Figure 13 shows the plot with the transformed coordinates (p, s), where Brownian-like trajectories indicate that the dynamics of system Equation ( 18) is undergoing chaos.

Conclusions
(i) We proposed a discrete complex interaction model about industrial production and environmental quality in a closed area, which can help us understand the above dynamical interaction mechanism from the environmental sustainability perspective.

(y n+1 − y n ) ∝ −y n Assumption 3
The payment on the pollution emission right has a negative linear impact on the production decision.(y n+1 − y n ) ∝ −b Assumption 4 The production amount has a negative linear impact on the change of environmental quality.(x n+1 − x n ) ∝ −y n Assumption 5 The change of environmental quality is affected simultaneously by the last environmental quality and deviations from its low and high thresholds.

Figure 1 .
Figure 1.Bifurcation slices of x with parameters δ and .
Critical value curve of Neimark−Sacker bifurcation

Figure 2 .
Figure 2. The critical value curve of Neimark-Sacker bifurcation with parameters δ and .

Figure 8 .
Figure 8. Plot of in the original state space (x, y) for system Equation (18) with = 0.4 and δ = 1.36.

Figure 12 .
Figure 12.Plot of in the original state space (x, y) for system Equation (18) with δ = 1.44 and = 0.32.