Some remarks on an environmental defensive expenditures model

Manuscript submitted to Website: http://AIMsciences.org AIMS’ Journals Volume 00, Number 0, Xxxx XXXX pp. 000–000 SOME REMARKS ON AN ENVIRONMENTAL DEFENSIVE EXPENDITURES MODEL Tom´ as Caraballo Departamento de Ecuaciones Diferenciales y An´alisis Num´erico Universidad de Sevilla Apdo. de Correos 1160, 41080 Sevilla, Spain Renato Colucci Departament of Engineering, University Niccol`o Cusano, Via Don Carlo Gnocchi, 3 00166, Roma, Italy. Luca Guerrini Department of Management, Polytechnic University of Marche, Piazza Martelli 8, 60121, Ancona (AN), Italy (Communicated by the associate editor name) Abstract. In this paper, we consider the environmental defensive expenditures model with delay proposed by Russu in [15] and obtain different results about stability of equilibria in the case of absence of delay. Moreover we provide a more detailed analysis of the stability for equilibria and Hopf bifurcation in the case with delay. Finally, we discuss possible modifications of the model in order to make it more accurate and realistic. 1. Introduction. We consider a model introduced in [15] that deals with management of tourism in protected areas (PAs). The model is based on the interactions among visitors V, quality of ecosystem goods E, and capital K, intended as accommodation and entertainment facilities in PAs, which is described by the system:        · V=m1E+m2K−aV 2, · E=r¯ P−E−(b−cη)Vd, · K= (1 −η)Vd−δK, (1) where Vd=V(t−τ), τ≥0 is the time delay and ¯ Pis the maximum tolerable pollution stock P; m1, m2, a, b, c, η, δ, ¯ Pare strictly positive constants, 0 ≤η < 1 and 0 < r < 1. A detailed discussion of the model can be found in [15] but it is worth mentioning that a part, η, of total revenues is invested to defend the environmental resources while the other part,(1 −η), is used to increase the capital stock K. Also, we point out that the quality of ecosystem goods is measured as E=¯ P−P. The idea is to study the stability of the model, that is the sustainability of the tourism in PAs, This work has been partially supported by FEDER and the Spanish Ministerio de Econom´ıa y Competitividad project MTM2015-63723-P and the Consejer´ıa de Innovaci´on, Ciencia y Empresa (Junta de Andaluc´ıa) under grant 2010/FQM314 and Proyecto de Excelencia P12-FQM-1492. 1 2 SOME REMARKS ON AN ENVIRONMENTAL DEFENSIVE EXPENDITURES MODEL depending on the parameters involved (see also [1], [2], [6]-[10], for more general discussions on this the topic). For this reason the existence and stability of positive equilibria play a crucial role in this context. Our analysis of the equilibria is different from that in [15], in particular we find a necessary and sufficient condition for the existence of a positive equilibrium (V∞, E∞, K∞) inside R+×(0,¯ P)×R+that is (see next section for details): b−cη > 0,¯ Paδr −m2(1 −η)(b−cη)>0, which is in contrast with (a) of Proposition 1 in [15]. We also consider interesting to revise (in section 3 below) the bifurcation analysis for τ > 0 providing more detailed results with respect to those in [15]. Moreover we prove the existence of stability switches for the fixed point in the model with delay. We did not prove that the hypotheses are compatible with the positiveness of E∞, but still we consider the result interesting enough. In fact, in an attempt of using the model in real cases, the existence of stability switches makes the strategy of the managers of PAs more complicate in order to stabilize the system (see section 3 below). In section 5 of this work we discuss some critical issues about the model and propose some changes in order to improve it. 2. Fixed point and stability analysis with τ= 0.The steady states of system (1) in absence of delay are obtained as non-negative solutions of the algebraic system m1E+m2K−aV 2= 0, r ¯ P−E−(b−cη)V= 0,(1 −η)V−δK = 0.(2) As 0 ≤P≤Pand E=P−P, one has that P−E≥0.Hence, when b−cη ≤0,it is immediately seen that system (2) has no solution (we look for interior steady states, so we exclude the solution E=Pobtained for b−cη = 0). Let b−cη > 0.Then, the existence of fixed points for (1) is related to existence of solutions for the following second order equation in V aδrV 2−[(1 −η)rm2−(b−cη)δm1]V−m1δrP= 0.(3) A direct calculation shows that Eq. (3) possesses only one positive solution, independently of the sign of the coefficient of V. Lemma 1. Let b−cη > 0.System (1) has a unique fixed point (V∞, E∞, K∞),where V∞=(1 −η)rm2−(b−cη)δm1+q[(1 −η)rm2−(b−cη)δm1]2+ 4aδ2r2m1P 2aδr , E∞=P−(b−cη) rV∞and K∞=(1 −η) δV∞. If we want E∞>0 and b−cη > 0 we obtain the following necessary condition: ¯ Paδr −m2(1 −η)(b−cη)>0. The study of local stability of equilibrium solution is based on the localization on the complex plane of the eigenvalues of the Jacobian matrix of (1) evaluated at (V∞, E∞, K∞).The characteristic equation is det  −2aV∞−λ m1m2 −b+cη −r−λ0 1−η0−δ−λ = 0, namely λ3+Aλ2+Bλ +C= 0,(4) where A= 2aV∞+δ+r > 0, B = 2a(δ+r)V∞+δr −(1 −η)m2+ (b−cη)m1, SOME REMARKS ON AN ENVIRONMENTAL DEFENSIVE EXPENDITURES MODEL 3 and C= 2aδrV∞−(1 −η)rm2+ (b−cη)δm1 =q[(1 −η)rm2−(b−cη)δm1]2+ 4aδ2r2m1P > 0.(5) Proposition 2. Let b−cη > 0. The fixed point (V∞, E∞, K∞)of system (1) is locally asymptotically stable in absence of delay if and only if the following condition holds [2a(δ+r)V∞−(1 −η)m2+ (b−cη)m1] 2aV∞+ (δ+r)δr + (δ+r)22aV∞ −(1 −η)δm2+ (b−cη)rm1>0.(S) Proof. By Routh–Hurwitz’s condition, all roots of Eq. (4) have negative real parts if and only if A > 0, B > 0, C > 0 and AB > C. Since A > 0 and C > 0,it is clear that AB > C implies B > 0. Expliciting the terms A, B and Cin AB > C gives 4a2(δ+r)V2 ∞+ 2a(δ+r)2V∞+ (δ+r)δr + [−(1 −η)m2+ (b−cη)m1] 2aV∞ −(1 −η)δm2+ (b−cη)rm1>0. The statement follows by rewriting the previous inequality. Remark 3. From (5) we have (b−cη)m1>[(1 −η)rm2]/δ −2arV∞.Therefore, 4a2(δ+r)V2 ∞+ [−(1 −η)m2+ (b−cη)m1] 2aV∞ = [2a(δ+r)V∞−(1 −η)m2+ (b−cη)m1] 2aV∞ >2aδV∞+r−δ δ(1 −η)m22aV∞ and 2a(δ+r)2V∞−(1 −η)δm2+ (b−cη)rm1>2aδ2V∞+ 4aδrV∞+r2−δ2 δ(1 −η)m2. Consequently, if r≥δ, then condition (S) holds true and the fixed point is always stable. In [15] the claim (without proof) is that AB > C and as a consequence the fixed point is always stable. However, we only were able to prove the case r≥δleaving the case r < δ unsolved. 3. Stability analysis and Hopf bifurcation with τ > 0.The appearance of constant time delay does not affect equilibria. Hence, under the above conditions (V∞, E∞, K∞) is still the unique positive equilibrium of system (1). In this case the associated characteristic equation of the linearization of system (1) at the equilibrium is given by det  −2aV∞−λ m1m2 (cη −b)e−λτ −r−λ0 (1 −η)e−λτ 0−δ−λ = 0, i.e. λ3+a2λ2+a1λ+a0+ (b1λ+b0)e−λτ = 0,(6) where a2= 2aV∞+r+δ > 0, a1= 2a(δ+r)V∞+δr > 0, a0= 2aδrV∞>0 4 SOME REMARKS ON AN ENVIRONMENTAL DEFENSIVE EXPENDITURES MODEL and b1=−(1 −η)m2+ (b−cη)m1, b0=−(1 −η)rm2+ (b−cη)δm1. The stability of the equilibrium point (V∞, E∞, K∞) will change when the system under consideration has zero or a pair of purely imaginary eigenvalues. The former occurs when λ= 0 in (6), i.e. if a0+b0= 2aδrV∞−(1 −η)rm2+ (b−cη)δm1= 0,but we know this value to be positive from (5). Thus, λ= 0 is not a root of the characteristic equation (6). The latter deals with assuming that (6) has a purely imaginary root λ=iω, with ω > 0. Then, it follows from (6) that −ω3i−a2ω2+a1ωi +a0+ (b1ωi +b0) (cos ωτ −isin ωτ) = 0.(7) By separating real and imaginary parts in (7), we obtain ω3−a1ω=b1ωcos ωτ −b0sin ωτ, a2ω2−a0=b1ωsin ωτ +b0cos ωτ (8) Adding up the squares of Eqs. (8), we obtain ω6+a2 2−2a1ω4+a2 1−2a0a2−b2 1ω2+a2 0−b2 0= 0.(9) Setting z=ω2,Eq. (9) becomes h(z) = z3+pz2+qz +s= 0,(10) where p=a2 2−2a1= 4a2V2 ∞+ (r+δ)2+δr > 0, q=a2 1−2a0a2−b2 1= 4a2δ2+r2V2 ∞+δ2r2−[−(1 −η)m2+ (b−cη)m1]2,(11) s=a2 0−b2 0= 4a2δ2r2V2 ∞−[−(1 −η)rm2+ (b−cη)δm1]2.(12) Lemma 4. 1) Let q≥0and s≥0.Then Eq. (10) has no positive roots. 2) Let q≥0and s < 0or q < 0and s≤0.Then Eq. (10) has one positive root, say z0. 3) Let q < 0and s > 0.Set z∗=−p+pp2−3q 3.(13) Then Eq. (10) has no positive roots if h(z∗)>0,Eq. (10) has one positive root z∗if h(z∗) = 0, and it has two positive roots, say z−and z+,with z−< z+,if h(z∗)<0.Furthermore, h′(z∗) = 0, h′(z−)<0and h′(z+)>0. Proof. 1) is straightforward. For 2) and 3), we apply Descartes’ rule of signs that states that the possible number of positive roots (10) is equal to the number of sign changes in the sequence formed by the polynomial’s coefficients, or less than the sign changes by a multiple of 2.Finally, notice that h(0) = s, h(+∞) = +∞, h′(z) = 3z2+ 2pz +q, h′′(z) = 6z+ 2p > 0; h(z) is a convex function, and it has a unique minimum z∗= (−p+pp2−3q)/3 if q < 0 in the positive semiaxes. By using (8), we obtain sin ωτ =(b1a2−b0)ω3+ (a1b0−a0b1)ω b2 1ω2+b2 0 ,cos ωτ =b1ω4+ (a2b0−a1b1)ω2−a0b0 b2 1ω2+b2 0 . SOME REMARKS ON AN ENVIRONMENTAL DEFENSIVE EXPENDITURES MODEL 5 If Eq. (9) has the unique positive root ω0=√z0,then we can determine the sequence (j= 0,1,2, ...) τ0 j=         1 ω0 cos−1b1ω4 0+ (a2b0−a1b1)ω2 0−a0b0 b2 1ω2 0+b2 0+2jπ ω0 ,if M0≥0, 2(j+π) ω0−1 ω0 cos−1b1ω4 0+ (a2b0−a1b1)ω2 0−a0b0 b2 1ω2 0+b2 0,if M0<0, (14) with M0= (b1a2−b0)ω2 0+a1b0−a0b1,at which the characteristic equation (6) has a pair of purely imaginary roots λ=±iω0. Similarly, if Eq. (9) has the two positive roots ω±=√z±,with ω−< ω+, then Eq. (6) has purely imaginary roots λ=±iω±when τtakes the critical values (j= 0,1,2, ...) τ± j=           1 ω± cos−1b1ω4 ±+ (a2b0−a1b1)ω2 ±−a0b0 b2 1ω2 ±+b2 0+2jπ ω± ,if M±≥0, 2(j+π) ω±−1 ω± cos−1b1ω4 ±+ (a2b0−a1b1)ω2 ±−a0b0 b2 1ω2 ±+b2 0,if M±<0, (15) with M±= (b1a2−b0)ω2 ±+a1b0−a0b1. Proposition 5. Let λ(τ)be the root of (6) satisfying ℜ(τ0 j) = 0 (resp. Re(τ± j) = 0) and ℑ(τ0 j) = ω0 (resp. ℑ(τ± j) = ω±). Then λ=±iω0(resp. λ=±iω±)are simple roots of (6) when τ=τ0 j(resp. τ=τ± j)and the following transversality conditions hold dℜ(λ) dτ τ=τ0 j,ω=ω0 >0, resp. dℜ(λ) dτ τ=τ+ j,ω=ω+ >0and dℜ(λ) dτ τ=τ− j,ω=ω− <0!. (16) Proof. From (6), taking the derivative with respect to τ, we have 3λ2+ 2a2λ+a1+ [b1−τ(b1λ+b0)] e−λτ dλ dτ = (b1λ+b0)λe−λτ .(17) Let λ=iω, with ω∈ {ω0, ω+, ω−}, be a root of (6) that is not simple. Then this leads us to conclude from (17) that (b1iωk+b0)e−iωkτk j= 0.By separating real and imaginary parts, b1ωkcos ωkτk j−b0sin ωkτk j= 0, b1ωksin ωkτk j+b0cos ωkτk j= 0.(18) From (8) and as ωk>0,we see that (18) reduces to ω2 k−a1= 0, a2ω2 k−a0= 0.Thus, a1a2=a0 and so the contradiction 0 <[2a(δ+r)V∞+δr] (δ+r) + 4a2(δ+r)V2 ∞= 0.This proves the first part of the statement. Next, from (17) dλ dτ −1 =3λ2+ 2a2λ+a1 (b1λ+b0)λe−λτ +b1 (b1λ+b0)λ−τ λ. 6 SOME REMARKS ON AN ENVIRONMENTAL DEFENSIVE EXPENDITURES MODEL Using (6), a direct calculation yields sign (d(ℜλ) dτ τ=τk j,ω=ωk)=sign (ℜdλ dτ −1 τ=τk j,ω=ωk) =sign (ℜ3λ2+ 2a2λ+a1 (b1λ+b0)λe−λτ +b1 (b1λ+b0)λ−τ λτ=τk j,ω=ωk) =sign (3ω6 k+ 2 a2 2−2a1ω4 k+a2 1−2a0a2−b2 1ω2 k b2 1ω4 k+b2 0ω2 k) =sign 3ω4 k+ 2pω2 k+q =sign {h′(zk)}. Finally, notice that sign {h′(zk)}= +1 if ω=ω0or ω=ω+, while sign {h′(zk)}=−1 if ω=ω−. This completes the proof. Remark 6. If Eq. (9) has the unique positive root ω∗=√z∗, then sign {h′(z∗)}= 0 and the transversality condition does not hold. In case Eq. (6) has two positive roots ω−, ω+,with ω−< ω+, then crossing from left to right with increasing τoccurs whenever τassumes a value corresponding to ω+,and crossing from right to left occurs for values of τcorresponding to ω−.This implies that there exists a finite number of delayed intervals in which the equilibrium point is locally asymptotically stable, while unstable for the outside of the delayed ranges. Hence, the system dynamics switches from stable to unstable, and then back to stable when delay increases and crosses the critical delayed values. Based on the transversality conditions (16) and the Hopf bifurcation theorem (see [12]), we have the following results on the stability of the equilibrium of system (1). Theorem 7. Let b−cη > 0and assume condition (S)holds. Let h(z), p, q, z∗, τ0 0and τ± 0be defined as in (10),(11),(12),(13),(14) and (15),respectively. 1) If q≥0and s≥0or if q < 0, s > 0and h(z∗)>0,then the fixed point (V∞, E∞, K∞)of system (1) is locally asymptotically stable for all τ≥0. 2) If q≥0and s < 0or if q < 0and s≤0,then the fixed point (V∞, E∞, K∞)of system (1) is locally asymptotically stable for τ < τ0 0and unstable for τ > τ0 0.Furthermore, system (1) undergoes a Hopf bifurcation at (V∞, E∞, K∞)when τ=τ0 0. 3) If q < 0, s > 0and h(z∗)<0, then there exist stability switches for τ > 0.Furthermore, system (1) undergoes a Hopf bifurcation at (V∞, E∞, K∞)for a certain values of τ=τ± j (j= 0,1,2, ...)where a stability switch occurs. Proof. The proof of the theorem follows from Lemma 4 and Proposition 5. The expressions of the bifurcating values for the delays is given by (14) and (15). Remark 8. We note that the bifurcation analysis (together with the study of roots of eq. (10)) we propose here is more detailed and different with respect to that contained in [15]. 4. Numerical Simulations. In this section we will provide some numerical experiments to illustrate the results of the previous sections. SOME REMARKS ON AN ENVIRONMENTAL DEFENSIVE EXPENDITURES MODEL 7 4.1. A case of instability. We start with an example in which the fixed point is unstable. The parameters are fixed as follows: a= 0.1; b= 0.85; c= 1; δ= 1; η= 0.75; m1= 0.5; m2= 2; ¯ P= 25; r= 0.001, while the initial point is V(0) = E(0) = K(0) = 1. This choice satisfies the condition: b−cη = 0.1>0, while does not satisfy the condition of positiveness of E∞: ¯ Paδr −m2(1 −η)(b−cη) = −0.0475 <0. The fixed point is: V∞= 0.2524, E∞=−0.2397, K∞= 0.0631. In order to check the stability we need to compute the following quantities: A= 1.0515, B =−0.3985, C = 0.0496. Then, the stability condition is not fulfilled because AB −C=−0.4685 <0. The solution is represented in Fig. 1 below: the fixed point is not stable, the solution diverges to −∞. 0 2 4 6 8 10 12 14 16 18 -20 -15 -10 -5 0 5 Figure 1. Figure of experiment 1: the fixed point is not stable, the solution diverges 4.2. A limit cycle. We now consider an example in which there exists a limit cycle. In this experiment the parameters are set as follows: a= 0.1; b= 0.85; c= 1; δ= 0.2; η= 0.75; m1= 5; m2= 0.001; ¯ P= 25; r= 0.000001, while the initial point is V(0) = E(0) = K(0) = 0.1. This choice satisfies the condition b−cη = 0.1>0, 8 SOME REMARKS ON AN ENVIRONMENTAL DEFENSIVE EXPENDITURES MODEL while does not satisfy the positiveness of E∞: ¯ Paδr −m2(1 −η)(b−cη) = −2.45 ·10−5<0. The fixed point is: V∞= 2.5·10−4, E∞=−4.0978 ·10−6, K∞= 3.125 ·10−4 In order to check the stability we need to compute the following quantities: A= 0.2001, B = 0.4998, C = 0.1. Then, the stability condition is not fulfilled either because AB −C=−2.2472 ·10−5. The solution is represented in figures 2 below: the fixed point is not stable, there exists an attractive limit cycle. 1 0.5 0 -0.5 -1 -0.1 -0.05 0 0.05 0.1 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 0.15 0 100 200 300 400 500 600 700 800 900 1000 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 Figure 2. Figures for Experiment 2: stable limit cycle. 4.3. A case of stability with r < δ.We consider an example in which the interior fixed point is stable. The parameters are fixed as follows: a= 1; b= 0.73; c= 0.8; δ= 0.2; η= 0.9; m1= 5; m2= 0.001; ¯ P= 1; r= 0.1, and the initial point is V(0) = E(0) = K(0) = 1. This choice satisfies the conditions of existence of a positive fixed point: b−cη = 0.01 >0,¯ Paδr −m2(1 −η)(b−cη) = 0.02 >0, and the stability condition (S). This simulation suggests that stability is possible also for r < δ. In this case the fixed point is stable (see figure 3 below): V∞= 2.2361, E∞= 0.7764, K∞= 1.1180. and the solutions converges to it very fast. SOME REMARKS ON AN ENVIRONMENTAL DEFENSIVE EXPENDITURES MODEL 9 0 5 10 15 20 25 30 35 40 45 50 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 Figure 3. Stability of the fixed point for r < δ. 4.4. Delay model: stability change. In this example we illustrate the case in which the delay affects the stability of the fixed point. We consider the following values of the parameters: a= 0.5; b= 1; c= 0.8; δ= 0.2; η= 0.9; m1= 1; m2= 0.2; ¯ P= 4; r= 0.1; thus, the fixed point is V∞= 1.2102, E∞= 0.6113, K∞= 0.6051. The initial data has been chosen as E, V, K = (1,1,1) for t∈[−τ, 0], while the critical value of the delay parameter for stability is τ0 0= 9.7592. The hypotheses of 2) in theorem 7 are fulfilled, then we expect that the fixed point (V∞, E∞, K∞) is stable for τ < τ0 0and unstable otherwise. In Fig. 4 below we represent the solutions for τ= 8.6 and τ= 9.8 respectively: the fixed point looses stability and the solution diverges to −∞. time t 0 200 400 600 800 1000 1200 solutions -0.5 0 0.5 1 1.5 2 time t 0 10 20 30 40 50 60 70 solutions -4 -3 -2 -1 0 1 2 3 Figure 4. The solution of the system with delay and for τ= 8.6 and τ= 9.8 respectively.