# Bifurcation and Chaos in a Price Game of Irrigation Water in a Coastal Irrigation District

Authors: Baogui Xin; Yuting Li
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

## Abstract

We propose a price game model of irrigation water in a coastal irrigation district. Then, we discuss the stability and codimension-two period-doubling (flip) bifurcation. Then, the MATLAB package Cl_MatContM is employed to illustrate its numerical bifurcations-based continuation methods. Lastly, the 0-1 test algorithm is used to compute the median value of correlation coefficient which can indicate whether the underlying dynamics is regular or chaotic.


---

## Body

## 1. Introduction

Water scarcity is one of the key problems affecting most countries in the world. With a burgeoning population, food price volatility, and climate change, water scarcity would also fuel future global conflict. Water scarcity is exacerbated by the indiscriminate discharge of industrial and municipal wastewater and is likely to affect the supply and demand of grain in the years ahead. Irrigation water availability is decreasing in many places where crop and plant production is taking place. Not only is there no set of efficient technique that can suddenly eliminate water scarcity, but also there is no optimal institutional arrangement for water, and rather it is critical to understand the potential contributions, facilitating conditions, and limitations of each [1, 2]. One of the important causes of water scarcity is that the demand exceeds a finite supply. All over the world, water regulations have historically focused on supply management. In fact, pricing mechanism may turn the tide against water scarcity by improving the water use efficiency [3–6]. Thereinto, the price game between water oligopolies is an important pricing mechanism. One of the simplest the price games is price game of irrigation water in a coastal irrigation district because there is little demand diversity of the irrigation water type among farmers. It will be considered in this paper.In recent years, a lot of research works [7–14] have shown that the game theory plays an important role in the economics and management field. Ji et al., Son et al., and Skoulidas et al. [15–17] studied the game model in an electric power market. Mu et al. [18, 19] analyzed the game model in a real estate market. Liu et al. [20] discussed the minority game in a financial market. Gkonis and Psaraftis [21] proposed a game model in the LNG market. Sun and Ma [22] presented a game model in Chinese cold rolled steel market. Sugawara and Omori [23] considered the duopoly in the Japanese airline market. Chung et al. [24] applied the game model into pollution permit markets. Ma and Zhang [25] build a price game in a property insurance market.Some references [26–32] have reported the complex dynamics of game model, such as bifurcation and chaos. Analyzing bifurcation and chaos is not an easy task for most of researchers. Fortunately, there are many powerful methods for us to study bifurcation and chaos, such as 0-1 test algorithm for chaos [33–38], MATLAB package MatCont series [39–43] for the bifurcation of discrete, and continuous dynamical systems.This paper is organized as follows. In Section2, a price game model of irrigation water in a coastal irrigation district is presented. In Section 3, the fixed points and their stabilities are studied. In Section 4, codimension-two period-doubling (flip) bifurcation is discussed. In Section 5, the 0-1 test algorithm and continuation methods are employed to validate the main results. Finally, conclusions in Section 6 close the paper.


## 2. A Pricing Game Model

In a coastal irrigation district, water mainly exists in ocean, rivers, lakes, or subterranean streams. But seawater and wastewater cannot directly be used for irrigation because the salt and pollutants will not allow the crops to grow. Generally speaking, groundwater can be directly used as irrigation water, but seawater and wastewater need to be pumped to a desalination plant and a wastewater treatment plant, respectively, and to be treated to be suitable for irrigation. Taking such factors as the freshwater scarcity, the high cost of wastewater treatment, and seawater desalination into consideration, wastewater treatment and seawater desalination have to rely on the help of government subsidies or tax breaks.Assumption 1.
As shown in Figure1, firms X, Y, and Z are the three water oligopolies of the irrigation water market in the coastal irrigation district. Firm X supplies irrigation water by directly pumping from rivers, lakes, or subterranean streams, and firm Y supplies irrigation water by wastewater treatment, and firm Z supplies irrigation water by seawater desalination.Figure 1
Schematic diagram of water supply in the coastal irrigation district.Assumption 2.
Firms X, Y, and Z compete with making different price of irrigation water in discrete-time periodst=0,1,2,….  Consider that pxt, pyt, and pzt represent, respectively, the irrigation water price of firms X, Y, and Z during period t=0,1,2,….Assumption 3.
The quantities, in which firms X, Y, and Z sell, respectively,Qxt, Qyt, and Qzt, are linear inverse demand functions determined by the following equations:
(1)Qxt=a-bpxt+d(pyt+pzt),Qyt=a-bpyt+d(pxt+pzt),Qzt=a-bpzt+d(pxt+pyt),
where a,b,d>0. The parameter d denotes the extent to which a firm's irrigation water is substituted by its rivals' water.Assumption 4.
The cost functions of firms X, Y, and Z have the following linear forms:(2)Cxt=c1Qxt,Cyt=c2Qyt,Czt=c3Qzt,
where parameters ci>0(i=1,2,3) are marginal costs of the firms X, Y, and Z, respectively.Assumption 5.
The profit functions of firms X, Y, and Z have the following forms:(3)Πxt=Qxtpxt-Cxt=(pxt-c1)(a-bpxt+d(pyt+pzt)),Πyt=Qytpyt+Qyts2-Cyt=(pyt+s2-c2)(a-bpyt+d(pxt+pzt)),Πzt=Qztpzt+Qzts3-Czt=(pzt+s3-c3)(a-bpzt+d(pxt+pyt)),
where s2 and s3>0 are the intensity coefficients of government support (such as subsidies and tax breaks) for the firms Y and Z, respectively.Assumption 6.
Firms X, Y, and Z always make the optimal price decision for the maximal marginal profit in every single period.The water prices of firms X, Y, and Z in period(t+1) are decided by solving the following optimization problem:
(4)pxt+1=argmaxpx∏x(pxt,pyt+1ex,pzt+1ex),pyt+1=argmaxpy∏x(pxt+1ey,pyt,pzt+1ey),pzt+1=argmaxpz∏x(pxt+1ez,pyt+1ez,pzt),
where pyt+1ex represents the expectation of firm X about the water price of firm Y during period t+1. Consider that pzt+1ex, pxt+1ey, pzt+1ey, pxt+1ez, and pyt+1ez may be explained by analogy.Assumption 7.
Each firm expects that its rivals' water price in period(t+1) will remain the same as in period (t).Thus,(5)pyt+1ex=pyt+1ez=pyt,pzt+1ex=pzt+1ey=pzt,pxt+1ey=pxt+1ez=pxt.The margin profits of firms X, Y, and Z in periodt are given, respectively, by
(6)∂Πxt∂pxt=a-2bpxt+d(pyt+pzt)+bc1,∂Πyt∂pyt=a-2bpyt+d(pxt+pzt)+b(c2-s2),∂Πzt∂pzt=a-2bpzt+d(pxt+pyt)+b(c3-s3).Let (6) be equal to 0; one can obtain the reaction functions of firms X, Y, and Z; that is, the optimal water prices are as follows:
(7)px*=a+d(py+pz)+bc12b,py*=a+d(px+pz)+b(c2-s2)2b,pz*=a+d(px+py)+b(c3-s3)2b.Assumption 8.
Firm X uses bounded rationality to make its price decisions with local information based on the marginal profits∂Πx/∂px and increase (decrease) its water prices in period (t+1) if the marginal profit is positive (negative) [8, 29, 44].The above adjustment mechanism of firm X has been called myopic by Dixit [45]. The dynamical adjustment mechanism of firm X can be written as follows:
(8)pxt+1=pxt+αpxt∂Πxt∂pxt,
where α>0 represents the adjustment speed of firm X.Assumption 9.
Firm Y is an adaptive decision maker and has adaptive expectations. Thus, its price decision in period(t+1) is mainly based on its reaction function and price in period (t).So, the price adjustment mechanism of firm Y can be written as follows:(9)pyt+1=(1-β)pyt+βpyt*,
where β>0 represents the adjustment speed of firm Y.Assumption 10.
Firm Z has simple rationality; that is, its price decision in period(t+1) is mainly based on its optimal reaction function in period (t).Thus, the price adjustment mechanism of firm Z can be expressed as follows:(10)pzt+1=pzt*.
So, the repeated price game of irrigation water in the coastal irrigation district has the following nonlinear form:
(11)pxt+1=pxt+αpxt(a-2bpxt+d(pyt+pzt)+bc1),pyt+1=(1-β)pyt+β2b(a+d(pxt+pzt)+b(c2-s2)),pzt+1=12b(a+d(pxt+pyt)+b(c3-s3)).
In what follows, we will focus on how the government supports intensities s2 and s3>0 which have an effect on the complex dynamics of the irrigation water price game.


## 3. Stability of Fixed Points

The fixed points of the system (11) satisfy the following algebraic system:
(12)px+αpx(a-2bpx+d(py+pz)+bc1)=0,(1-β)py+β2b(a+d(px+pz)+b(c2-s2))=0,12b(a+d(px+py)+b(c3-s3))=0.
By simple computation, one can obtain two fixed points E0=(px0,py0,pz0) and E1=(px1,py1,pz1), where
(13)px0=0,py0=ad+bdc3-bds3+2ab+2b2c2-2b2s24b2-d2,pz0=2ab+2b2c3-2b2s3+ad+bdc2-bds24b2-d2,px1=ad+2b(a+bc1)-bd(c1-c2-c3+s2+s3)2(2b+d)(b-d),py1=2b2(c2-s2)+bd(c1-s3-c2+c3+s2)+a(2b+d)2(2b+d)(b-d),pz1=2b2(c3-s3)+bd(c1+c2-c3-s2+s3)+a(2b+d)2(2b+d)(b-d).If the characteristic polynomial of a 3-order square matrix can be written as(14)P(λ)=(λ)3+a2λ2+a1λ+a0=0,
then one can get directly the following Lemma 11 from the Jury stability criterion.Lemma 11.
A necessary and sufficient condition that the characteristic polynomial of a matrixA=(aij)3×3 has all of its roots inside the unit circle is that
(15)1+a0+a1+a2>0,1-a0+a1-a2>0,1-a02>|a1-a0a2|,|a0|<1.


### 3.1. Stability ofE0

The Jacobian matrix of system (11) at the point E0 can be written as
(16)A(E0)=[1+α(a+d(y0+z0)+bc1)00βd2b1-ββd2bd2bd2b0].
Its characteristic polynomial can be written as
(17)P(λ)=(λ)3+A2λ2+A1λ+A0=0,
where A2=β-α(a+dy0+dz0+bc1)-2, A1=α(1-β)bc1+((y0+z0)(1-β)d+a(1-β))α+1-β-(βd2/4b2), and A0=(-1-αβ(a+d(y0+z0)+bc1))d2/4b2.From Lemma11, one can get the locally asymptotically stable region ΩE0(s2,s3) with respect to parameters (s2,s3) as follows:
(18)ΩE0(s2,s3)={1A02(s2,s3):1+A0+A1+A2>0,1-A0+A1-A2>0,1-A02>|A1-A0A2|,|A0|<1}.


### 3.2. Stability ofE1

The Jacobian matrix of system (11) at the point E1 has the following form:(19)A(E1)=[1+aα-4αbpx1+αdpy1+αdpz1+αbc1αdpx1αdpx1βd2b1-ββd2bd2bd2b0].Its characteristic polynomial can be written as
(20)P(λ)=(λ)3+a2λ2+a1λ+a0=0,
where a0=W5(W4(s2+s3)+W1+W2+W3), a1=W5(W11(s2+s3)+W6+W7+W8+W9+W10), a2=W5(W13(s2+s3)+W12), W0=c1-c2-c3, W1=β(αbW0-aα-1)d4, W2=(((1-4β)aα-β)-((2c2+2c3)β+W0)αb)bd3, W3=(2(1-2β)αb3c1+((1-2β)2aα+2β)b2)d2, W4=2αβb2d3-αb2d3+αβbd4, W5=1/(8b4-4b3d-4d2b2), W6=8α(β-1)b5c1, W7=(α(4-β)dW0+8(β-1)(aα-1))b4, W8=2d((2(β-1)(aα+1)-α(1+β)c1d)b3, W9=((W0(β+1)αd+2(β-(β+1)aα-2))b2d2, W10=d4β+(β-α(β+1))abd3, W11=d3αb2+4αb4d+αb2d3β-4αb4dβ, W12=8αb5c1+4(2aα+2β-αdW0-4)b4+4(2+aα-β)db3+4(2-β)b2d2, and W13=-4αb4d.It is obvious that the fixed pointE1 is locally asymptotically stable if and only if Lemma 11 holds. One can get the locally asymptotically stable region ΩE1(s2,s3) with respect to parameters (s2,s3) as follows:
(21)ΩE1(s2,s3)={a02(s2,s3):1+a0+a1+a2>0,1-a0+a1-a2>0,1-a02>|a1-a0a2|,|a0|<1}.


### 3.3. Parameter Basin with respect to(s2,s3)

Letα=0.36, β=0.2, a=6, b=2.5, c1=0.1, c2=0.3, and c3=0.4; a parameter basin with respect to the parameters (s2,s3) is shown in Figure 2, in which the two red regions correspond to ΩE0(s2,s3) and ΩE1(s2,s3), respectively, which are asymptotically stable state, the blue region denotes stable cycles of period two, the yellow region denotes chaotic state, and the white region denotes divergence state, as shown in Table 1.Table 1
Legend of the color coding for Figure2.Figure 2
Parameter basin versus the parameters(s2,s3).The regionsΩE0(s2,s3) and ΩE1(s2,s3) show that the price game for irrigation water will reach the Nash equilibrium by modulating limited times with random initial prices. Obviously, E0=(px0,py0,pz0) is a bounded equilibrium point [46], which indicates that px0=0; that is, free supply of irrigation water is an optimal strategy of firm X. But, in fact, it will never happen in the real word. Thus, E0 is not considered in the paper. And E1=(px1,py1,pz1) is a Nash equilibrium point, which is practical and feasible. So, we will continue to discuss E1 in the sections below.


## 3.1. Stability ofE0

The Jacobian matrix of system (11) at the point E0 can be written as
(16)A(E0)=[1+α(a+d(y0+z0)+bc1)00βd2b1-ββd2bd2bd2b0].
Its characteristic polynomial can be written as
(17)P(λ)=(λ)3+A2λ2+A1λ+A0=0,
where A2=β-α(a+dy0+dz0+bc1)-2, A1=α(1-β)bc1+((y0+z0)(1-β)d+a(1-β))α+1-β-(βd2/4b2), and A0=(-1-αβ(a+d(y0+z0)+bc1))d2/4b2.From Lemma11, one can get the locally asymptotically stable region ΩE0(s2,s3) with respect to parameters (s2,s3) as follows:
(18)ΩE0(s2,s3)={1A02(s2,s3):1+A0+A1+A2>0,1-A0+A1-A2>0,1-A02>|A1-A0A2|,|A0|<1}.


## 3.2. Stability ofE1

The Jacobian matrix of system (11) at the point E1 has the following form:(19)A(E1)=[1+aα-4αbpx1+αdpy1+αdpz1+αbc1αdpx1αdpx1βd2b1-ββd2bd2bd2b0].Its characteristic polynomial can be written as
(20)P(λ)=(λ)3+a2λ2+a1λ+a0=0,
where a0=W5(W4(s2+s3)+W1+W2+W3), a1=W5(W11(s2+s3)+W6+W7+W8+W9+W10), a2=W5(W13(s2+s3)+W12), W0=c1-c2-c3, W1=β(αbW0-aα-1)d4, W2=(((1-4β)aα-β)-((2c2+2c3)β+W0)αb)bd3, W3=(2(1-2β)αb3c1+((1-2β)2aα+2β)b2)d2, W4=2αβb2d3-αb2d3+αβbd4, W5=1/(8b4-4b3d-4d2b2), W6=8α(β-1)b5c1, W7=(α(4-β)dW0+8(β-1)(aα-1))b4, W8=2d((2(β-1)(aα+1)-α(1+β)c1d)b3, W9=((W0(β+1)αd+2(β-(β+1)aα-2))b2d2, W10=d4β+(β-α(β+1))abd3, W11=d3αb2+4αb4d+αb2d3β-4αb4dβ, W12=8αb5c1+4(2aα+2β-αdW0-4)b4+4(2+aα-β)db3+4(2-β)b2d2, and W13=-4αb4d.It is obvious that the fixed pointE1 is locally asymptotically stable if and only if Lemma 11 holds. One can get the locally asymptotically stable region ΩE1(s2,s3) with respect to parameters (s2,s3) as follows:
(21)ΩE1(s2,s3)={a02(s2,s3):1+a0+a1+a2>0,1-a0+a1-a2>0,1-a02>|a1-a0a2|,|a0|<1}.


## 3.3. Parameter Basin with respect to(s2,s3)

Letα=0.36, β=0.2, a=6, b=2.5, c1=0.1, c2=0.3, and c3=0.4; a parameter basin with respect to the parameters (s2,s3) is shown in Figure 2, in which the two red regions correspond to ΩE0(s2,s3) and ΩE1(s2,s3), respectively, which are asymptotically stable state, the blue region denotes stable cycles of period two, the yellow region denotes chaotic state, and the white region denotes divergence state, as shown in Table 1.Table 1
Legend of the color coding for Figure2.Figure 2
Parameter basin versus the parameters(s2,s3).The regionsΩE0(s2,s3) and ΩE1(s2,s3) show that the price game for irrigation water will reach the Nash equilibrium by modulating limited times with random initial prices. Obviously, E0=(px0,py0,pz0) is a bounded equilibrium point [46], which indicates that px0=0; that is, free supply of irrigation water is an optimal strategy of firm X. But, in fact, it will never happen in the real word. Thus, E0 is not considered in the paper. And E1=(px1,py1,pz1) is a Nash equilibrium point, which is practical and feasible. So, we will continue to discuss E1 in the sections below.


## 4. Codimension-Two Period-Doubling (Flip) Bifurcation

There are many bifurcation theories [47–49] that can be used into system (11), but the Kuznetsov bifurcation theory [50] is more effective to discuss the bifurcation in system (11). In what follows, we let α=0.36, β=0.2, a=6, b=2.5, c1=0.1, c2=0.3, and c3=0.4. The system (11) can be rewritten as follows:
(22)pxt+1=pxt+1.8pxt(1.25-pxt+0.11(pyt+pzt)),pyt+1=0.8pyt+0.022(pxt+pzt)-0.1s2+0.27,pzt+1=0.11(pxt+pyt)-0.5s3+1.4.Its Nash equilibriumE1= (1.096671926,1.186762016-0.4504504505s2, 2.580653598+0.4504504505s2). The Jacobian matrix of system (22) at the point E1 is
(23)Apb(E1)=[-0.9740.2170.2170.0220.80.0220.110.110],
which has a simple real eigenvalue λ1=-1 and other two eigenvalues λ2=0.0195 and λ3=0.807. From Figures 2, 3, and 4, one can find that a period-doubling bifurcation occurs when a simple real eigenvalue λ1=-1 crosses the boundary BE of the stability region ΩE1(s2,s3). That is, the critical parameters values s2 and s3 satisfy s2+s3=8.463662665 at the boundary BE of the stability region ΩE1(s2,s3).Figure 3
The Lyapunov exponent spectrum (blue, green, and red) and bifurcation ofpx (pinkish red) versus the parameters s2 when s3=0.4.Figure 4
The Lyapunov exponent spectrum (blue, green, and red) and bifurcation ofpx (pinkish red) versus the parameters s3 when s2=0.2.When the fixed pointE1 loses stability via a period-doubling bifurcation point, the restriction of system (22) to a one-dimensional center manifold at the critical parameter value can be transformed to the normal form as follows:
(24)Xn+1=-Xn+16b1Xn3+O(Xn4),Xn∈ℝ1,
where b1≠0 is called normal form coefficient [50], which is given by
(25)b1=16〈p,C(q,q,q)+3B(q,(I3-A)-1B(q,q))〉,
where I3 is the unit 3×3 matrix, Aq=-q, ATp=-p, A=Apb(E1), 〈q,q〉=〈p,q〉=1, 〈,〉 denotes the inner product, and the multilinear functions B and C are defined, respectively, by
(26)Bi(x,y)=∑j,k=1n‍∂2Xi(ξ,0)∂ξj∂ξk|ξ=0xjyk,i=1,2,Ci(x,y,z)=∑j,k,l=1n‍∂3Xi(ξ,0)∂ξj∂ξk∂ξl|ξ=0xjykzl,i=1,2.
For the system (22),
(27)q=(-0.994075101,0.01082788553,0.1081571934)T,p=(-0.9818661631,0.1055591755,0.2108811391)T,B(ξ,η)=(0.198(ξ1η2+ξ1η3+ξ2η1+ξ3η1)-3.6ξ1η100),C(ξ,η,ζ)=(000).
One can obtain
(28)B(q,q)=(-3.60430602600),B(q,(I3-A)-1B(q,q))=(-6.66976309200),C(q,q,q)=(000).
So, the critical normal form coefficient
(29)b1=3.274406>0,
which means that the period-doubling bifurcation at the fixed point E1 is supercritical.


## 5. Numerical Simulation

From Figures3 and 4, it can be observed that there is a very good agreement between the bifurcation diagram and the Lyapunov exponent spectrum. What is more, it can be find that the Lyapunov exponent spectrum and the bifurcation diagrams in Figures 3 and 4 well coincide with the parameter basin diagram in Figure 2, respectively. In this section, the numerical bifurcation and chaos will be employed to verify the above main results.


### 5.1. Numerical Bifurcation

In this subsection, based on continuation methods [51], we will discuss numerical bifurcations by using the MATLAB package Cl_MatContM [39–43].Firstly, we consider thatE1=(1.596211596,1.596211596,1.551166551) which is in the stable region BCDE of Figure 2. We do a numerical continuation of E1 with s2 free, and s3=0.4 fixed, as shown in Figure 5 and Table 2. Switchings at PD points of the second and fourth iterates are given in Figure 6.Table 2
Numerical continuation ofE1 with control parameter s2 as shown in Figure 5.
No.
Label
x
Normal form coefficient
Critical point
1
PD
(1.096672 −2.445180  1.051627  8.063663)
3.274406
e
+
00
2-cycle
2
PD
(0.830878 −0.392361  1.333279  4.045975)
1.868612
e
+
01
4-cycle
3
PD
(1.612503 −0.400875  1.248237  4.045975)
9.352670
e
+
01
4-cycle
4
PD
(0.702644  0.038991  1.393175  3.208909)
5.559248
e
+
02
8-cycle
5
PD
(1.019443  0.033503  1.378263  3.208909)
1.160791
e
+
02
8-cycle
6
PD
(1.594166  0.026410  1.281580  3.208909)
5.498331
e
+
02
8-cycle
7
PD
(1.727479  0.028661  1.315824  3.208909)
3.265258
e
+
03
8-cycleFigure 5
Continuation ofE1 in (s2,x)-space.Figure 6
Curves of fixed points of the 1st, 2nd, and 4th iterates.Secondly, from the fixed pointE1=(1.596211596,1.596211596,1.551166551), we do a numerical continuation of E1 with s2=0.2 fixed and s3 free, as shown in Figure 7 and Table 3.Table 3
Numerical continuation ofE1 with control parameter s3 as shown in Figure 6.
No.
Label
x
Normal form coefficient
Critical point
1
PD
(1.096672   1.096672 −2.490564  8.263663)
3.274406
e
+
00
2-cycle
2
PD
(0.830878  1.340060 −0.399142  4.245975)
1.868612
e
+
01
4-cycle
3
PD
(1.612503  1.331546 −0.484184  4.245975)
9.352670
e
+
01
4-cycle
4
PD
(0.702644  1.394355  0.037811  3.408909)
5.559248
e
+
02
8-cycle
5
PD
(1.019443  1.388868  0.022899  3.408909)
1.160791
e
+
02
8-cycle
6
PD
(1.594166  1.381774 −0.073784  3.408909)
5.498331
e
+
02
8-cycle
7
PD
(1.727479  1.384026 −0.039540  3.408909)
3.265258
e
+
03
8-cycleFigure 7
Continuation ofE1 in (s3,x)-space.In Tables2 and 3, the first three entries of x are the coordinate values of the fixed point E1, and the last entry of x is the value of the free parameters s2 or s3 at the corresponding bifurcation point. It is obvious that the normal form coefficient of the PD point is 3.274406, confirming (29). What is more, the detected bifurcation points in Figures 5 and 6 are in accordance with the statement in Figure 2. In addition, the label ‘‘PD’’ means that a period-doubling (flip) bifurcation occurs. In Table 2, the critical point ‘‘2-cycle’’ means that a stable 2-cycle is born when s2<8.063663, ‘‘4-cycle’’ means that a stable 4-cycle is born when s2<4.045975, and so on.


### 5.2. Numerical Chaos

In the above section, the Wolf algorithm [52] is employed to calculate the Lyapunov exponent spectrum shown in Figures 3 and 4, by which one can find chaos when a largest Lyapunov exponent is greater than 0. In this section, we will use a reliable and efficient binary test for the chaos (called ‘‘0-1 test’’) to detect chaotic attractors.


#### 5.2.1. The 0-1 Test Algorithm

The 0-1 test algorithm [33–38] can be described as follows.Consider a discrete set of measurement dataϕ(n) sampled at times n=1,2,3,…,N, where N is the amount of the data.Step 1.
Choose a random numberc∈(π/5,4π/5), and define the following new coordinates (pc(n),sc(n)):
(30)pc(n)=∑j=1n‍ϕ(j)cos(θ(j)),sc(n)=∑j=1n‍ϕ(j)sin(θ(j)),
where
(31)θ(j)=jc+∑i=1j‍ϕ(j),j=1,2,3,…,n.Step 2.
Define the mean square displacementMc(n) as follows:
(32)Mc(n)=limN→∞1N∑j=1N‍(pc(j+n)-pc(j))2+(sc(j+n)-sc(j))2,n∈[1,N10].Step 3.
Define the modified mean square displacementDc(n) as follows:
(33)Dc(n)=Mc(n)-(limN→∞1N∑j=1N‍ϕ(j))21-cosnc1-cosc,Step 4.
Define the median value of correlation coefficientK as follows:
(34)K=median(Kc),
where
(35)Kc=cov(ξ,Δ)var(ξ)var(Δ)∈[-1,1],
in which ξ=(1,2,3,…,ncut), Δ=(Dc(1),Dc(2),…,Dc(ncut)), ncut=round(N/10), and the covariance and variance are defined with vectors x,y of length q as follows:
(36)cov(x,y)=1q∑j=1q‍(x(j)-x-)(y(j)-y-),x-=1q∑j=1q‍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 quasiperiodic), 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 the Brownian-like (unbounded) trajectories imply that the underlying dynamics is chaotic.


#### 5.2.2. Application

We use the data setpx of the system (22) to implement the 0-1 test with s2 and s2, respectively. The new coordinates variables (p,s) are shown in Figures 8, 9, 10, 11, 12, and 13.Figure 8
K and the largest Lyapunov exponent versus s3∈[0,25] and s2=0.2.Figure 9
Plots versuss2=0.2 and s3=1 in new coordinates (p,s) space.Figure 10
Plots versuss2=0.1 and s3=1 in the original state space.Figure 11
K and the largest Lyapunov exponent versus s2∈[0,25] and s3=0.4.Figure 12
Plots versuss2=0.1 and s3=0.4 in new coordinates (p,s) space.Figure 13
Plots versuss2=0.1 and s3=0.4 in the original state space.Whens2=0.2 is fixed and s3 varies from 0 to 25 in increments of 0.05, one can get the diagram of K value as shown in Figure 8, which is consistent with the numerical simulation in the above sections. If s3=1, then K≈1, which means that the system is chaotic as shown in Figures 9 and 10.Similarly, fixings3=0.4 and varying s2 from 0 to 25 in increments of 0.05, one can get the diagram of K value as shown in Figure 11, which well coincides with the numerical simulation in the above sections. If s2=0.1, then K≈1, which means that the system is chaotic as shown in Figures 12 and 13.From Figures8 and 11, it can be observed that there is a very good agreement between the largest Lyapunov exponent and the median value of correlation coefficient K.


## 5.1. Numerical Bifurcation

In this subsection, based on continuation methods [51], we will discuss numerical bifurcations by using the MATLAB package Cl_MatContM [39–43].Firstly, we consider thatE1=(1.596211596,1.596211596,1.551166551) which is in the stable region BCDE of Figure 2. We do a numerical continuation of E1 with s2 free, and s3=0.4 fixed, as shown in Figure 5 and Table 2. Switchings at PD points of the second and fourth iterates are given in Figure 6.Table 2
Numerical continuation ofE1 with control parameter s2 as shown in Figure 5.
No.
Label
x
Normal form coefficient
Critical point
1
PD
(1.096672 −2.445180  1.051627  8.063663)
3.274406
e
+
00
2-cycle
2
PD
(0.830878 −0.392361  1.333279  4.045975)
1.868612
e
+
01
4-cycle
3
PD
(1.612503 −0.400875  1.248237  4.045975)
9.352670
e
+
01
4-cycle
4
PD
(0.702644  0.038991  1.393175  3.208909)
5.559248
e
+
02
8-cycle
5
PD
(1.019443  0.033503  1.378263  3.208909)
1.160791
e
+
02
8-cycle
6
PD
(1.594166  0.026410  1.281580  3.208909)
5.498331
e
+
02
8-cycle
7
PD
(1.727479  0.028661  1.315824  3.208909)
3.265258
e
+
03
8-cycleFigure 5
Continuation ofE1 in (s2,x)-space.Figure 6
Curves of fixed points of the 1st, 2nd, and 4th iterates.Secondly, from the fixed pointE1=(1.596211596,1.596211596,1.551166551), we do a numerical continuation of E1 with s2=0.2 fixed and s3 free, as shown in Figure 7 and Table 3.Table 3
Numerical continuation ofE1 with control parameter s3 as shown in Figure 6.
No.
Label
x
Normal form coefficient
Critical point
1
PD
(1.096672   1.096672 −2.490564  8.263663)
3.274406
e
+
00
2-cycle
2
PD
(0.830878  1.340060 −0.399142  4.245975)
1.868612
e
+
01
4-cycle
3
PD
(1.612503  1.331546 −0.484184  4.245975)
9.352670
e
+
01
4-cycle
4
PD
(0.702644  1.394355  0.037811  3.408909)
5.559248
e
+
02
8-cycle
5
PD
(1.019443  1.388868  0.022899  3.408909)
1.160791
e
+
02
8-cycle
6
PD
(1.594166  1.381774 −0.073784  3.408909)
5.498331
e
+
02
8-cycle
7
PD
(1.727479  1.384026 −0.039540  3.408909)
3.265258
e
+
03
8-cycleFigure 7
Continuation ofE1 in (s3,x)-space.In Tables2 and 3, the first three entries of x are the coordinate values of the fixed point E1, and the last entry of x is the value of the free parameters s2 or s3 at the corresponding bifurcation point. It is obvious that the normal form coefficient of the PD point is 3.274406, confirming (29). What is more, the detected bifurcation points in Figures 5 and 6 are in accordance with the statement in Figure 2. In addition, the label ‘‘PD’’ means that a period-doubling (flip) bifurcation occurs. In Table 2, the critical point ‘‘2-cycle’’ means that a stable 2-cycle is born when s2<8.063663, ‘‘4-cycle’’ means that a stable 4-cycle is born when s2<4.045975, and so on.


## 5.2. Numerical Chaos

In the above section, the Wolf algorithm [52] is employed to calculate the Lyapunov exponent spectrum shown in Figures 3 and 4, by which one can find chaos when a largest Lyapunov exponent is greater than 0. In this section, we will use a reliable and efficient binary test for the chaos (called ‘‘0-1 test’’) to detect chaotic attractors.


### 5.2.1. The 0-1 Test Algorithm

The 0-1 test algorithm [33–38] can be described as follows.Consider a discrete set of measurement dataϕ(n) sampled at times n=1,2,3,…,N, where N is the amount of the data.Step 1.
Choose a random numberc∈(π/5,4π/5), and define the following new coordinates (pc(n),sc(n)):
(30)pc(n)=∑j=1n‍ϕ(j)cos(θ(j)),sc(n)=∑j=1n‍ϕ(j)sin(θ(j)),
where
(31)θ(j)=jc+∑i=1j‍ϕ(j),j=1,2,3,…,n.Step 2.
Define the mean square displacementMc(n) as follows:
(32)Mc(n)=limN→∞1N∑j=1N‍(pc(j+n)-pc(j))2+(sc(j+n)-sc(j))2,n∈[1,N10].Step 3.
Define the modified mean square displacementDc(n) as follows:
(33)Dc(n)=Mc(n)-(limN→∞1N∑j=1N‍ϕ(j))21-cosnc1-cosc,Step 4.
Define the median value of correlation coefficientK as follows:
(34)K=median(Kc),
where
(35)Kc=cov(ξ,Δ)var(ξ)var(Δ)∈[-1,1],
in which ξ=(1,2,3,…,ncut), Δ=(Dc(1),Dc(2),…,Dc(ncut)), ncut=round(N/10), and the covariance and variance are defined with vectors x,y of length q as follows:
(36)cov(x,y)=1q∑j=1q‍(x(j)-x-)(y(j)-y-),x-=1q∑j=1q‍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 quasiperiodic), 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 the Brownian-like (unbounded) trajectories imply that the underlying dynamics is chaotic.


### 5.2.2. Application

We use the data setpx of the system (22) to implement the 0-1 test with s2 and s2, respectively. The new coordinates variables (p,s) are shown in Figures 8, 9, 10, 11, 12, and 13.Figure 8
K and the largest Lyapunov exponent versus s3∈[0,25] and s2=0.2.Figure 9
Plots versuss2=0.2 and s3=1 in new coordinates (p,s) space.Figure 10
Plots versuss2=0.1 and s3=1 in the original state space.Figure 11
K and the largest Lyapunov exponent versus s2∈[0,25] and s3=0.4.Figure 12
Plots versuss2=0.1 and s3=0.4 in new coordinates (p,s) space.Figure 13
Plots versuss2=0.1 and s3=0.4 in the original state space.Whens2=0.2 is fixed and s3 varies from 0 to 25 in increments of 0.05, one can get the diagram of K value as shown in Figure 8, which is consistent with the numerical simulation in the above sections. If s3=1, then K≈1, which means that the system is chaotic as shown in Figures 9 and 10.Similarly, fixings3=0.4 and varying s2 from 0 to 25 in increments of 0.05, one can get the diagram of K value as shown in Figure 11, which well coincides with the numerical simulation in the above sections. If s2=0.1, then K≈1, which means that the system is chaotic as shown in Figures 12 and 13.From Figures8 and 11, it can be observed that there is a very good agreement between the largest Lyapunov exponent and the median value of correlation coefficient K.


## 5.2.1. The 0-1 Test Algorithm

The 0-1 test algorithm [33–38] can be described as follows.Consider a discrete set of measurement dataϕ(n) sampled at times n=1,2,3,…,N, where N is the amount of the data.Step 1.
Choose a random numberc∈(π/5,4π/5), and define the following new coordinates (pc(n),sc(n)):
(30)pc(n)=∑j=1n‍ϕ(j)cos(θ(j)),sc(n)=∑j=1n‍ϕ(j)sin(θ(j)),
where
(31)θ(j)=jc+∑i=1j‍ϕ(j),j=1,2,3,…,n.Step 2.
Define the mean square displacementMc(n) as follows:
(32)Mc(n)=limN→∞1N∑j=1N‍(pc(j+n)-pc(j))2+(sc(j+n)-sc(j))2,n∈[1,N10].Step 3.
Define the modified mean square displacementDc(n) as follows:
(33)Dc(n)=Mc(n)-(limN→∞1N∑j=1N‍ϕ(j))21-cosnc1-cosc,Step 4.
Define the median value of correlation coefficientK as follows:
(34)K=median(Kc),
where
(35)Kc=cov(ξ,Δ)var(ξ)var(Δ)∈[-1,1],
in which ξ=(1,2,3,…,ncut), Δ=(Dc(1),Dc(2),…,Dc(ncut)), ncut=round(N/10), and the covariance and variance are defined with vectors x,y of length q as follows:
(36)cov(x,y)=1q∑j=1q‍(x(j)-x-)(y(j)-y-),x-=1q∑j=1q‍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 quasiperiodic), 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 the Brownian-like (unbounded) trajectories imply that the underlying dynamics is chaotic.


## 5.2.2. Application

We use the data setpx of the system (22) to implement the 0-1 test with s2 and s2, respectively. The new coordinates variables (p,s) are shown in Figures 8, 9, 10, 11, 12, and 13.Figure 8
K and the largest Lyapunov exponent versus s3∈[0,25] and s2=0.2.Figure 9
Plots versuss2=0.2 and s3=1 in new coordinates (p,s) space.Figure 10
Plots versuss2=0.1 and s3=1 in the original state space.Figure 11
K and the largest Lyapunov exponent versus s2∈[0,25] and s3=0.4.Figure 12
Plots versuss2=0.1 and s3=0.4 in new coordinates (p,s) space.Figure 13
Plots versuss2=0.1 and s3=0.4 in the original state space.Whens2=0.2 is fixed and s3 varies from 0 to 25 in increments of 0.05, one can get the diagram of K value as shown in Figure 8, which is consistent with the numerical simulation in the above sections. If s3=1, then K≈1, which means that the system is chaotic as shown in Figures 9 and 10.Similarly, fixings3=0.4 and varying s2 from 0 to 25 in increments of 0.05, one can get the diagram of K value as shown in Figure 11, which well coincides with the numerical simulation in the above sections. If s2=0.1, then K≈1, which means that the system is chaotic as shown in Figures 12 and 13.From Figures8 and 11, it can be observed that there is a very good agreement between the largest Lyapunov exponent and the median value of correlation coefficient K.


## 6. Conclusion

In this paper, we have proposed a nonlinear discrete price game model of irrigation water in a coastal irrigation district. Its stability and codimension-two period-doubling (flip) bifurcation are emphatically discussed. Based on continuation methods, its numerical bifurcations are analyzed by using the MATLAB package Cl_MatContM. Its numerical chaos is shown by means of the 0-1 test algorithm.



---

*Source: 408904-2013-05-12.xml*