Interactions Between Wealth and Natural Resources: A Nonlinear Ode Model

1. Introduction

This paper continues a research that we have carried out in our three previous works [1,2,3]. In all these papers we introduced low-dimension ODE systems intended to model the interactions between society and natural resources. The models in the quoted papers deal with four variables: x is the population, y the level of renewable resources, z the level of non-renewable resources, w the wealth produced using and depleting the resources. In this paper we introduced a model dealing only with $y,z,w$. This derives from the effort to study a model able to capture some characteristics of modern wealthy societies. In this kind of societies population is less or more stationary (or slowly decreasing) and the dynamics of interaction between society and resources is led by other variables. In [3] we introduced an ODE system in which the dynamics is led by non-renewable resources, while in the present paper we study a model in which the evolution is driven by the wealth w. The system is the following

$$\left\{\begin{array}{cc}\hfill y^{\prime }& =\gamma \lambda y-\gamma y^2-{\delta }_{1}yw,\hfill \\ \hfill z^{\prime }& =-{\delta }_{2}zw+k{\displaystyle \frac{{\delta }_3{w}^2}{{\delta }_3{w}^2+1}}z,\hfill \\ \hfill w^{\prime }& ={\delta }_{1}yw+{\delta }_{2}zw-{\delta }_3{w}^2-sw.\hfill \end{array}\right.$$

(1)

Here, $y,z,w$ are as we said above, while $\gamma ,\lambda ,k,{\delta }_i$ are positive coefficients. The equations for $y^{\prime }$ is the usual logistic equation, with a depletion term ${\delta }_{1}yw$. This means that the rate of depletion is given by the level of wealth in the society. In the same way, the equation for $z^{\prime }$ has a depletion term with rate of depletion ${\delta }_{2}w$. There is also, in the equation for z, a replenishment term, which says that human ingenuity is able to find new sources of non-renewable resources. The replenishment term depends only on w and is given through the function $f\left(t\right)=k{\displaystyle \frac{t^2}{t^2+1}}$, which is a saturation function of the type used in population dynamics. This term means that there is a replenishment of this kind of resources, but its rate is bounded. Finally, the equation for $w^{\prime }$ states that wealth increases thanks to the use of resources (renewable and non-renewable) and decreases because of consumption ($sw$) and because of the investment needed to replenish non-renewable resources (${\delta }_3{w}^2$). In the present model, differently from our previous papers, the rate of consumption depends on w, to be coherent with the idea to study a model in which the dynamics is driven by w.

We now briefly describes the main result of this paper (see also the conclusions, section 8). As we are interested on the trajectories with positive or non negative component, we first have proved that relevant trajectories remain positive. We then prove that none of the variables can diverge to $+\infty$, which means, in particular, that there can’t be indefinite growth in wealth. We then compute the equilibrium points of the system, and we study their stability. In some cases, due to the difficulty of this study, we get stability or instability results only when the depletion coefficients ${\delta }_i\to 0^{+}(i=1,2,3)$. We obtain that there are stable equilibria with positive values for all the variables, at least for some range of the parameters.

The search for mathematical models to investigate global social evolution and highlight the risks of collapse has nowadays a long tradition, starting at least with the celebrated work of the MIT group on the "Limits to growth": see the books [4,5,6], the update [7] and also [8]. The models introduced by the MIT group are high dimensional, and are studied through numerical computations. More recently some low dimensional model has been proposed for the same purposes. Our work started, in [1], from the HANDY model, which has been introduced in [9,10], and has then been studied by various authors (see [11,12,13,14,15,16,17,18]). After [1] we have changed in various way our models, so the problems studied in [2,3] and in the present paper can not be anymore considered examples of HANDY models. See also [19] for another model of interactions between society and nature. The main contribution of our research (including the present paper and the quoted ones) to this debate lies in the fact that, even assuming a moderate optimism on the replenishment of non-renewable resources, the aim of an indefinite growth in wealth seems impossible to attain, and also disruptive for the society. During our research, we have obtained this kind of results with different hypotheses on the main drivers of societal evolution, so they seem to have some kind of robustness (using the word in a non rigorous sense).

The paper is organized as follows: after the introduction, in section 2 we prove that the relevant trajectories are defined and non negative for all $t\ge 0$, and that they cannot go to $+\infty$. In section 3 we compute all the equilibria of system (1). In the following sections 4, 5, 6, we analyze the stability of the equilibria. In section 7 we show some numerical simulations and in section 8 we present a summary of the main results obtained.

2. General Results

We are interested in non-negative solutions of (1) hence, as in our previous papers, we will work in the positive cone $\mathcal{C}$ defined as

$$\mathcal{C}=\left\{(y,z,w)\in {\mathbb{R}}^3|y>0,z>0,w>0\right\},$$

or in its closure

$$\overline{\mathcal{C}}=\left\{(y,z,w)\in {\mathbb{R}}^3|y\ge 0,z\ge 0,w\ge 0\right\}.$$

Let $F:{\mathbb{R}}^3\to {\mathbb{R}}^3$ be the vector field of the system (1). Note that $F\in C^{\infty }\left({\mathbb{R}}^3\right)$. Thus, we have a standard result of existence and uniqueness:

Proposition 1.

For any $t_0\in \mathbb{R}$ and $X_0\in {\mathbb{R}}^3$, there exists a unique solution to the Cauchy problem

$$\left\{\begin{array}{cc}\hfill & X^{\prime }\left(t\right)=F\left(X\right),\hfill \\ \hfill & X\left(t_0\right)=X_0,\hfill \end{array}\right.$$

(2)

defined on a maximal interval $J=(\alpha ,\beta )$ with $-\infty \le \alpha <t_0<\omega \le +\infty$.

If $X_0=(y_0,z_0,w_0)$, we assume $y_0>0,z_0>0,w_0\ge 0$, that is, we assume that the social-nature evolution we study starts with a positive amount of renewable and non-renewable resources, while the starting social wealth is positive or zero. However, the next proposition shows that the case $w_0=0$ can be dealt with in a very easy way.

Proposition 2.

Let $X_0=(y_0,z_0,0)$ with $y_0>0,z_0>0$. Let $y\left(t\right)$ be the solution of the logistic equation

$$\left\{\begin{array}{cc}\hfill y^{\prime }\left(t\right)& =\gamma y(\lambda -y)\hfill \\ \hfill y\left(t_0\right)& =y_0.\hfill \end{array}\right.$$

(3)

Then the function $X\left(t\right)=(y\left(t\right),z_0,0)$ is the unique solution of the Cauchy problem:

$$\left\{\begin{array}{cc}\hfill & X^{\prime }\left(t\right)=F\left(X\right),\hfill \\ \hfill & X\left(t_0\right)=X_0.\hfill \end{array}\right.$$

(4)

Proof. 

The proof follows from a straightforward computation.    □

Proposition 2 shows that, when the initial value $w_0$ is zero, the solution is known, and no further analysis is necessary. For this reason, from now on, we assume $w_0>0$, that is, we will study system (1) with initial condition $X_0=(y_0,z_0,w_0)\in \mathcal{C}$. We will also assume $t_0=0$.

Proposition 3.

Let $X\left(t\right)=\left(y\right(t),z(t),w(t\left)\right)$ be a solution of (2) with $t_0=0$ and $X\left(0\right)\in \mathcal{C}$. Then $X\left(t\right)\in \mathcal{C}$ for all $t\in (\alpha ,\omega )$.

Proof. 

From Proposition 2, we know that $w\left(t\right)>0$ for all $t\in (\alpha ,\omega )$. On the other hand, the two equations for $y^{\prime }\left(t\right)$ and $z^{\prime }\left(t\right)$ have the zero constant $\xi \left(t\right)=0$ as a solution for all $t\in \mathbb{R}$, hence this constant cannot be crossed by another solution, ensuring that $y\left(t\right)$ and $z\left(t\right)$ remain positive for all $t\in (\alpha ,\omega )$.    □

Proposition 4.

Let $X_0\in \mathcal{C}$ and $X\left(t\right)=\left(y\right(t),z(t),w(t\left)\right)$ be the solution of (2) with $t_0=0$. Then $y\left(t\right)$ is bounded on $[0,\omega )$.

Proof. 

We distinguish two cases

$i)$$y\left(t\right)\ge \lambda$ for all $t\in [0,\omega )$: in this case $y^{\prime }\left(t\right)<0$ for all t, hence $\lambda \le y\left(t\right)\le y_0$ for all t, hence $y\left(t\right)$ is bounded.

$ii)$$y\left(t_0\right)<\lambda$ for some $t_0\in [0,\omega )$. Let us define

$${\tau }_0=sup\{t\in [t_0,\omega )|y\left(s\right)<\lambda \forall s\in [t_0,t)\}.$$

If ${\tau }_0=\beta$ we have $y\left(t\right)<\lambda$ for all $t\in [t_0,omega)$, so $y\left(t\right)$ is obviously bounded in $[0,\omega )$. If we assume ${\tau }_0<\omega$ then , by continuity, it is easy to obtain ${\tau }_0>t_0$, $y\left(t\right)<\lambda$ for all $t\in [t_0,{\tau }_0)$, $y\left({\tau }_0\right)=\lambda$, so that $y^{\prime }\left({\tau }_0\right)\ge 0$. On the other hand, using the equation for $y^{\prime }$, we get

$$y^{\prime }\left({\tau }_0\right)=\gamma y\left({\tau }_0\right)(\lambda -y\left({\tau }_0\right))-{\delta }_{1}y\left({\tau }_0\right)w\left({\tau }_0\right)=-{\delta }_1\lambda w\left({\tau }_0\right)<0.$$

This is a contradiction, so this case is ruled out, hence ${\tau }_0=\omega$ and y is bounded.

   □

Proposition 5.

Under the same assumptions, we have $\omega =+\infty$.

Proof. 

We argue by contradiction. Assume $\omega <+\infty$. Consider the variable $z\left(t\right)$, which satisfies $z^{\prime }\left(t\right)<kz\left(t\right)$ for all $t\in [0,\omega )$, hence $z\left(t\right)\le z_0{e}^{kt}$ for all $t\in [0,\omega )$. Consequently, $z\left(t\right)\le z_0{e}^{k\omega }$ for all $t\in [0,\omega )$, so it is bounded in $[0,\omega )$. Similarly, for $w\left(t\right)$, recalling that $y\left(t\right)$ is bounded, we have $w^{\prime }\left(t\right)=w({\delta }_{1}y+{\delta }_{2}z-{\delta }_{3}w-s)<w({\delta }_{1}y+{\delta }_{2}z)<Mw$, $M>0$ hence $w\left(t\right)\le w_0{e}^{Mt}$ on $[0,\omega )$. Therefore, $w\left(t\right)\le w_0{e}^{M\omega }$ on $[0,\omega )$, implying that $w\left(t\right)$ is bounded on $[0,\omega )$. Thus, $X\left(t\right)$ is bounded on $[0,\omega )$, which contradicts the standard ODE theory. Therefore, $\omega =+\infty$.    □

Proposition 6.

It cannot happen that ${lim}_{t\to +\infty }w\left(t\right)=+\infty$.

Proof. 

Assume by contradiction that $w\left(t\right)\to +\infty$. Then we have ${\displaystyle \frac{w}{{\delta }_3{w}^2+1}\to 0}$, and given that

$$z^{\prime }=zw\left(-{\delta }_2+k{\displaystyle \frac{{\delta }_{3}w}{{\delta }_3{w}^2+1}}\right),$$

we obtain that there exists $T>0$ such that $z^{\prime }\left(t\right)<0$ for all $t>T$. Therefore, $z\left(t\right)$ is decreasing on $(T,+\infty )$ and has a limit ${lim}_{t\to +\infty }z\left(t\right)=\ell \ge 0$. If $\ell >0$, then $z^{\prime }\left(t\right)\to -\infty$, which contradicts well-known theorems. Hence $\ell =0$. Now, consider $w^{\prime }\left(t\right)=w({\delta }_{1}y+{\delta }_{2}z-{\delta }_{3}w-s)$. Since $y\left(t\right)$ is bounded, $z\left(t\right)\to 0$, and $w\left(t\right)\to +\infty$, there exists $T_1>0$ such that $w^{\prime }\left(t\right)<0$ for all $t\ge T_1$. This implies $w\left(t\right)\le w\left(T_1\right)$ for all $t\ge T_1$, which contradicts the assumption $w\left(t\right)\to +\infty$. The thesis is thus proved.    □

Proposition 7.

It cannot happen that ${lim}_{t\to +\infty }z\left(t\right)=+\infty$.

Proof. 

We will prove that $z\left(t\right)\to +\infty$ implies $w\left(t\right)\to +\infty$, and by Proposition 6 this is a contradiction, hence the thesis. So, assume $z\left(t\right)\to +\infty$ as $t\to +\infty$. Fix $M>0$ and $t_M>0$ such that ${\displaystyle \frac{{\delta }_1}{{\delta }_3}z\left(t\right)>M+\frac{s}{{\delta }_3}}$ for all $t\ge t_M$.

We consider the following cases:

• Assume ${\displaystyle w\left(t\right)\le \frac{{\delta }_1}{{\delta }_3}y\left(t\right)+\frac{{\delta }_2}{{\delta }_3}z\left(t\right)-\frac{s}{{\delta }_3}}$ for all $t\ge t_M$. This implies $w^{\prime }\left(t\right)\ge 0$ on $[t_M,+\infty )$, hence there exists a limit $\ell ={lim}_{t\to +\infty }w\left(t\right)$. If $\ell <+\infty$, recalling that $\ell \ge w\left(t_M\right)>0$, from $w^{\prime }\left(t\right)=w({\delta }_{1}y+{\delta }_{2}z-{\delta }_{3}w-s)$ we would obtain $w^{\prime }\left(t\right)\to +\infty$, a contradiction with well-known theorems. Hence, it must be $\ell =+\infty$, that is, $w\left(t\right)\to +\infty$, but this contradicts Proposition 6. Thus, we can rule out this case, meaning that it is not possible that

  $$w\left(t\right)\le {\displaystyle \frac{{\delta }_1}{{\delta }_3}}y\left(t\right)+{\displaystyle \frac{{\delta }_2}{{\delta }_3}}z\left(t\right)-{\displaystyle \frac{s}{{\delta }_3}},\forall t\ge t_M.$$
• Therefore, there exists $s_M\ge t_M$ such that

  $$w\left(s_M\right)>{\displaystyle \frac{{\delta }_1}{{\delta }_3}}y\left(s_M\right)+{\displaystyle \frac{{\delta }_2}{{\delta }_3}}z\left(s_M\right)-{\displaystyle \frac{s}{{\delta }_3}}.$$
  Now, let us pick $t>s_M$. If

  $$w\left(t\right)\ge {\displaystyle \frac{{\delta }_1}{{\delta }_3}}y\left(t\right)+{\displaystyle \frac{{\delta }_2}{{\delta }_3}}z\left(t\right)-{\displaystyle \frac{s}{{\delta }_3}},$$

  we obtain $w\left(t\right)>M$ because $t>s_M\ge t_M$ hence ${\displaystyle \frac{{\delta }_1}{{\delta }_3}z\left(t\right)>M+\frac{s}{{\delta }_3}}$. On the other hand, assume

  $$w\left(t\right)<{\displaystyle \frac{{\delta }_1}{{\delta }_3}}y\left(t\right)+{\displaystyle \frac{{\delta }_2}{{\delta }_3}}z\left(t\right)-{\displaystyle \frac{s}{{\delta }_3}}.$$
  Define

  $${\tau }_M=inf\left\{\tau \in (s_M,t)|w\left(\sigma \right)<{\displaystyle \frac{{\delta }_1}{{\delta }_3}}y\left(\sigma \right)+{\displaystyle \frac{{\delta }_2}{{\delta }_3}}z\left(\sigma \right)-{\displaystyle \frac{s}{{\delta }_3}},\forall \sigma \in (\tau ,t)\right\}.$$
  By continuity, we get ${\tau }_M\in (s_M,t)$, ${\displaystyle w\left(\sigma \right)<\frac{{\delta }_1}{{\delta }_3}y\left(\sigma \right)+\frac{{\delta }_2}{{\delta }_3}z\left(\sigma \right)-\frac{s}{{\delta }_3}}$ for all $\sigma \in ({\tau }_M,t)$, and ${\displaystyle w\left({\tau }_M\right)=\frac{{\delta }_1}{{\delta }_3}y\left({\tau }_M\right)+\frac{{\delta }_2}{{\delta }_3}z\left({\tau }_M\right)-\frac{s}{{\delta }_3}}$. Hence, $w^{\prime }\left(t\right)>0$ on $({\tau }_M,t)$, so $w\left(t\right)>w\left({\tau }_M\right)>M$. We have thus proved that for all $t>s_M$, it holds that $w\left(t\right)>M$. Therefore, for any $M>0$, there exists $s_M>0$ such that $w\left(t\right)>M$ for all $t>s_M$, which means

  $$\underset{t\to +\infty }{lim}w\left(t\right)=+\infty .$$
  This is ruled out by Proposition 6, so we have obtained a contradiction, and the thesis is proved.

   □

3. Equilibrium Points

Equilibrium points are the solutions of the following system:

$$\begin{array}{cc}\hfill & y(\gamma \lambda -\gamma y-{\delta }_{1}w)=0\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill & zw\left(-{\delta }_2+k{\displaystyle \frac{{\delta }_{3}w}{{\delta }_3{w}^2+1}}\right)=0\hfill \end{array}$$

(5b)

$$\begin{array}{cc}\hfill & w\left({\delta }_{1}y+{\delta }_{2}z-{\delta }_{3}w-s\right)=0\hfill \end{array}$$

(5c)

Let us consider the different possible cases:

1)

$y=0,w=0$, no conditions on z. We obtain a family of equilibria $(0,Z,0)$ with $Z\ge 0$.

2)

$y=0$, $w\ne 0$, $z=0$. In this case, we obtain the equilibrium $\left({\displaystyle 0,0,-\frac{s}{{\delta }_3}}\right)$, which is not in $\overline{\mathcal{C}}$, so it is not of interest.

3)

$y=0$, $z\ne 0$, $w\ne 0$. In this case, from (), we obtain the equation for w:

$${\delta }_2{\delta }_3{w}^2-k{\delta }_{3}w+{\delta }_2=0,$$

which has the solutions

$$w_{1,2}={\displaystyle \frac{k{\delta }_3\pm \sqrt{k^2{\delta }_3^2-4{\delta }_2^2{\delta }_3}}{2{\delta }_2{\delta }_3}}={\displaystyle \frac{k\pm \sqrt{k^2-4\frac{{\delta }_2^2}{{\delta }_3}}}{2{\delta }_2}},$$

(6)

and we assume of course $k\ge {\displaystyle \frac{2{\delta }_2}{\sqrt{{\delta }_3}}}$ For z we obtain, from equation (),

$$z={\displaystyle \frac{{\delta }_3}{{\delta }_2}}w+{\displaystyle \frac{s}{{\delta }_2}}.$$

Hence, we have two equilibria:

$$P_i=\left(0,{\displaystyle \frac{{\delta }_3}{{\delta }_2}}{w}_i+{\displaystyle \frac{s}{{\delta }_2}},w_i\right),i=1,2.$$

(7)

4)

$y\ne 0$, $w=0$. We obtain $y=\lambda$ and no conditions on z, so we have a family of equilibria:

$$(\lambda ,Z,0),Z\ge 0.$$

(8)

5)

$y\ne 0$, $w\ne 0$, $z=0$. In this case, we obtain the linear system:

$$\left\{\begin{array}{cc}\hfill \gamma y+{\delta }_{1}w& =\gamma \lambda ,\hfill \\ \hfill {\delta }_{1}y-{\delta }_{3}w& =s,\hfill \end{array}\right.$$

whose solutions are

$$y_0={\displaystyle \frac{\gamma \lambda {\delta }_3+{\delta }_{1}s}{\gamma {\delta }_3+{\delta }_1^2}},w_0=\gamma {\displaystyle \frac{\lambda {\delta }_1-s}{\gamma {\delta }_3+{\delta }_1^2}}.$$

(9)

Thus, we obtain an equilibrium:

$$P_0=(y_0,0,w_0),{\delta }_1\lambda \ge s.$$

(10)

6)

$y\ne 0$, $z\ne 0$, $w\ne 0$. As before, we obtain the two values $w_i,i=1,2$, given in (6) for the variable w, while the values for y and z are now given by

$$y_i=\lambda -{\displaystyle \frac{{\delta }_1}{\gamma }}{w}_i,z_i=-{\displaystyle \frac{{\delta }_1}{{\delta }_2}}{y}_i+{\displaystyle \frac{{\delta }_3}{{\delta }_2}}{w}_i+{\displaystyle \frac{s}{{\delta }_2}},i=1,2,$$

(11)

giving rise to two steady states:

$$P_i=(y_i,z_i,w_i),i=1,2.$$

(12)

4. Stability of Equilibria: Cases 1 and 3

The Jacobian matrix of the vector field $F\left(X\right)$ at a point $X=(y,z,w)$ is given by:

$$JF\left(X\right)=\left(\begin{array}{ccc}\gamma \lambda -2\gamma y-{\delta }_{1}w& 0& -{\delta }_{1}y\\ \\ 0& {\displaystyle -{\delta }_{2}w+k\frac{{\delta }_3{w}^2}{{\delta }_3{w}^2+1}}& {\displaystyle -{\delta }_{2}z+2k{\delta }_3\frac{zw}{{({\delta }_3{w}^2+1)}^2}}\\ \\ {\delta }_{1}w& {\delta }_{2}w& {\delta }_{1}y+{\delta }_{2}z-2{\delta }_{3}w-s\end{array}\right).$$

(13)

Let us now analyze the stability of the equilibria previously identified. We will not consider case 2) because, as we have seen above, the equilibrium is not in $\overline{\mathcal{C}}$.

4.1. Case 1

In case 1), we have $y=0$, so the first row of the Jacobian matrix is the vector $(\gamma \lambda ,0,0)$. This means that ${\rho }_1=\gamma \lambda >0$ is an eigenvalue, hence all these equilibria are unstable.

4.2. Case 3

In case 3), we also have $y=w=0$, so the first row of the Jacobian matrix is the vector $(\gamma \lambda -{\delta }_1{w}_i,0,0)$. This means that ${\rho }_1=\gamma \lambda -{\delta }_1{w}_i$ is an eigenvalue. We will obtain asymptotic results as ${\delta }_i\to 0^{+}$ for $i=1,2,3$. In this section we will assume, for simplicity, linear relations between the coefficients ${\delta }_i$. Specifically, we assume ${\delta }_3=\delta$, ${\delta }_1=a\delta$, and ${\delta }_2=b\delta$ with $a,b>0$, and we will study the asymptotic behavior of equilibria as $\delta \to 0^{+}$. Now we want to obtain the asymptotic developments of w. Recall that

$$\begin{array}{cc}\hfill w_1& ={\displaystyle \frac{k-\sqrt{k^2-4\frac{{\delta }_2^2}{{\delta }_3}}}{2{\delta }_2}}={\displaystyle \frac{k-\sqrt{k^2-4b^2\delta }}{2b\delta }},\hfill \\ \hfill w_2& ={\displaystyle \frac{k+\sqrt{k^2-4b^2\delta }}{2b\delta }}.\hfill \end{array}$$

Let us now find what happens as $\delta \to 0^{+}$. First, we have

$$\sqrt{k^2-4b^2\delta }=k-{\displaystyle \frac{2b^2}{k}}\delta -{\displaystyle \frac{2b^4}{k^3}}{\delta }^2+\mathcal{O}\left({\delta }^3\right),$$

hence

$$\begin{array}{cc}\hfill w_1& ={\displaystyle \frac{b}{k}}+{\displaystyle \frac{b^3}{k^3}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ \hfill w_2& ={\displaystyle \frac{k}{b}}{\displaystyle \frac{1}{\delta }}-{\displaystyle \frac{b}{k}}-{\displaystyle \frac{b^3}{k^3}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \end{array}$$

For the equilibrium $P_1=(0,z_1,w_1)$ we obtain the eigenvalue ${\rho }_1=\gamma \lambda +\mathcal{O}\left(\delta \right)>0$, indicating that this is an unstable equilibrium, as $\delta \to 0^{+}$. For the equilibrium $P_2$ the eigenvalue becomes

$${\rho }_1=\gamma \lambda -{\displaystyle \frac{a}{b}}k-{\displaystyle \frac{ab}{k}}\delta +\mathcal{O}\left({\delta }^2\right).$$

Therefore, if ${\displaystyle \gamma \lambda -\frac{a}{b}k>0}$, the equilibrium $P_2$ is also unstable for $\delta \to 0^{+}$. Let us now analyze the case ${\displaystyle \gamma \lambda -\frac{a}{b}k<0}$. We have

$$k{\displaystyle \frac{{\delta }_{3}w}{{\delta }_3{w}^2+1}}={\delta }_2,$$

so the entry $J_{23}$ of the Jacobian matrix (13), computed at the points (12), becomes

$$\begin{array}{cc}& {\displaystyle -{\delta }_{2}z+2k{\delta }_3\frac{zw}{{({\delta }_3{w}^2+1)}^2}=-{\delta }_{2}z+2{\delta }_2\frac{z}{{\delta }_3{w}^2+1}}\hfill \\ & {\displaystyle =-{\delta }_{2}z+2{\delta }_2^2\frac{z}{k{\delta }_{3}w}={\delta }_{2}z\left(2\frac{{\delta }_2}{k{\delta }_{3}w}-1\right)=\delta bz\left(2\frac{b}{kw}-1\right).}\hfill \end{array}$$

Hence, the Jacobian matrix in this case becomes

$$JF\left(P_2\right)=\left(\begin{array}{ccc}\gamma \lambda -a\delta w_2& 0& 0\\ \\ 0& 0& {\displaystyle \delta b{z}_2\left(2\frac{b}{k{w}_2}-1\right)}\\ \\ a\delta w_2& b\delta w_2& -\delta w_2\end{array}\right).$$

(14)

We have

$$z_2={\displaystyle \frac{{\delta }_3}{{\delta }_2}}{w}_2+{\displaystyle \frac{s}{{\delta }_2}}={\displaystyle \frac{1}{b}}{w}_2+{\displaystyle \frac{s}{b\delta }}=\left({\displaystyle \frac{k}{b^2}}+{\displaystyle \frac{s}{b}}\right){\displaystyle \frac{1}{\delta }}-{\displaystyle \frac{1}{k}}+\mathcal{O}\left(\delta \right).$$

We obtain the following developments:

$$\begin{array}{cc}& \delta z_2=\left({\displaystyle \frac{k}{b^2}}+{\displaystyle \frac{s}{b}}\right)-{\displaystyle \frac{1}{k}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ & {\displaystyle \frac{1}{w_2}}={\displaystyle \frac{b}{k}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ & 2{\displaystyle \frac{b}{k{w}_2}}-1=-1+{\displaystyle \frac{2b^2}{k^2}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ & \delta b{z}_2\left(2{\displaystyle \frac{b}{k{w}_2}}-1\right)=-s-{\displaystyle \frac{k}{b}}+\mathcal{O}\left(\delta \right).\hfill \end{array}$$

As a consequence, we have

$$\begin{array}{cc}\hfill J_{23}\left(P_2\right)& =\delta b{z}_2\left(2{\displaystyle \frac{b}{k{w}_2}}-1\right)=-s-{\displaystyle \frac{k}{b}}+\mathcal{O}\left(\delta \right),\hfill \\ \hfill J_{31}\left(P_2\right)& =a\delta w_2={\displaystyle \frac{ak}{b}}+\mathcal{O}\left(\delta \right),\hfill \\ \hfill J_{32}\left(P_2\right)& =b\delta w_2=k+\mathcal{O}\left(\delta \right),\hfill \\ \hfill J_{33}\left(P_2\right)& =-\delta w_2=-{\displaystyle \frac{k}{b}}+\mathcal{O}\left(\delta \right).\hfill \end{array}$$

Given these asymptotic expansions, we obtain

$$JF\left(P_2\right)=\left(\begin{array}{ccc}{\displaystyle \gamma \lambda -\frac{a}{b}k+\mathcal{O}\left(\delta \right)}& 0& 0\\ \\ 0& 0& {\displaystyle -s-\frac{k}{b}+\mathcal{O}\left(\delta \right)}\\ \\ {\displaystyle \frac{ak}{b}+\mathcal{O}\left(\delta \right)}& {\displaystyle k+\mathcal{O}\left(\delta \right)}& {\displaystyle -\frac{k}{b}+\mathcal{O}\left(\delta \right)}\end{array}\right).$$

(15)

This matrix has an eigenvalue ${\rho }_1=\gamma \lambda -{\displaystyle \frac{a}{b}}k-{\displaystyle \frac{ab}{k}}\delta +\mathcal{O}\left({\delta }^2\right)$ and we are assuming $\gamma \lambda -{\displaystyle \frac{a}{b}}k<0$, hence ${\rho }_1<0$ as $\delta \to 0^{+}$. The other eigenvalues are obtained from the submatrix

$$SJ=\left(\begin{array}{cc}{\displaystyle 0}& {\displaystyle -s-\frac{k}{b}+\mathcal{O}\left(\delta \right)}\\ \\ {\displaystyle k+\mathcal{O}\left(\delta \right)}& {\displaystyle -\frac{k}{b}+\mathcal{O}\left(\delta \right)}\end{array}\right)$$

(16)

We have Trace($SJ$)=$-{\displaystyle \frac{k}{b}}+\mathcal{O}\left(\delta \right)<0$ and Det($SJ$)=$k(s+{\displaystyle \frac{k}{b}})+\mathcal{O}\left(\delta \right)>0$, so the eigenvalues of $SJ$ have negative real part as $\delta \to 0^{+}$. Hence all the eigenvalues of $JF\left(P_2\right)$ have negative real part, and $P_2$ is an asymptotically stable equilibrium. The following proposition summarizes the results of this section

Proposition 8.

All the equilibria of Case 1 are unstable (for any value of the parameters). As to Case 3, the equilibrium $P_1$ is unstable as $\delta \to 0^{+}$, while equilibrium $P_2$ is unstable if ${\displaystyle \gamma \lambda -\frac{a}{b}k>0}$ and is asymptotically stable if ${\displaystyle \gamma \lambda -\frac{a}{b}k<0}$ (and $\delta \to 0^{+}$).

5. Stability of Equilibria: Case 4

In case 4 the equilibria are given by $P=(\lambda ,Z,0)$ with $Z\ge 0$, and the Jacobian matrix is

$$JF\left(P\right)=\left(\begin{array}{ccc}-\gamma \lambda & 0& -{\delta }_1\lambda \\ \\ 0& 0& -{\delta }_{2}Z\\ \\ 0& 0& {\delta }_1\lambda +{\delta }_{2}Z-s\end{array}\right).$$

(17)

The eigenvalues are ${\rho }_1=-\gamma \lambda <0$, ${\rho }_2=0$, and ${\rho }_3={\delta }_1\lambda +{\delta }_{2}Z-s$. If ${\delta }_1\lambda +{\delta }_{2}Z-s>0$, then ${\rho }_3>0$ and the equilibria are unstable. Let us study the case ${\delta }_1\lambda +{\delta }_{2}Z-s<0$, that is, ${\displaystyle Z<\frac{s}{{\delta }_2}-\frac{{\delta }_1}{{\delta }_2}\lambda }$. In this case the study of the Jacobian matrix $JF\left(P\right)$ is not useful, because two eigenvalues are negative and one is zero. We will however obtain a stability result with an ad hoc argument. We will deal only with trajectories starting in $\mathcal{C}$, because in our model these are the relevant trajectories. Let us start by introducing the rectangular neighborhoods $I_{ϵ}$ of P as follows, for $ϵ>0$:

$$I_{ϵ}=(\lambda -ϵ,\lambda +ϵ)\times (Z-ϵ,Z+ϵ)\times (-ϵ,ϵ).$$

After this we introduce the constants ${ϵ}_1,{ϵ}_2,\varrho$ as follows: since ${\displaystyle Z<\frac{s}{{\delta }_2}-\frac{{\delta }_1}{{\delta }_2}\lambda }$, it is possible to choose two positive numbers $\varrho ,{ϵ}_1>0$ such that for all $ϵ\in (0,{ϵ}_1)$, the following holds:

$${\delta }_1(\lambda +ϵ)+{\delta }_2(Z+ϵ)-s<-\varrho .$$

Also we pick ${ϵ}_2>0$ such that ${\displaystyle -{\delta }_2+k{\delta }_3{ϵ}_2<-\frac{{\delta }_2}{2}}$, that is, ${\displaystyle {ϵ}_2<\frac{{\delta }_2}{2k{\delta }_3}}$.

We then state the result we are going to prove:

Proposition 9.

There is $\overline{ϵ}>0$ such that for all $0<ϵ<\overline{ϵ}$ there is $\eta >0$ such that, if we pick $Q=(y_0,z_0,w_0)\in I_{\eta }\cap \mathcal{C}$ and we call $X(t,Q)=\left(y\right(t),z(t),w(t\left)\right)$ the trajectory starting from Q at $t=0$, then $X(t,Q)\in I_{ϵ}$ for all $t\ge 0$. Also, it holds:

$i)$

$z,w$ are decreasing function, and $0<w\left(t\right)<w_0{e}^{-\varrho t}$ for all $t\ge 0$.

$ii)$

${lim}_{t\to +\infty }y\left(t\right)=\lambda$.

$iii)$

There exists a positive constant $c_1>0$, independent of ϵ and $y_0,z_0,w_0$, such that the following holds:

$$z_0{e}^{-c_1{w}_0}\le z\left(t\right)\le z_0,\forall t\ge 0.$$

Notice that Proposition 9 states that P is Lyapunov-stable and give some information on the asymptotic behavior of the variables as $t\to \infty$, but does not imply asymptotic stability. Indeed, if $z_0<Z$, $z\left(t\right)$ does not converge to Z.

The proof of Proposition 9 will be done in several steps.

Step 1: For any $ϵ>0$, if $y_0<\lambda +ϵ$, then $y\left(t\right)<\lambda +ϵ$ for all $t\ge 0$.

Proof. 

By the arguments of Proposition 4 it is to get that $y\left(t\right)\le max\{y_0,\lambda \}$, and if $y_0<\lambda +ϵ$, then obviously $max\{y_0,\lambda \}<\lambda +ϵ$.    □

Step 2: Assume $ϵ<min\{{ϵ}_1,{ϵ}_2\}$ and $y_0<\lambda +ϵ$, $z_0<Z+ϵ$, $w_0<ϵ$. Then it holds $z^{\prime }\left(t\right)<0$ and $w^{\prime }\left(t\right)<0$ for all $t\ge 0$.

Proof. 

As $ϵ<{ϵ}_2$ it is easy to get

$$-{\delta }_2+k{\displaystyle \frac{{\delta }_3{w}_0}{{\delta }_3{w}_0^2+1}}<-{\delta }_2+k{\delta }_3{w}_0<-{\displaystyle \frac{{\delta }_2}{2}}.$$

Hence we get

$$\begin{array}{cc}\hfill z^{\prime }\left(0\right)& =z_0{w}_0\left(-{\delta }_2+k{\displaystyle \frac{{\delta }_3{w}_0}{{\delta }_3{w}_0^2+1}}\right)<-{\displaystyle \frac{{\delta }_2}{2}}{z}_0{w}_0<0,\hfill \\ \hfill w^{\prime }\left(0\right)& =w_0\left({\delta }_1{y}_0+{\delta }_2{z}_0-{\delta }_3{w}_0-s\right)<w_0\left({\delta }_1(\lambda +ϵ)+{\delta }_2(Z+ϵ)-s\right)<-\varrho w_0<0.\hfill \end{array}$$

As a consequence, there exists $t_1>0$ such that for all $t\in [0,t_1]$, it holds $z^{\prime }\left(t\right)<0$ and $w^{\prime }\left(t\right)<0$.

Now we define

$$\begin{array}{cc}\hfill {\tau }_1& =sup\left\{t>0|z^{\prime }\left(s\right)<0,\forall s\in [0,t]\right\},\hfill \\ \hfill {\tau }_2& =sup\left\{t>0|w^{\prime }\left(s\right)<0,\forall s\in [0,t]\right\}.\hfill \end{array}$$

We know that ${\tau }_i>t_1$. We want to prove that ${\tau }_1={\tau }_2=+\infty$, which gives the thesis. Let us consider all other possible cases and show that they lead to a contradiction.

i)

If ${\tau }_1<{\tau }_2$ then ${\tau }_1<+\infty$; we have $w^{\prime }\left(t\right)<0$ for all $t\in [0,{\tau }_2)$, hence $w\left(t\right)<w_0$. Therefore,

$$-{\delta }_2+k{\displaystyle \frac{{\delta }_{3}w\left(t\right)}{{\delta }_3{w}^2\left(t\right)+1}}<-{\delta }_2+k{\delta }_{3}w\left(t\right)<-{\delta }_2+k{\delta }_3{w}_0<-{\displaystyle \frac{{\delta }_2}{2}},\forall t\in [0,{\tau }_2),$$

implying ${\displaystyle z^{\prime }\left(t\right)<-\frac{{\delta }_2}{2}z\left(t\right)w\left(t\right)<0}$ for $t\in [0,{\tau }_2)$. Since ${\tau }_1<{\tau }_2$, we obtain a contradiction with the definition of ${\tau }_1$.

ii)

If ${\tau }_2<{\tau }_1$ we have $z^{\prime }\left(t\right)<0$ for all $t\in [0,{\tau }_1)$, hence $z\left(t\right)<z_0$. Therefore,

$${\delta }_{1}y\left(t\right)+{\delta }_{2}z\left(t\right)-s<{\delta }_1(\lambda +ϵ)+{\delta }_2{z}_0-s<{\delta }_1(\lambda +ϵ)+{\delta }_2(Z+ϵ)-s<-\varrho ,$$

hence $w^{\prime }\left(t\right)<-\varrho w\left(t\right)<0$ for all $t\in [0,{\tau }_1)$. Since ${\tau }_2<{\tau }_1$, this is a contradiction with the definition of ${\tau }_2$.

iii)

If ${\tau }_1={\tau }_2<+\infty$, by continuity, we easily obtain $z^{\prime }\left(t\right)<0$ and $w^{\prime }\left(t\right)<0$ on $[0,{\tau }_1)$, with $z^{\prime }\left({\tau }_1\right)=w^{\prime }\left({\tau }_1\right)=0$. Therefore, as above, $z^{\prime }\left(t\right)<z\left(0\right)<Z+ϵ$ for all $t\in [0,{\tau }_1)$, and $w^{\prime }\left(t\right)<-\varrho w\left(t\right)$ on $[0,{\tau }_1)$, hence $w^{\prime }\left({\tau }_1\right)\le -\varrho w\left({\tau }_1\right)<0$, which is a contradiction.

Thus, the only possibility left is ${\tau }_1={\tau }_2=+\infty$, which proves the thesis.

   □

Step 3: Hypotheses as in Step 2. Then, it is $0<w\left(t\right)<w_0{e}^{-\varrho t}$ for all $t\ge 0$.

Proof. 

By the previous steps we know that $y\left(t\right)<\lambda +ϵ$ and $z\left(t\right)<Z+ϵ$ for all $t\ge 0$. Using the equation for $w^{\prime }$ we obtain, as in Step 2:

$$w^{\prime }=w\left({\delta }_{1}y+{\delta }_{2}z-{\delta }_{3}w-s\right)<w\left({\delta }_1(\lambda +ϵ)+{\delta }_2(Z+ϵ)-s\right)<-\varrho w,$$

and a simple integration gives the result.

   □

Step 4: Hypotheses as in Step 2. Then there exist a positive constants $c_1>0$, independent of $ϵ$ and $y_0,z_0,w_0$, such that the following holds:

$$z_0{e}^{-c_1{w}_0}\le z\left(t\right)\le z_0,\forall t\ge 0.$$

Proof. 

The inequality $z\left(t\right)\le z_0$ is obvious from Step 2. Now, using Step 3 and the equation for $z^{\prime }$ in the model (1), we have

$$z^{\prime }\left(t\right)>-{\delta }_{2}z\left(t\right)w\left(t\right)\ge -{\delta }_2{w}_0{e}^{-\varrho t}z\left(t\right),\forall t\ge 0,$$

hence ${\displaystyle \frac{z^{\prime }\left(t\right)}{z\left(t\right)}>-{\delta }_2{w}_0{e}^{-\varrho t}}$, so that

$$log{\displaystyle \frac{z\left(t\right)}{z_0}}\ge -{\delta }_2{w}_0\left[-{\displaystyle \frac{1}{\varrho }}(e^{-\varrho t}-1)\right]={\displaystyle \frac{{\delta }_2}{\varrho }}{w}_0(e^{-\varrho t}-1)\ge -{\displaystyle \frac{{\delta }_2}{\varrho }}{w}_0,$$

so that $z\left(t\right)\ge z_0{e}^{{\displaystyle -\frac{{\delta }_2}{\varrho }{w}_0}}=z_0{e}^{-c_1{w}_0}$ with ${\displaystyle c_1=\frac{{\delta }_2}{\varrho }}$.

   □

Step 5: There are ${ϵ}_3<min\{{ϵ}_1,{ϵ}_2\}$ and $0<{\mu }_1<1$ sucht that, if $0<ϵ<{ϵ}_3$, $0<\mu <{\mu }_1$, $Z-\mu ϵ<z_0<Z+\mu ϵ$, $0<w_0<\mu ϵ$ and $y_0<\lambda +ϵ$, then $Z-ϵ<z\left(t\right)<Z+ϵ$ for all $t\ge 0$.

Proof. 

Notice that the thesis is obvious if $Z=0$, because $z\left(t\right)$ is decreasing and positive. Hence we assume $Z>0$. Let us take $ϵ<min\{{ϵ}_1,{ϵ}_2\}$. We have $z\left(t\right)<Z+ϵ$ for all $t\ge 0$, by the hypotheses and the previous Steps. Now, assume $0<ϵ<Z$ and define

$${\mu }_1={\displaystyle \frac{1}{1+c_{1}Z}}<1.$$

Let us pick $0<\mu <{\mu }_1$ and consider the inequality

$$\mu ϵ<{\displaystyle \frac{1}{c_1}}log\left({\displaystyle \frac{Z-\mu ϵ}{Z-ϵ}}\right).$$

(18)

By the standard asymptotic development of $log(1+x)$ this is equivalent to

$$\mu ϵ<{\displaystyle \frac{1}{c_1}}{\displaystyle \frac{(1-\mu )ϵ}{Z-ϵ}}+\mathcal{O}\left({ϵ}^2\right),$$

and by simple computations this is equivalent to

$$\mu <{\displaystyle \frac{1}{1+c_{1}Z}}+\mathcal{O}\left({ϵ}^2\right).$$

By our choice of $\mu <{\mu }_1$, we get that there is ${ϵ}_3$ such that for all $0<ϵ<{ϵ}_3$ the inequality (18) holds, and of course we can choose ${ϵ}_3<min\{{ϵ}_1,{ϵ}_2,Z\}$. Now from (18) we deduce, as $w_0<\mu ϵ$ and $z_0>Z-ϵ$,

$$w_0<{\displaystyle \frac{1}{c_1}}log\left({\displaystyle \frac{Z-\mu ϵ}{Z-ϵ}}\right)\mathrm{hence}-c_1{w}_0>log\left({\displaystyle \frac{Z-ϵ}{Z-\mu ϵ}}\right),$$

so that

$$e^{-c_1{w}_0}>{\displaystyle \frac{Z-ϵ}{Z-\mu ϵ}}>{\displaystyle \frac{Z-ϵ}{z_0}}\mathrm{hence}{z}_0{e}^{-c_1{w}_0}>Z-ϵ,$$

and we get the result thanks to Step 4.

   □

For the next Step 6, we introduce ${ϵ}_4=\gamma \lambda /{\delta }_1>0$, and we notice that for all $0<ϵ<{ϵ}_4$ and all $\mu \in [0,1]$ it is $\lambda -{\delta }_1\mu ϵ/\gamma >0$.

Step 6: Pick $ϵ<min\{{ϵ}_1,{ϵ}_2,{ϵ}_3,{ϵ}_4\}$ and $\mu <min\{{\mu }_1,\gamma /{\delta }_1\}$. Assume $0<w_0<\mu ϵ$, $z_0<Z+\mu ϵ$, $\lambda -ϵ<y_0<\lambda +ϵ$. Then $y\left(t\right)>\lambda -ϵ$ for all $t\ge 0$.

Proof. 

By the equation for $y^{\prime }$ and the previous Steps we get

$$y^{\prime }=\gamma y\left(\lambda -y\right)-{\delta }_{1}yw>\gamma y\left(\lambda -y\right)-{\delta }_{1}y\mu ϵ=\gamma y\left(\lambda -y-{\delta }_1\mu ϵ\right)=\gamma y\left({\lambda }_1-y\right)$$

where ${\lambda }_1=\lambda -{\delta }_1\mu ϵ/\gamma >0$. Let us consider the function $\overline{y}$ which is the solution of the following logistic Cauchy problem:

$$\left\{\begin{array}{cc}\hfill {\overline{y}}^{\prime }\left(t\right)& =\gamma \overline{y}({\lambda }_1-\overline{y}),\hfill \\ \hfill \overline{y}\left(0\right)& =y_0\hfill \end{array}\right.$$

(19)

By standard comparison principles we obtain $y\left(t\right)>\overline{y}\left(t\right)$ for all $t>0$. Now, the well-known properties of logistic equation give that, if $y_0<{\lambda }_1$ then $\overline{y}\left(t\right)>y_0$ for all $t>0$, hence $y\left(t\right)\ge y_0>\lambda -ϵ$ for all $t\ge 0$. On the other hand, if $y_0\ge {\lambda }_1$, we know that $\overline{y}\left(t\right)\ge {\lambda }_1$ for all t, so also $y\left(t\right)>{\lambda }_1=\lambda -{\delta }_1\mu ϵ/\gamma >\lambda -ϵ$, by our choice of $\mu$. The proof of Step 6 is now complete.

   □

Step 7: Assume the hypotheses of Step 6. Then it is ${lim}_{t\to +\infty }y\left(t\right)=\lambda$.

Proof. 

We will study ${lim\; sup}_{t\to +\infty }y\left(t\right)$ and ${lim\; inf}_{t\to +\infty }y\left(t\right)$, and we will prove that they are equal. Let us start by proving that ${lim\; sup}_{t\to +\infty }y\left(t\right)\le \lambda$. We use the arguments of Proposition 4, and we distinguish two cases, as we did there:

$i)$ $y\left(t\right)\ge \lambda$ for all $t\ge 0$: in this case we know that $y^{\prime }\left(t\right)<0$ for all t, so there exists $L={lim}_{t\to +\infty }y\left(t\right)$ and it holds $L\ge \lambda$. If $L>\lambda$ then, recalling Step 3,

$$\underset{t\to +\infty }{lim}{y}^{\prime }\left(t\right)=\underset{t\to +\infty }{lim}\gamma y\left(t\right)(\lambda -y\left(t\right))-{\delta }_{1}y\left(t\right)w\left(t\right)=\gamma L(\lambda -L)<0$$

which is a contradiction with standard result. Hence $L=\lambda$, which trivially implies ${lim\; sup}_{t\to +\infty }y\left(t\right)=\lambda$.

$ii)$ $y\left(t_0\right)<\lambda$ for some $t_0\ge 0$: in this case we have proved, in Proposition 4, that $y\left(t\right)<\lambda$ for all $t\ge t_0$, which of course implies ${lim\; sup}_{t\to +\infty }y\left(t\right)\le \lambda$. So we get, in any case, ${lim\; sup}_{t\to +\infty }y\left(t\right)\le \lambda$. Let us now prove that ${lim\; inf}_{t\to +\infty }y\left(t\right)\ge \lambda$. Fix any $0<\sigma <\gamma \lambda /{\delta }_1$, and pick $t_{\sigma }\ge 0$ such that for all $t\ge t_{\sigma }$ it holds $w\left(t\right)<\sigma$. From the equation for $y^{\prime }$ we get

$$y^{\prime }\left(t\right)=\gamma y\left(t\right)\left(\lambda -y\left(t\right)-{\displaystyle \frac{{\delta }_1}{\gamma }}w\left(t\right)\right)>\gamma y\left(t\right)\left(\lambda -y\left(t\right)-{\displaystyle \frac{{\delta }_1}{\gamma }}\sigma \right)=\gamma y\left(t\right)\left({\lambda }_1-y\left(t\right)\right),$$

where ${\lambda }_1=\lambda -{\displaystyle \frac{{\delta }_1}{\gamma }}\sigma >0$ and the inequality holds for all $t\ge t_{\sigma }$. Now we pick the value ${\overline{y}}_{\sigma }={\displaystyle \frac{1}{2}}min\{y\left(t_{\sigma }\right),{\lambda }_1\}>0$ and consider the function $\overline{y}\left(t\right)$ which is the solution of the following logistic Cauchy problem

$$\left\{\begin{array}{cc}\hfill {\overline{y}}^{\prime }\left(t\right)& =\gamma \overline{y}({\lambda }_1-\overline{y}),\hfill \\ \hfill \overline{y}\left(t_{\sigma }\right)& ={\overline{y}}_{\sigma }\hfill \end{array}\right.$$

(20)

As $0<\overline{y}\left(t_{\sigma }\right)<{\lambda }_1$, standard results give that $\overline{y}$ is increasing for $t\ge t_{\sigma }$ and ${lim}_{t\to +\infty }\overline{y}\left(t\right)={\lambda }_1$. On the other hand, an easy application of comparison principles gives that $y\left(t\right)>\overline{y}\left(t\right)$ for all $t\ge t_{\sigma }$, hence ${lim\; inf}_{t\to +\infty }y\left(t\right)\ge {\lambda }_1=\lambda -{\displaystyle \frac{{\delta }_1}{\gamma }}\sigma$. As this hold for all $\sigma <\gamma \lambda /{\delta }_1$, we get ${lim\; inf}_{t\to +\infty }y\left(t\right)\ge \lambda$.

We can now conclude the argument: we have obtained ${lim\; inf}_{t\to +\infty }y\left(t\right)\ge \lambda \ge {lim\; sup}_{t\to +\infty }y\left(t\right)$, and this obviously implies

$$\underset{t\to +\infty }{lim\; inf}y\left(t\right)=\underset{t\to +\infty }{lim\; sup}y\left(t\right)=\lambda =\underset{t\to +\infty }{lim}y\left(t\right).$$

The thesis is proved.

   □

Step 8: We now end the proof of Proposition 9: it is enough to choose $\overline{ϵ}=min\{{ϵ}_1,{ϵ}_2,{ϵ}_3,{ϵ}_4\}$ and $0<\mu <min\{{\mu }_1,\gamma /{\delta }_1\}$ and, for any $0<ϵ<\overline{ϵ}$, to define $\eta =\mu ϵ$. The previous Steps give all the result stated in Proposition 9.

Let us now collect the result of this section.

Proposition 10.

Let $P=(\lambda ,Z,0)$ with $Z\ge 0$. If ${\delta }_1\lambda +{\delta }_{2}Z-s>0$, then the equilibrium point P is unstable. If ${\delta }_1\lambda +{\delta }_{2}Z-s<0$, then P is Lyapunov-stable and the asymptotic result of Proposition 9 hold.

Remark 1.

We notice that the results of this section are not asymptotic results, but hold for any value of the parameters, provided the hypotheses are satisfied.

6. Stability of Equilibria: Cases 5 and 6

We now deal with cases 5 and 6. As in section 3 we will study the stability and instability properties of equilibria when ${\delta }_i\to 0^{+}$ and we will assume, for simplicity, linear relations between the coefficients ${\delta }_i$. So, as above, we set ${\delta }_3=\delta$, ${\delta }_1=a\delta$, and ${\delta }_2=b\delta$ with $a,b>0$, and we will study the asymptotic behavior of equilibria as $\delta \to 0^{+}$. This approach rules out case 5) because, as ${\delta }_1\to 0^{+}$, from (9) we get $w_0<0$, so the steady state is not in $\mathcal{C}$. We can repeat the same computations of Case 3, obtaining that the Jacobian matrix (13), computed at the points (12), is as follows:

$$JF\left(P_i\right)=\left(\begin{array}{ccc}-\gamma y_i& 0& -a\delta y_i\\ \\ 0& 0& {\displaystyle \delta b{z}_i\left(2\frac{b}{k{w}_i}-1\right)}\\ \\ a\delta w_i& b\delta w_i& -\delta w_i\end{array}\right),i=1,2.$$

(21)

Now we want to obtain the asymptotic developments of y, z, and w. For the values $w_i$ ($i=1,2$) the same computations of Case 3 hold, while for $y_i,z_i$ we use 11. Hence we have, for $P_1$,

$$\begin{array}{cc}\hfill w_1& ={\displaystyle \frac{b}{k}}+{\displaystyle \frac{b^3}{k^3}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ \hfill y_1& =\lambda -{\displaystyle \frac{ab}{\gamma k}}\delta -{\displaystyle \frac{a{b}^3}{\gamma k^3}}{\delta }^2+\mathcal{O}\left({\delta }^3\right),\hfill \\ \hfill z_1& ={\displaystyle \frac{s}{b\delta }}-{\displaystyle \frac{a\lambda }{b}}+{\displaystyle \frac{1}{k}}+\left({\displaystyle \frac{a^2}{k\gamma }}+{\displaystyle \frac{b^2}{k^3}}\right)\delta +\mathcal{O}\left({\delta }^2\right).\hfill \end{array}$$

Looking at the Jacobian matrix (21), we obtain in this case

$$\begin{array}{cc}\hfill J_{11}\left(P_1\right)& =-\gamma y_1=-\gamma \lambda +{\displaystyle \frac{ab}{k}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ \hfill J_{13}\left(P_1\right)& =-a\delta y_1=-a\lambda \delta +\mathcal{O}\left({\delta }^2\right).\hfill \end{array}$$

Considering

$$\begin{array}{cc}& \delta z_1={\displaystyle \frac{s}{b}}+\mathcal{O}\left(\delta \right),\hfill \\ & {\displaystyle \frac{1}{w_1}}={\displaystyle \frac{k}{b}}+\mathcal{O}\left(\delta \right),\hfill \\ & 2{\displaystyle \frac{b}{k{w}_1}}-1=1+\mathcal{O}\left(\delta \right),\hfill \\ & \delta b{z}_1\left(2{\displaystyle \frac{b}{k{w}_1}}-1\right)=s+\mathcal{O}\left(\delta \right),\hfill \end{array}$$

we have

$$\begin{array}{cc}\hfill J_{23}\left(P_1\right)& =\delta b{z}_1\left(2{\displaystyle \frac{b}{k{w}_1}}-1\right)=s+\mathcal{O}\left(\delta \right),\hfill \\ \hfill J_{31}\left(P_1\right)& =a\delta w_1={\displaystyle \frac{ab}{k}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ \hfill J_{32}\left(P_1\right)& =b\delta w_1={\displaystyle \frac{b^2}{k}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ \hfill J_{33}\left(P_1\right)& =-\delta w_1=-{\displaystyle \frac{b}{k}}\delta +\mathcal{O}\left({\delta }^2\right).\hfill \end{array}$$

Hence, the Jacobian matrix is given by

$$JF\left(P_1\right)=\left(\begin{array}{ccc}{\displaystyle -\gamma \lambda +\frac{ab}{k}\delta +\mathcal{O}\left({\delta }^2\right)}& 0& -a\lambda \delta +\mathcal{O}\left({\delta }^2\right)\\ \\ 0& 0& s+\mathcal{O}\left(\delta \right)\\ \\ {\displaystyle \frac{ab}{k}\delta +\mathcal{O}\left({\delta }^2\right)}& {\displaystyle \frac{b^2}{k}\delta +\mathcal{O}\left({\delta }^2\right)}& {\displaystyle -\frac{b}{k}\delta +\mathcal{O}\left({\delta }^2\right)}\end{array}\right).$$

(22)

Its determinant is given by Det$JF\left(P_1\right)={\displaystyle \frac{b^2\gamma \lambda s}{k}}\delta +\mathcal{O}\left({\delta }^2\right)$. Denoting by $p_3\left(\rho \right)$ = Det$\left(\rho I-JF\left(P_1\right)\right)$ its characteristic polynomial, we obtain $p_3\left(0\right)=$ Det$\left(-JF\left(P_1\right)\right)=-$Det$JF\left(P_1\right)=-{\displaystyle \frac{b^2\gamma \lambda s}{k}}\delta +\mathcal{O}\left({\delta }^2\right)<0$ as $\delta \to 0^{+}$. As $p_3\left(\rho \right)={\rho }^3+$ lower order terms, we obtain that $p_3\left(\rho \right)$ has a real positive root, which is a real positive eigenvalue of $JF\left(P_1\right)$. We then deduce that the equilibrium $P_1=(y_1,z_1,w_1)$ is unstable as $\delta \to 0^{+}$. Let us now see what happens for $P_2$. We have:

$$\begin{array}{cc}\hfill w_2& ={\displaystyle \frac{k}{b}}{\displaystyle \frac{1}{\delta }}-{\displaystyle \frac{b}{k}}-{\displaystyle \frac{b^3}{k^3}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ \hfill y_2& =\lambda -{\displaystyle \frac{ak}{b\gamma }}+{\displaystyle \frac{ab}{k\gamma }}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ \hfill z_2& =\left({\displaystyle \frac{k}{b^2}}+{\displaystyle \frac{s}{b}}\right){\displaystyle \frac{1}{\delta }}-{\displaystyle \frac{a}{b}}\left(\lambda -{\displaystyle \frac{ak}{b\gamma }}\right)-{\displaystyle \frac{1}{k}}+\mathcal{O}\left(\delta \right).\hfill \end{array}$$

We note that $y_2>0$ if ${\displaystyle \lambda -\frac{ak}{b\gamma }>0}$, which is true only if the previously encountered condition ${\displaystyle \gamma \lambda -\frac{a}{b}k>0}$ is satisfied. From now on, we will consider this condition valid so that $P_2$ remains in $\mathcal{C}$.

Looking at the Jacobian matrix (21), we obtain in this case:

$$\begin{array}{cc}\hfill J_{11}\left(P_2\right)& =-\gamma y_2=-\gamma \lambda +{\displaystyle \frac{ak}{b}}+\mathcal{O}\left(\delta \right),\hfill \\ \hfill J_{13}\left(P_2\right)& =-a\delta y_2=-{\displaystyle \frac{a}{\gamma }}\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)\delta +\mathcal{O}\left({\delta }^2\right).\hfill \end{array}$$

Now we notice that

$$\begin{array}{cc}& \delta z_2=\left({\displaystyle \frac{k}{b^2}}+{\displaystyle \frac{s}{b}}\right)-\left[{\displaystyle \frac{a}{b}}\left(\lambda -{\displaystyle \frac{ak}{b\gamma }}\right)+{\displaystyle \frac{1}{k}}\right]\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ & {\displaystyle \frac{1}{w_2}}={\displaystyle \frac{b}{k}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ & 2{\displaystyle \frac{b}{k{w}_2}}-1=-1+{\displaystyle \frac{2b^2}{k^2}}\delta +\mathcal{O}\left({\delta }^2\right),\hfill \\ & \delta b{z}_2\left(2{\displaystyle \frac{b}{k{w}_2}}-1\right)=-s-{\displaystyle \frac{k}{b}}+\mathcal{O}\left(\delta \right),\hfill \end{array}$$

so that we have:

$$\begin{array}{cc}\hfill J_{23}\left(P_2\right)& =\delta b{z}_2\left(2{\displaystyle \frac{b}{k{w}_2}}-1\right)=-s-{\displaystyle \frac{k}{b}}+\mathcal{O}\left(\delta \right),\hfill \\ \hfill J_{31}\left(P_2\right)& =a\delta w_2={\displaystyle \frac{ak}{b}}+\mathcal{O}\left(\delta \right),\hfill \\ \hfill J_{32}\left(P_2\right)& =b\delta w_2=k+\mathcal{O}\left(\delta \right),\hfill \\ \hfill J_{33}\left(P_2\right)& =-\delta w_2=-{\displaystyle \frac{k}{b}}+\mathcal{O}\left(\delta \right).\hfill \end{array}$$

Hence, the Jacobian matrix is given by:

$$JF\left(P_2\right)=\left(\begin{array}{ccc}{\displaystyle -\gamma \lambda +\frac{a}{b}k+\mathcal{O}\left(\delta \right)}& 0& {\displaystyle -\frac{a}{\gamma }\left(\gamma \lambda -\frac{a}{b}k\right)\delta +\mathcal{O}\left({\delta }^2\right)}\\ \\ 0& 0& {\displaystyle -s-\frac{k}{b}+\mathcal{O}\left(\delta \right)}\\ \\ {\displaystyle \frac{a}{b}k+\mathcal{O}\left(\delta \right)}& {\displaystyle k+\mathcal{O}\left(\delta \right)}& {\displaystyle -\frac{k}{b}+\mathcal{O}\left(\delta \right)}\end{array}\right).$$

(23)

Denoting by $p_3\left(\rho \right)={\rho }^3+a_1{\rho }^2+a_2\rho +a_3$ its characteristic polynomial, we obtain:

$$\begin{array}{cc}\hfill a_1& ={\displaystyle \frac{1}{b}}(b\gamma \lambda -ak+k)+\mathcal{O}\left(\delta \right)={\displaystyle \frac{k}{b}}+\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)+\mathcal{O}\left(\delta \right),\hfill \\ \hfill a_2& ={\displaystyle \frac{k}{b^2}}(b^{2}s+b\gamma \lambda -ak+bk)+\mathcal{O}\left(\delta \right)=ks+{\displaystyle \frac{k^2}{b}}+{\displaystyle \frac{k}{b}}\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)+\mathcal{O}\left(\delta \right),\hfill \\ \hfill a_3& =k\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)\left(s+{\displaystyle \frac{k}{b}}\right)+\mathcal{O}\left(\delta \right).\hfill \end{array}$$

We apply now the Routh-Hurwitz criterion (see [20]). As we are assuming ${\displaystyle \gamma \lambda -\frac{a}{b}k>0}$, we have $a_i>0$.So we need only to check if $a_1{a}_2>a_3$. We have

$$\begin{array}{cc}\hfill a_1{a}_2& ={\displaystyle \frac{k^{2}s}{b}}+{\displaystyle \frac{k^3}{b^2}}+{\displaystyle \frac{k^2}{b^2}}\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)+ks\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)+{\displaystyle \frac{k^2}{b}}\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)+{\displaystyle \frac{k}{b}}{\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)}^2+\mathcal{O}\left(\delta \right)\hfill \\ & ={\displaystyle \frac{k^{2}s}{b}}+{\displaystyle \frac{k^3}{b^2}}+{\displaystyle \frac{k^2}{b^2}}\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)+k\left(s+{\displaystyle \frac{k}{b}}\right)\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)+{\displaystyle \frac{k}{b}}{\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)}^2+\mathcal{O}\left(\delta \right)\hfill \\ & =a_3+{\displaystyle \frac{k^{2}s}{b}}+{\displaystyle \frac{k^3}{b^2}}+{\displaystyle \frac{k^2}{b^2}}\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)+{\displaystyle \frac{k}{b}}{\left(\gamma \lambda -{\displaystyle \frac{a}{b}}k\right)}^2+\mathcal{O}\left(\delta \right)>a_3\hfill \end{array}$$

as $\delta \to 0^{+}$, so we have asymptotic stability.

The following proposition summarizes the results of this section.

Proposition 11.

When $\delta \to 0^{+}$ we have the following results:

i)

The equilibrium $P_0$ of Case 5 is not in $\overline{\mathcal{C}}$.

ii)

The equilibrium $P_1$ of Case 6 is always unstable.

iii)

The equilibrium $P_2$ of Case 6 is asymptotically stable if ${\displaystyle \gamma \lambda -\frac{a}{b}k>0}$ and it is outside $\overline{\mathcal{C}}$ if ${\displaystyle \gamma \lambda -\frac{a}{b}k<0}$.

Remark 2.

If $a=b=1$, that is ${\delta }_1={\delta }_2={\delta }_3$, we obtain the condition $\gamma \lambda -k>0$ that we found in our previous works.

7. Simulations

We now present the results obtained from simulations conducted using Matlab’s ode45 solver, which utilizes an explicit adaptive Runge-Kutta method. The goal of these simulations was to verify our theoretical analyses. We display our findings in two formats: Cartesian graphs in scenarios where instability is more pronounced and phase diagrams in cases of stability where these representations provide clearer insights.

In the Cartesian graphs shown, the evolution of the variables is depicted by introducing small perturbations to the initial conditions around the equilibrium point. Different colors represent the variables: green for renewable resources, black for non-renewable resources and blue for accumulated wealth. Dashed lines indicate the coordinates of the equilibrium point under consideration, with the color corresponding to each specific component (e.g., blue for wealth). The x-axis is labeled with a generic time unit that does not correspond to specific intervals such as days, months, or years, but instead represents a general time scale that may not cover extended periods. The simulation is terminated once the system’s behavior becomes evident.

The phase diagrams illustrate the behavior of renewable or non-renewable resources and wealth in proximity to the steady states discussed in our study. These diagrams focus on two main variables, plotting three distinct trajectories (depicted in magenta, dark green and red) generated by applying three different perturbations around the critical point within a two-dimensional space. The graphs highlight stable points, where the dynamics naturally lead the system towards equilibrium. Initial points are marked with small circles, while steady states are indicated by small stars.

7.1. Case 3

We wanted to test a stable situation with ${\displaystyle \gamma \lambda -\frac{a}{b}k<0}$. We considered the following parameter values: $k=1,s=1,\gamma =0.01,\lambda =10,{\delta }_{\ell }=0.01,\ell =1,2,3$, that is $a=b=1$, resulting in the stable point $P_2=(0,198.9898,98.9898)$ with ${\displaystyle \gamma \lambda -\frac{ak}{b}=-0.9}$. In Figure 1 we have represented the phase diagrams obtained by considering non-renewable resources versus wealth on the left, and renewable resources versus wealth on the right. The Jacobian matrix calculated at $P_2$ has eigenvalues ${\rho }_{1,2}=-0.4949\pm 1.2981i,{\rho }_3=-0.8899$.

We also tested, without displaying the graph, some cases of instability, noting that the trajectories, in this situation, stabilize at the point denoted as $P_2$ in Case 6. Specifically, by varying the parameter $\lambda =150$ in Case 3, we obtain ${\displaystyle \gamma \lambda -\frac{ak}{b}=0.5}$ with a critical point $P_2=(0,1999.0,999.0)$ and convergence to $(50.1,1948.9,999.0)$, which corresponds to the point $P_2$ with the same parameters as in Case 6.

7.2. Case 4

In the simulations for the scenario described in Case 4, we considered the parameters $\gamma =0.01,k=1,\lambda =10,{\delta }_1={\delta }_2={\delta }_3=1$, starting from the critical point $P=(10,Z,0)$ with $Z=5$; since we want to study ${\displaystyle Z<\frac{s}{{\delta }_2}-\frac{{\delta }_1}{{\delta }_2}\lambda }$, we have set $s=20$.

In Figure 2, we show the simulations related to $P=(10,5,0)$; on the left, we started with a positive perturbation on the coordinate $w=0$, while applying a negative perturbation to $y=\lambda$ and $Z=5$. In the Figure 2 on the right, all perturbations are positive. As predicted by the theory, the trajectories of the variables $z\left(t\right)$ and $w\left(t\right)$ decrease, while $y\left(t\right)\to \lambda$. We note that for both graphs, we have convergence $y\left(t\right)\to \lambda ,w\left(t\right)\to 0^{+}$, while the trajectory of $z\left(t\right)$ behaves differently. Furthermore, the convergence of the variable $z\left(t\right)$ varies depending on the magnitude of the input perturbation and its relationship with Z. The Jacobian matrix calculated at the point has the following eigenvalues: ${\rho }_1=-0.1,{\rho }_2=0,{\rho }_3=-5$.

7.3. Case 6

Figure 3 shows the simulations related to the point $P_2=(50.1,1948.9,999.0)$ obtained with $k=1,s=1,\gamma =0.01,\lambda =150,{\delta }_{\ell }=0.001,\ell =1,2,3$, i.e. $a=b=1$ and ${\displaystyle \gamma \lambda -\frac{ak}{b}=0.5}$, which is a stable point. The eigenvalues of the Jacobian matrix calculated here are ${\rho }_1=-0.4866,{\rho }_{2,3}=-0.5067\pm 1.3205i$. We note that this is a stable situation where all variables are strictly positive.

8. Conclusions

We now summarize the main results of this paper. First, we have obtained that for the trajectories in $\mathcal{C}$ it can’t happen $y\left(t\right)\to +\infty$ or $z\left(t\right)\to +\infty$ or $w\left(t\right)\to +\infty$, as $t\to +\infty$. Of course, this does not mean that this variables remain bounded, but that there is a positive constant C such that if, for example, $w>C$ at a time $T>0$, then w must go below C at a time $T_1>T$. So an excessive growth in wealth generates oscillations, and these are phenomena that generally have negative social consequences. Another relevant result is the following: although we have obtained that many equilibria of the system are unstable (or are outside $\overline{\mathcal{C}}$), there are still, at least for some range of the parameters, stable equilibria with positive values for all the variables (see Proposition 11 and Figure 3). We note that all numerical simulations confirm the theoretical results obtained. From these results a general indication seems to emerge: assuming the hypotheses of the model, an indefinite growth in wealth is not a reasonable aim, while it seems possible to drive the society towards a steady state with positive values of all the variables, at least for small rates of exploitation of resources. Hence, the results of our works (the present one and the previous papers [1,2,3]) are in line with the ideas of those who support the need to abandon the idea of infinite growth and instead to search a path towards stationary economic situations, possibly through a period of degrowth.

These results are similar to those that we have obtained in our previous papers, using similar but different models. This seems to us a relevant achievement in itself, because it means that these results are very stable with respect to changes in the model. Despite the fact that, in any case, all the models that we have studied, in the present and the previous papers, are of course too simple to allow precise predictions on social evolution, this stability of the qualitative ideas that they suggest may be an interesting contribution to the debate on the future of our societies.