On an Extended Time-Varying Beverton–Holt Equation Subject to Harvesting Monitoring and Population Excess Penalty

Abstract

This paper considers a more general eventually time-varying Beverton–Holt equation for species evolution which can include a harvesting action and a penalty for overpopulation numbers. The harvesting action may be positive (typically consisting of hunting or fishing) or negative which refers to repopulation within the environment. One considers also a penalty of quadratic type on the overpopulation and the introduction of a term related to Allee effect to take account of small levels of population. The intrinsic growth rate is assumed either to exceed unity or to be under unity. In the second case, the extinction point is a locally stable attractor while the other positive equilibrium point is unstable contrarily to the commonly studied case of intrinsic growth rate exceeding unity where the above roles are inverted. This consequence implies that the extinction point is also globally asymptotically stable for any given finite initial condition. In the case when the eventual overpopulation is penalized with a sufficiently large coefficient which exceeds a prescribed threshold, to quantify such an excess, only a globally asymptotically stable extinction attractor is present and no other positive equilibrium points exist. In the case of a positive moderate quadratic evaluation term for such an overpopulation, one or two positive equilibrium points coexist with the extinction one. The smaller one is unstable contrarily to the extinction equilibrium which is locally asymptotically stable. If it exists a second largest positive equilibrium point, being distinct to the above-given one, then it can be unstable or locally stable depending on the parameterization. Also, some methods of monitoring the population evolution through control laws on the harvesting action are discussed.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Beverton–Holt equation is an useful discrete equation for modelling the evolution of species which reproduce by eggs such as birds, fishes and insects [1]. It can be considered the counterpart of the logistic equation in the Verhulst´s continuous-time model. The basic Beverton–Holt model is parameterized by two positive sequences, namely, the carrying capacity of the environment which depends on resources availability, temperature, humidity, etc., and the intrinsic growth rate which is associated with the species reproduction capability, the survivorship chance etc. The intrinsic growth rate has typically to exceed unity to avoid extinction. In the most general case, those parameters can be changed to sequences to describe potential different behaviors of the population evolution in different periods, for instance, seasonality. There are other two typical parameters to be eventually considered which generalize the model such as the independent consumption which describes recruitment variations due to unforeseen disturbances and eventual repopulation or interchange of population with neighbour environments and the harvesting process associated with fishing or hunting and which depends on regulation based on the available spawning stock and foreseen recruitment. The basic time-invariant Beverton–Holt equation has two equilibrium points which are the extinction point, which is locally unstable, and the carrying capacity level which is globally asymptotically stable. The so-called Cushing–Henson conjecture established that, if the equation is modelled by periodic parameterizing sequences of carrying capacities and intrinsic growth rates, then the averaged periodic sequence of population lies below the average of the corresponding average of the carrying capacities. The conjecture has been rigorously proved to be true by Stevic [2]. Some extensions of the basic model concerning the Cushing–Henson conjectures for the Beverton–Holt q-difference equation have been discussed in [3]. A control theory point of view on the Beverton–Holt equation has been adopted and discussed in [4–6] while giving a design procedure of the environment carrying capacity for monitoring the suitable sequence of values to follow of the population evolution. The applicability of the proposed method is claimed for semiopen environments such like certain fisheries. It is discussed in [7] how, in practice, the intrinsic growth rate can be dependent on the environment carrying capacity. Also, it is discussed in [8] an impulsive extended competition Beverton–Holt model between species from the stability point of view. The usefulness of the Beverton–Holt and other mathematical models in maritime biology is described in [9]. In [10–13], the harvesting action is investigated in an extended Beverton–Holt model. Normally, harvesting refers to fishing or hunting which is subjected to authorities regulation but it can also be total or partially illegal while associated to furtive uncontrolled actions. See also some references therein and [6, 8]. Other related studies consider alternative generalizations concerned with periodic behaviors, associated, for instance, to seasonality [14, 15], global dynamics analysis of some extended equation versions [16], presence of bifurcations [17, 18] or resonances [19], and perturbations of the basic model. See, for instance [18]. On the other hand, an extended Beverton model on isolated time scales is analyzed in detail in [20]. Also, an extension of the Beverton–Holt model including discrete delays in the evolution dynamics has been investigated in [21]. A Beverton–Holt model extension including discrete delays in the evolution dynamics has been investigated in [21]. On the other hand, it can be pointed out that Beverton–Holt-based models are used also by biologists when monitoring fishing stock availability and fishing migrations to evaluate the recommended maximum number of captures (or recommended harvesting action) to avoid the environment degradation and species extinction. See, for instance [22, 23], and some of the references therein.

In this paper, we focus on a generalized Beverton–Holt equation which assumes a quadratic-type penalty for the population excess describing the potential internal competence between the individuals for food, refuge, etc. The harvesting action is considered jointly with eventually present independent consumption if necessary. It is seen that the presence of such a term can translate into the presence of two other equilibrium points. The paper also designs species evolution control laws by monitoring the harvesting action and the influence in the results of considering a modelling function of Allee’s effect which makes difficult growing or even can cause extinction for small numbers of reproductive individuals.

The paper is organized as follows. Section 2 deals with the equilibrium points in the presence and absence of harvesting action, considered together with eventual independent consumption, in the case when the intrinsic growth rate elements exceed unity. The harvesting sequence can have positive, negative or null elements. The local asymptotic stability of each feasible (that is being real and non-negative) equilibrium points is characterized in the case when the parameterizing sequences converge to limits. Section 3 develops two methods to derive control laws for the harvesting actions again if the intrinsic growth rate exceeds unity for all time. The first proposed method is based on the convergence of the solution sequence of the population to a prescribed targeted equilibrium point of the population value by choosing the harvesting sequence. Classical criteria for convergence of sequences, such as D’Alembert, Cauchy, and Raabe criteria, are involved in the respective monitoring rules of the harvesting action. The second method relies on a sample-to-sample monitoring of the solution sequence to target a prescribed evolution pattern by designing the harvesting sequence. Next, Section 4 relies on introducing Allee’s effect to modify the basic Beverton–Holt equation to describe the situation arising under small numbers of individuals which make difficult the reproductive action and can lead to extinction even the intrinsic growth rate exceeds unity for all time. The equilibrium points, their stability conditions as well as extinction conditions are investigated if the intrinsic growth rate exceeds unity for all time. The second part of this section proposes a penalty term in the Beverton–Holt equation for high levels of population in the absence of harvesting. The resulting equilibrium points and their stability issues are also investigated if the intrinsic growth rate exceeds unity for all time. On the other hand, Section 5 relies on extinction conditions and the local asymptotic stability of the extinction equilibrium point if the intrinsic growth sequence has elements being less than unity in the absence of harvesting by considering the modified model with quadratic penalty time for high levels of population. The local asymptotic stability of the positive equilibrium points is also investigated. Section 6 is devoted to discuss some numerical examples and, finally, some conclusions end the paper. In the following, the subsequent notation is used:

  ${\mathbf{Z}}_{0+}$ and ${\mathbf{Z}}_{0-}$ denote, respectively, the sets of non-negative and nonpositive integer numbers

  ${\mathbf{Z}}_{+}$ and ${\mathbf{Z}}_{-}$ denote, respectively, the sets of positive and negative integer numbers

  ${\mathrm{R}}_{0+}$ and ${\mathrm{R}}_{0-}$ denote, respectively, the sets of non-negative and nonpositive real numbers

  ${\mathrm{R}}_{+}$ and ${\mathrm{R}}_{-}$ denote, respectively, the sets of positive and negative real numbers

2. Equilibrium Points and Their Local Asymptotic Stability under Positive, Null, or Negative Harvesting

“Harvesting” is referred to fishing and hunting subject to administrative regulation which depends on the species population. On the other hand, “independent consumption” refers to positive of negative supplies of extra populations due to migrations from outside of the environment under consideration [10, 13, 14]. In the sequel, we consider both effects together integrated in the same additive perturbation sequence to that of the standard population evolution sequence. Consider the following, in general, time-varying Beverton–Holt equation subject to harvesting action eventually combined with independent consumption:$$\begin{array}{}\text{(1)}& \begin{array}{cc}{y}_{n+1}=\frac{a_n{K}_n{y}_n}{K_n+a_n-1y_n}-h_n\text{,}& \forall n\in {\mathbf{Z}}_{0+}\end{array}\text{,}\end{array}$$with initial condition $y_0\ge 0$, where ${a_n}_{n=0}^{\infty }$ is the intrinsic growth rate of the species and ${K_n}_0^{\infty }$ is the environment carrying capacity sequence, both being positive real sequences. If a constraint $y_{n+1}\ge z_{n+1}\ge 0$; $\forall n\in {\mathbf{Z}}_{0+}$ is prefixed for some sequence ${z_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}$ so that the harvesting sequence ${h_n}_{n=0}^{\infty }$ has to fulfil:$$\begin{array}{}\text{(2)}& \begin{array}{cc}{h}_n\in -\infty ,\frac{a_n{K}_n{y}_n}{K_n+a_n-1y_n}-z_{n+1}\text{,}& \forall n\in {\mathbf{Z}}_{0+}\end{array}\text{.}\end{array}$$

The harvesting sequence defines the population amount which is not related to the dynamics evolution within the environment because of natural reproduction and dead concerns. It is related to an increase or decrease of individuals due to population flux either from or to the habitat plus eventual decrease of population due to hunting or fishing. In this way, the sequence can take negative values at a particular sampling instant because of the sign in (1), this situation will correspond to an increase of the amounts of individuals), positive values (that is, decrease of population), or zero (that is, the population is just modified by the natural reproduction and dead within the considered habitat). According to that philosophy, the harvesting is considered in this paper as the eventual combination of an eventual traditional harvesting (that is, hunting/fishing) and eventual migrations in both senses from or to the habitat under study. Also, the hunting or fishing includes, in general, legal or illegal actions (poaching).

It turns out that any equilibrium point needs to be real non-negative in order to be feasible (that is, real and positive) as it is addressed in the subsequent result.

Theorem 1.

Assume that ${a_n}_0^{\infty }\subset 1,\overline{a}⟶a$, ${K_n}_0^{\infty }\subset 0,\overline{K}⟶K$, ${z_n}_0^{\infty }\subset 0,\overline{z}⟶z$ and ${h_n}_0^{\infty }⟶h$ with $h_n\in -\infty ,a_n{K}_n{y}_n/K_n+a_n-1y_n-z_{n+1}$; $\forall n\in {\mathbf{Z}}_{0+}$. Then, the solution of (1) has:

(i) A real positive equilibrium point $\overline{y}=K>0$ and a null equilibrium point (extinction) at ${\overline{y}}_0=0$ if $h=0$.

(ii) If $h\ge K$ then there is no nonextinction feasible equilibrium point and the nonextinction equilibrium points are feasible if $h<K$.

(iii) If $h\ne 0$ then the potential nonextinction equilibrium points ${\overline{y}}_1>0$ and ${\overline{y}}_2\ge {\overline{y}}_1>0$ are given by$$\begin{array}{}\text{(3)}& {\overline{y}}_{\mathrm{1,2}}=a-1K-h\mp \frac{\sqrt{{a-1}^2{K-h}^2-4a-1Kh}}{2a-1}\text{,}\end{array}$$

  which are real if and only if $h\in -\infty ,a+1-2\sqrt{a}/a-1K\cup a+1+2\sqrt{a}/a-1K,+\infty$. The equilibrium point ${\overline{y}}_2$ is feasible if and only if $h\in -\infty ,K$ which restricts the above-given inequality for realness. Also, the equilibrium point ${\overline{y}}_1$ is feasible if and only if $h\in 0,K$. If $h=Ka+1\mp 2\sqrt{a}/a-1$ then ${\overline{y}}_1={\overline{y}}_2>0$. In terms of the intrinsic growth rate, the nonextinction equilibrium points ${\overline{y}}_1$ and ${\overline{y}}_2$ are both feasible if $h\in 0,K$ and $a-1{K-h}^2\ge 4Kh$, equivalently, if $a\ge {K+h/K-h}^2$ ($a>1$ if $h=0$). Also, ${\overline{y}}_1$ is not feasible for $h<0$ and ${\overline{y}}_2$ is feasible for $h<0$ irrespective of $a>1$ and $K>0$.

Proof.

Note directly that the extinction level ${\overline{y}}_0=0$ is an equilibrium point. Also, by replacing in (1) the limits of the various sequences, one gets a single root $\overline{y}=K$ if $h=0$ and, if $h\ne 0$, then$$\begin{array}{}\text{(4)}& a-1{\overline{y}}^2+a-1h-K\overline{y}+Kh=0\text{.}\end{array}$$

Thus, since $a>1$, if $h\ge K$ then there is no real non-negative solution to (1) since (4) fails for $\overline{y}>0$ and $h\ge K$.

If $0<h<K$, or if $h<0$, then the roots of (2) are$$\begin{array}{c}\text{(5)}& {\overline{y}}_1=\frac{a-1K-h-\sqrt{{a-1}^2{K-h}^2-4a-1Kh}}{2a-1}\text{,}{\overline{y}}_2=\frac{a-1K-h+\sqrt{{a-1}^2{K-h}^2-4a-1Kh}}{2a-1}\text{.}\end{array}$$

Properties [(i)-(ii)] have been proved. On the other hand, the nonextinction equilibrium points ${\overline{y}}_1$ and ${\overline{y}}_2$ are both feasible for $h\in 0,K$ if $a-1{K-h}^2\ge 4Kh$, equivalently, if $a\ge K+h/{K-h}^2$ ($a>1$ if $h=0$). Also, ${\overline{y}}_1$ is not feasible for $h<0$ from (5) and ${\overline{y}}_2$ is feasible, also from (5), for $h<0$ irrespective of $a>1$ and $K>0$.

Note that, in order for the roots of (4) to be real, since $a>1$, it is needed that $\theta h\ge 0$, where$$\begin{array}{c}\text{(6)}& \theta h=a-1{K-h}^2-4Kh=a-1h^2+a-1K^2-2a-1Kh-4Kh=a-1h^2+a-1K^2-2aKh+2Kh-4Kh=a-1h^2+a-1K^2-2Kha+1\text{.}\end{array}$$

Note that, ${\overline{y}}_{\mathrm{1,2}}$ are real if and only if the zeros of $\theta h$, that is, $h_{\mathrm{1,2}}=Ka+1\mp 2\sqrt{a}/a-1$, are non-negative real since $a>1$. Since $\theta h$ is a convex parabola, $\theta h\ge 0$, and both equilibrium points ${\overline{y}}_1$ and ${\overline{y}}_2$ are real (and they can be positive) if and only if $h\in -\infty ,Ka+1-2\sqrt{a}/a-1\cup Ka+1+2\sqrt{a}/a-1,+\infty$. Contrarily, if $h\in Ka+1-2\sqrt{a}/a-1,Ka+1+2\sqrt{a}/a-1$ then the zeros are not real and the nonextinction equilibrium points ${\overline{y}}_{\mathrm{1,2}}$ never exist. Since, their feasibility implies that $h\le K$ then the equilibrium point ${\overline{y}}_2$ is feasible if and only if $h\in -\infty ,K$. Since $a+1-2\sqrt{a}/a-1<1$ for $a>1$. Also, the equilibrium point ${\overline{y}}_1$ is feasible if and only if $h\in 0,Ka+1-2\sqrt{a}/a-1$. It is obvious that $h=h_{\mathrm{1,2}}=Ka+1\mp 2\sqrt{a}/a-1$ then ${\overline{y}}_1={\overline{y}}_2$. Property (iii) is proved.

The use of the inverse sequence of that of a the population evolution sequence is of interest to derive easily some interesting results concerning the stability and the asymptotic boundedness of the solution as it is addressed in the subsequent result:

Theorem 2.

Define the inverse sequence of the solution of (1) as $x_n=y_n^{-1}$; $\forall n\in {\mathbf{Z}}_{0+}$. The following properties hold:

(i) The inverse sequence ${x_n}_{n=0}^{\infty }$ of the solution ${y_n}_{n=0}^{\infty }$ is given by the discrete equation:$$\begin{array}{}\text{(7)}& x_{n+1}={\mu }_n{x}_n+{\nu }_n+h_n^I\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$

  where ${\mu }_n=a_n^{-1}$; ${\nu }_n=1-{\mu }_n{K}_n^{-1}=a_n-1a_n^{-1}{K}_n^{-1}$, subject to $h_n\in -\infty ,a_n{K}_n{y}_n/K_n+a_n-1y_n$; $\forall n\in {\mathbf{Z}}_{0+}$, and $h_n^I$ is zero if $h_n=0$ for any $n\in {\mathbf{Z}}_{0+}$, with$$\begin{array}{}\text{(8)}& \begin{array}{cc}{h}_n=\frac{a_n{K}_n{h}_n^I}{21-a_n^{-1}+h_n^I+a_n^{-1}{x}_n{K}_n{x}_n+K_n^{-1}{a}_n+a_n^{-1}-2+h_n^I{a}_n-1}\text{;}& \forall n\in {\mathbf{Z}}_{0+}\end{array}\text{,}\end{array}$$

  for $x_0=y_0^{-1}>0$. The solution is equivalently expressed from given initial conditions as follows:$$\begin{array}{}\text{(9)}& x_{n+1}={\prod }_{i=0}^n{\mu }_i{x}_0+{\sum }_{i=0}^n{\prod }_{j=i+1}^n{\mu }_j{\nu }_i+h_i^I\text{.}\end{array}$$

(ii) ${x_n}_{n=0}^{\infty }$ is bounded, equivalently, ${y_n}_{n=0}^{\infty }$ does not vanish neither at any sample nor asymptotically (and then the population does not extinguish either in finite time or asymptotically) if ${\lim \sup }_{n⟶\infty }{\sum }_{i=0}^n{\prod }_{j=i+1}^n{\mu }_j{\nu }_i+h_i^I<\infty$. In particular, if ${a_n}_{n=0}^{\infty }\subset \underset{\_}{a},\overline{a}$ and ${\lim \inf }_{n⟶\infty }{a}_n\ge \underset{\_}{a}>1$ then$$\begin{array}{}\text{(10)}& \begin{array}{cc}{x}_n\le {\underset{\_}{a}}^{-n}{x}_0+\frac{1-{\underset{\_}{a}}^{-n}}{1-\underset{\_}{a}}\underset{0\le i\le n-1}{\sup }{\nu }_i+h_i^I\text{;\hspace{0.17em}}& \forall n\in {\mathbf{Z}}_{+}\text{,}\end{array}\end{array}$$

  is bounded for any given finite $x_0\ge 0$ for all $n\in {\mathbf{Z}}_{0+}$ if ${\sup }_{0\le i\le n-1}{\nu }_i+h_i^I<+\infty$. For any finite $x_0>0$ , so that $y_0=x_0^{-1}$, the sequence ${y_n=x_n^{-1}}_{n=0}^{\infty }$ satisfies:$$\begin{array}{}\text{(11)}& y_n\ge \frac{1-\underset{\_}{a}}{{\underset{\_}{a}}^{-n}1-\underset{\_}{a}{x}_0+1-{\underset{\_}{a}}^{-n}\underset{0\le i\le n-1}{\sup }{\nu }_i+h_i^I}>0\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{+}\text{.}\end{array}$$

(iii) Define $b_n=1-{\mu }_n=a_n-1/a_n$ by assuming that ${a_n}_{n=0}^{\infty }\subset 1,\infty$ subject to ${\sum }_{n=0}^{\infty }{b}_n={\sum }_{n=0}^{\infty }{a}_n-1/a_n=\infty$. Assume that $h_n^I+a_n-1/K_n{a}_n\le \epsilon a_n-1/a_n$ for any integer $n\ge n_0$, some $n_0\in {\mathbf{Z}}_{0+}$ and some $\epsilon \in {\mathbf{R}}_{+}$. Then, ${\lim \sup }_{n⟶\infty }{x}_n\le \epsilon \iff {\lim \inf }_{n⟶\infty }{y}_n\ge {\epsilon }^{-1}$. The constraint $h_n^I+a_n-1/K_n{a}_n\le \epsilon a_n-1/a_n$ is satisfied under any of the subsequent stipulations for each $n\in {\mathbf{Z}}_{0+}$:

(1) $0\le h_n^I\le \epsilon -K_n^{-1}{a}_n-1/a_n$ requiring that the harvesting sequence $h_n\ge 0$ and the carrying capacity $K_n\ge 1/\epsilon$

(2) $h_n^I<0$ and $a_n-1/a_n{K}_n^{-1}\ge h_n^I\ge K_n^{-1}-\epsilon a_n-1/a_n\le h_n^I<0$ requiring that $h_n<0$ and $K_n<1/\epsilon$

(3) $h_n^I<0$ and $a_n-1/K_n{a}_n<h_n^I\le a_n-1/a_n\epsilon +K_n^{-1}$ requiring that $h_n<0$

(iv) The extinction equilibrium point ${\overline{y}}_0=0$ is unstable. The two positive equilibrium points ${\overline{y}}_{\mathrm{1,2}}$ in (5) arising when the parametrical sequences ${a_n}_{n=0}^{\infty }⟶a$, ${K_n}_{n=0}^{\infty }⟶K$, ${h_n}_{n=0}^{\infty }⟶h$ fulfil the following theories:

(a) ${\overline{y}}_2>{\overline{y}}_1$ is jointly feasible and locally asymptotically stable if $-\infty <h<\sqrt{a}-1/\sqrt{a}+1K$

(b) If ${\overline{y}}_1<{\overline{y}}_2$ is feasible then it is not locally asymptotically stable

(c) If $h=\sqrt{a}-1/\sqrt{a}+1K$ then the unique nonextinction equilibrium point ${\overline{y}}_1={\overline{y}}_2=K-\mathrm{h}/2=K/\sqrt{a}+1$ is jointly feasible and locally asymptotically stable

Proof.

One obtains (7) directly from the following equivalent expression to (1):$$\begin{array}{c}\text{(12)}& x_{n+1}=\frac{1}{y_{n+1}}=\frac{1}{a_n{K}_n{y}_n/K_n+a_n-1y_n-h_n}=\frac{K_n+a_n-1y_n}{a_n{K}_n{y}_n-h_n{K}_n+a_n-1y_n}\begin{array}{cc}=a_n^{-1}{x}_n-a_n^{-1}-1K_n^{-1}+h_n^I\text{;\hspace{0.17em}}& \forall n\in {\mathbf{Z}}_{0+}\end{array}\text{,}\end{array}$$subject to $h_n\in -\infty ,a_n{K}_n{y}_n/K_n+a_n-1y_n$; $\forall n\in {\mathbf{Z}}_{0+}$ for keeping the sample-to-sample non-negativity of ${y_n}_{n=0}^{\infty }$, where $x_n=y_n^{-1}$; $\forall n\in {\mathbf{Z}}_{0+}$, and$$\begin{array}{}\text{(13)}& h_n^I=x_{n+1}-a_n^{-1}{x}_n+a_n^{-1}-1K_n^{-1}=\frac{K_n{x}_n+a_n-1}{a_n{K}_n-g_n}-a_n^{-1}{x}_n+a_n^{-1}-1K_n^{-1}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$with $g_n=h_n{K}_n{x}_n+a_n-1$; $\forall n\in {\mathbf{Z}}_{0+}$. One gets from (13) that$$\begin{array}{}\text{(14)}& K_n^{-1}+h_n^I-a_n^{-1}{K}_n^{-1}-x_n{g}_n=a_n{K}_n{h}_n^I\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$which leads to$$\begin{array}{}\text{(15)}& g_n=\frac{a_n{K}_n{h}_n^I}{K_n^{-1}+h_n^I-a_n^{-1}{K}_n^{-1}-x_n}=h_n{K}_n{x}_n+a_n-1\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$so that$$\begin{array}{c}\text{(16)}& h_n=\frac{a_n{K}_n{h}_n^I}{K_n^{-1}+h_n^I-a_n^{-1}{K}_n^{-1}-x_n{K}_n{x}_n+a_n-1}=\frac{a_n{K}_n{h}_n^I}{2x_n+K_n^{-1}{a}_n+a_n^{-1}-2+h_n^I{K}_n{x}_n+a_n-1-2a_n^{-1}{x}_n+a_n^{-1}{K}_n{x}_n^2}\\ =\frac{a_n{K}_n{h}_n^I}{21-a_n^{-1}+h_n^I+a_n^{-1}{x}_n{K}_n{x}_n+K_n^{-1}{a}_n+a_n^{-1}-2+h_n^I{a}_n-1}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$which leads directly to (8). Equation (9) follows directly from recursive calculations with (7). Property (i) has been proved. Property (ii) is a direct consequence of Property (i) since ${\prod }_{i=0}^{n+1}{\mu }_i<{\prod }_{i=0}^n{\mu }_i<1$; $\forall n\in {\mathbf{Z}}_{0+}$ and ${\lim }_{n⟶\infty }{\prod }_{i=0}^n{\mu }_i=0$ since ${\lim \inf }_{n⟶\infty }{a}_n\ge \underset{\_}{a}>1$.

To prove Property (iii), we rewrite an upper-bounding expression of (7) as$$\begin{array}{}\text{(17)}& \begin{array}{cc}{x}_{n+1}=1-b_n{x}_n+\frac{b_n}{K_n}+h_n^I\le 1-b_n{x}_n+\frac{b_n}{K_n}+h_n^I\le 1-b_n{x}_n+{\epsilon }_n{b}_n\text{;\hspace{0.17em}}& \forall n\in {\mathbf{Z}}_{0+}\end{array}\text{,}\end{array}$$since ${a_n}_{n=0}^{\infty }\subset 1,\infty$ implies that $b_n=1-{\mu }_n=a_n-1/a_n\in 0,1$; $\forall n\in {\mathbf{Z}}_{0+}$ subject to ${\sum }_{n=0}^{\infty }{b}_n={\sum }_{n=0}^{\infty }{a}_n-1/a_n=\infty$. Since it is also assumed that $h_n^I+b_n/K_n\le \epsilon b_n=\epsilon a_n-1/a_n$; $\forall n\in {\mathbf{Z}}_{0+}$, this constraint is achieved if $h_n^I+b_n/K_n=h_n^I+a_n-1/K_n{a}_n\le {\epsilon }_n{b}_n\le \epsilon a_n-1/a_n$; $\forall n\ge n_0\in {\mathbf{Z}}_{0+}$, and some $n_0\in {\mathbf{Z}}_{0+}$, where ${{\epsilon }_n}_{n=0}^{\infty }\subset 0,\epsilon$ for any integer $n\ge n_0$, some $n_0\in {\mathbf{Z}}_{0+}$ and some $\epsilon \in {\mathbf{R}}_{+}$ then it follows from [Lemma 1.2 (i), [24]] that $0\le \underset{n⟶\infty }{\lim \sup }{x}_n\le \epsilon$. The stipulation 1 is got from the constraints $h_n^I\ge 0$ (implying that $h_n\ge 0$) and then $h_n^I+a_n-1/K_n{a}_n=|h_n^I+a_n-1/K_n{a}_n|\le \epsilon a_n-1/a_n$, $\forall n\in {\mathbf{Z}}_{0+}$.

The stipulation 2 follows for $h_n^I<0$ (implying that $h_n<0$) and $h_n^i\le b_n{K}_n^{-1}$ so that the subsequent constraint holds $h_n^I+a_n-1/K_n{a}_n=a_n-1/K_n{a}_n-h_n^I$; $\forall n\in {\mathbf{Z}}_{0+}$.

The stipulation 3 follows for $h_n^I<0$ and $h_n^i\le b_n{K}_n^{-1}$ so that the following constraint holds $h_n^I+a_n-1/K_n{a}_n=h_n^I-a_n-1/K_n{a}_n$; $\forall n\in {\mathbf{Z}}_{0+}$.

Property (iii) has been proved. Property (iv) follows from a local perturbation analysis. It is assumed that the invariant equation (1) perturbed from any equilibrium point $\overline{y}$ as $y_n=\overline{y}+\delta y_n$; $\forall n\in {\mathbf{Z}}_{0+}$. The linearized perturbation transmitted to the next sample is $\overline{y}+\delta {\overline{y}}_{n+1}=aK/K/\overline{y}+\delta y_n+a-1-h$; $\forall n\in {\mathbf{Z}}_{0+}$ which leads to $\delta y_{n+1}\le K^{2}a/{K+\overline{y}a-1}^2\delta y_n+o\delta y_n$; $\forall n\in {\mathbf{Z}}_{0+}$ and, for sufficiently small $\delta y_n$, $\delta y_{n+1}/\delta y_n<1$, $\forall n\in {\mathbf{Z}}_{0+}$, so that $\overline{y}$ is locally asymptotically stable, if and only if $K\sqrt{a}/K+\overline{y}a-1<1$. Equivalently, if and only if $\overline{y}>K\sqrt{a}-1/a-1=K/\sqrt{a}+1$. Contrarily, the equilibrium point $\overline{y}$ is not locally asymptotically stable (that is, either critically stable or unstable) if and only if $K\sqrt{a}/K+\overline{y}a-1\ge 1$, that is, if and only if $\overline{y}\le K/\sqrt{a}+1$. In particular, it is unstable if $K\sqrt{a}/K+\overline{y}a-1>1$, that is, if $\overline{y}<K/\sqrt{a}+1$. The local asymptotic stability constraint fails and the instability constraint holds if $\overline{y}={\overline{y}}_0=0$ (extinction equilibrium point) since $a>1$.

Then, the extinction equilibrium point is unstable.

For addressing the local asymptotic stability of the other two equilibrium points $\overline{y}={\overline{y}}_1$ and $\overline{y}={\overline{y}}_2$, provided they are feasible and distinct (that is, the radicand of (5) is real positive), note that equilibrium points are locally asymptotically stable if and only if$$\begin{array}{}\text{(18)}& {\overline{y}}_{\mathrm{1,2}}=\frac{a-1K-h\mp \sqrt{{a-1}^2{K-h}^2-4a-1Kh}}{2a-1}>\frac{K}{\sqrt{a}+1}\text{.}\end{array}$$

Considering ${\overline{y}}_2$ under the above constraint, that one becomes equivalent to$$\begin{array}{}\text{(19)}& \sqrt{{a-1}^2{K-h}^2-4a-1Kh}>\sqrt{a}-12K-\sqrt{a}+1K-h\text{.}\end{array}$$

The above-given constraint (19) holds if $u>v$, where$$\begin{array}{c}\text{(20)}& u=a-1a-1{K-h}^2-4Kh=a-1a{K-h}^2-{K+h}^2\\ ={\alpha }_1{h}^2+{\beta }_{1}h+{\gamma }_1\text{,}\end{array}$$where$$\begin{array}{}\text{(21)}& {\alpha }_1={a-1}^2{\beta }_1=-2K{a}^2-1;{\gamma }_1={a-1}^2{K}^2\text{,}\end{array}$$and$$\begin{array}{c}\text{(22)}& v={\sqrt{a}-12K-\sqrt{a}+1K-h}^2=a+1-2\sqrt{a}4K^2+a+1+2\sqrt{a}{K}^2+h^2-2Kh-4K\sqrt{a}+1K-h\\ ={\alpha }_2{h}^2+{\beta }_{2}h+{\gamma }_2\text{,}\end{array}$$where$$\begin{array}{}\text{(23)}& {\alpha }_2={\alpha }_1={a-1}^2;{\beta }_2=2Ka-12\sqrt{a}-a-1;{\gamma }_2=K^2{a}^2+6a-4a\sqrt{a}-4\sqrt{a}+1\text{.}\end{array}$$

Thus, $u>v$, as necessary condition $a\ge {K+h/K-h}^2$, equivalently, $h\le \sqrt{a}-1/\sqrt{a}+1K$ in order for the radicand in the definition of ${\overline{y}}_{\mathrm{1,2}}$ to be non-negative (equilibrium feasibility condition), with $a>{K+h/K-h}^2$ if $h\ne K\sqrt{a}-1/\sqrt{a}+1$ and $h\ne 0$. For $h=0$, ${\overline{y}}_2=K$ is locally asymptotically stable from $u>v$ since $\sqrt{a}+1>\sqrt{a}-1$ holds trivially. Now (20)–(23), one gets that ${\overline{y}}_2$ is locally asymptotically stable, that is, $u>v$, if and only if:$$\begin{array}{}\text{(24)}& h<\frac{{\gamma }_1-{\gamma }_2}{{\beta }_2-{\beta }_1}=\frac{{a-1}^2-a^2+6a-4a\sqrt{a}-4\sqrt{a}+1}{2a-12\sqrt{a}-a-1+2a^2-1}K=\frac{a+1\sqrt{a}-2a}{a-1\sqrt{a}}K=\frac{\sqrt{a}-1}{\sqrt{a}+1}K\text{,}\end{array}$$which coincides with the feasibility condition. Thus, ${\overline{y}}_2$ is both feasible and locally asymptotically stable if and only if $-\infty <h<\sqrt{a}-1/\sqrt{a}+1K$.

For the local asymptotic stability of ${\overline{y}}_1$, equation (19) becomes modified as follows:$$\begin{array}{}\text{(25)}& -\sqrt{{a-1}^2{K-h}^2-4a-1Kh}>\sqrt{a}-12K-\sqrt{a}+1K-h\text{,}\end{array}$$so that, one gets the two associated constraints:

(1) $\sqrt{a}-12K-\sqrt{a}+1K-h\le 0$

(2) $\sqrt{u}=\sqrt{{a-1}^2{K-h}^2-4a-1Kh}<\sqrt{v}=-\sqrt{a}-12K-\sqrt{a}+1K-h\iff u<v$

The above-given first condition is equivalent to $h\le \sqrt{a}-1/\sqrt{a}+1K$ since $a>1$. This constraint is the feasibility constraint for harvesting also already needed for ${\stackrel{⟶}{y}}_2$. The above-given second condition is equivalent to just to reverse the equality in the constraint (24), that is,$$\begin{array}{}\text{(26)}& h>\frac{a+1\sqrt{a}-2a}{a-1\sqrt{a}}K=\frac{\sqrt{a}-1}{\sqrt{a}+1}K\text{,}\end{array}$$so that ${\overline{y}}_1\ne {\overline{y}}_2$ is both feasible and locally asymptotically stable if and only if$$\begin{array}{}\text{(27)}& \frac{\sqrt{a}-1}{\sqrt{a}+1}K\ge h>\frac{\sqrt{a}-1}{\sqrt{a}+1}K\text{,}\end{array}$$which is a contradiction. Thus, ${\overline{y}}_1$ is unstable if feasible and distinct of ${\overline{y}}_2$.

If ${\overline{y}}_1={\overline{y}}_2$ is feasible, that is, the radicand of (5) is null so that $h=\sqrt{a}-1/\sqrt{a}+1K$, then the equilibrium point is given by ${\overline{y}}_1={\overline{y}}_2=K-h/2=K/\sqrt{a}+1$ which satisfies trivially the above given local asymptotic stability condition ${\overline{y}}_1={\overline{y}}_2\ge K/\sqrt{a}+1$ so that the confluent nonextinction equilibrium point ${\overline{y}}_1={\overline{y}}_2=K-h/2=K/\sqrt{a}+1$ resulting with $h=\sqrt{a}-1/\sqrt{a}+1K$ is locally asymptotically stable. Property (iv) has been proved.

Concerning Theorem 2(i), note that the denominator in the right-hand-side of (8) cannot be zero at any sample since the the value of the sequence ${h_n}_{n=0}^{\infty }$ is bounded by hypothesis.

Remark 1.

Note that, the admissible harvesting sequence of Theorem 2(iii) can be generated from (8) by generating ${h_n^I}_{n=0}^{\infty }$ as follows by fulfilling one of the stipulations 1–3 for each $n\in {\mathbf{Z}}_{0+}$:

(a) Through the stipulation 1 in the proof of Theorem 2: $h_n^I=\epsilon -K_n^{-1}{a}_n-1/a_n-{\sigma }_n\ge 0$; $\forall n\in {\mathbf{Z}}_{0+}$, where $\epsilon \in {\mathbf{R}}_{+}$ is chosen such that $K_n\ge {\epsilon }^{-1}$; $\forall n\in {\mathbf{Z}}_{0+}$ and the sequence ${{\sigma }_n}_{n=0}^{\infty }$ is generated subject to $0\le {\sigma }_n\le \epsilon -K_n^{-1}{a}_n-1/a_n$; $\forall n\in {\mathbf{Z}}_{0+}$. Note that, $h_n^I\ge 0$ and $h_n\ge 0$; $\forall n\in {\mathbf{Z}}_{0+}$.

(b) Through the stipulation 2 in the proof of Theorem 2: $h_n^I=1-a_n/a_n{K}_n+{\sigma }_n$, where $0\le {\sigma }_n\le a_n-1/a_{n}2K_n^{-1}-\epsilon$; $\forall n\in {\mathbf{Z}}_{0+}$ and $\epsilon \in {\mathbf{R}}_{+}$ is chosen such that $K_n<{\epsilon }^{-1}$. Note that, $h_n^I<0$ and $h_n<0$.

(c) Through the stipulation 3 in the proof of Theorem 2: $h_n^I=1-a_n/a_{n}1/K_n+\epsilon +{\sigma }_n$, where $0\le {\sigma }_n\le a_n-1/a_n\epsilon$; $\forall n\in {\mathbf{Z}}_{0+}$ and $\epsilon \in {\mathbf{R}}_{+}$. Note that, $h_n^I<0$ and $h_n<0$.

Remark 2.

The local asymptotic stability of the equilibrium points addressed in Theorem 2 (iv) relies to the cases of absence of harvesting in the steady state dynamics ($h=0$) or in the cases of stationary fishing/hunting $h>0$ or stationary repopulation actions $h<0$. Those cases correspond to constant values of the harvesting sequence in finite time or asymptotically. In the paper, the dynamics is globally stable if the population solution sequence is bounded for any given finite initial condition. This circumstance might be compatible with the event that some of the equilibrium points be locally unstable, stable, or critically stable if there are more than one equilibrium points. An equilibrium point is said to be globally asymptotically stable if it globally stable and all solution converges asymptotically to such a point for any given finite initial conditions.

In the third case, the largest positive equilibrium point is larger under negative stationary harvesting (having a meaning of stationary repopulation and/or immigration to the habitat from outside), than the equilibrium point $K$ arising in the absence of harvesting. In the second case, the global stability condition leads to the conclusion that the larger equilibrium point ${\overline{y}}_2$ is locally asymptotically stable and the smaller one ${\overline{y}}_1$ is not locally asymptotically stable unless they are coincident for a stationary harvesting effort $h=\sqrt{a}-1/\sqrt{a}+1K$.

The following result proves that, if the harvesting action sequence has a limit $h$, then a limit point of the solution cannot exceed the amount $K-h$.

Proposition 1.

If ${a_n}_0^{\infty }\subset 1,\overline{a}⟶a$, ${K_n}_0^{\infty }\subset 0,\overline{K}⟶K$, ${h_n}_{n=0}^{\infty }⟶h$ and ${y_n}_{n=0}^{\infty }⟶y$ then$$\begin{array}{c}\text{(28)}& h=0⟺y=0\vee y=Kh>0⟺y<K-h\\ h<0⟺y>K+h\text{.}\end{array}$$

Proof.

On gets from (1) that$$\begin{array}{}\text{(29)}& h_n=\frac{a_n{K}_n{y}_n}{K_n+a_n-1y_n}-y_{n+1}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

If ${y_n}_{n=0}^{\infty }⟶y$, one gets by taking limits in (29) as $n⟶\infty$ that$$\begin{array}{}\text{(30)}& h=\frac{aK}{K+a-1y}-1y=\frac{a-1K-y}{K+a-1y}y\text{,}\end{array}$$and, equivalently,$$\begin{array}{}\text{(31)}& hK=a-1K-y-hy\text{,}\text{(32)}& hK=-hK=a-1y-h-K\text{\hspace{0.17em}if\hspace{0.17em}}h<0\text{.}\end{array}$$

Then, the given properties follow directly from (31) and .

The following result establishes the boundedness of the solution sequence under bounded non-negative harvesting.

Proposition 2.

Assume that ${a_n}_0^{\infty }\subset 1,\overline{a}$, ${K_n}_0^{\infty }\subset 0,\overline{K}$, ${h_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}$ If ${h_n}_{n=0}^{\infty }\subset 0,a_n{K}_n{y}_n/K_n+a_n-1y_n$; $\forall n\in {\mathbf{Z}}_{0+}$ then ${y_n}_0^{\infty }\subset {\mathbf{R}}_{0+}$ is bounded.

Proof.

Assume on the contrary that ${y_n}_{n=0}^{\infty }⟶+\infty$.Then, from L´Hopital rule for quotients with numerator and denominator tending to infinity, since the harvesting sequence is non-negative,$$\begin{array}{}\text{(33)}& \underset{y_n⟶\infty n⟶\infty }{\lim }\sup y_{n+1}-\frac{a_n{K}_n{y}_n}{K_n+a_n-1y_n}=\underset{n⟶\infty }{\lim }\sup y_{n+1}-\frac{a_n{K}_n}{a_n-1}\le 0\text{,}\end{array}$$so that$$\begin{array}{}\text{(34)}& \underset{n⟶\infty }{\lim \sup }{y}_{n+1}=\underset{n⟶\infty }{\lim }{y}_{n+1}\le \underset{n⟶\infty }{\lim \inf }\frac{a_n{K}_n}{a_n-1}<+\infty \text{.}\end{array}$$

A contradiction for the sequence ${y_n}_{n=0}^{\infty }$ to diverge, which completes the proof.

3. Control Laws for Monitoring the Harvesting Action

The first part of this section is addressed to derive harvesting control laws based on Theorem 2(iv) guaranteeing the convergence to a prescribed equilibrium point $x^{\ast }$ under the assumption ${\sum }_{n=0}^{\infty }{a}_n-1/a_n=\infty$ on the intrinsic growth sequence ${a_n}_{n=0}^{\infty }$. Known criteria for absolute convergence of series or for convergence of series of non-negative elements to a prescribed limit can be used to calculate the harvesting control sequence ${h_n}_{n=0}^{\infty }$ based on the previous calculation of ${h_n^I}_{n=0}^{\infty }$, which reflects the harvesting effect in the inverse of the solution sequence, so as to satisfy Theorem 2(iv). The last part of the section proposes harvesting control laws which make the solution sequence of the population evolution to sample-to-sample, rather than asymptotically, behave according to a prescribed suitable pattern.

Now, rewrite the population solution sequence ${x_n=y_n^{-1}}_{n=0}^{\infty }$ as an equilibrium perturbation in the form ${{\tilde{x}}_n+x^{\ast }}_{n=0}^{\infty }$, where $x^{\ast }$ is the suited equilibrium point and ${\tilde{x}}_n=x_n-x^{\ast }$; $\forall n\in {\mathbf{Z}}_{0+}$. Thus, one can rewrite from (6) the one-step ahead evolution of the incremental sequence ${\tilde{x}}_n$ in the form of Lemma 1.2 (iii) of [24] as follows:$$\begin{array}{}\text{(35)}& {\tilde{x}}_{n+1}=1-b_n{\tilde{x}}_n-b_n{x}^{\ast }+{\nu }_n+h_n^I=1-b_n{\tilde{x}}_n+{\omega }_n+c_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$where for each $n\in {\mathbf{Z}}_{0+}$,$$\begin{array}{}\text{(36)}& b_n=1-{\mu }_n=1-a_n^{-1}=\frac{a_n-1}{a_n}\text{,}\text{(37)}& {\omega }_n={\nu }_n+{\beta }_n=\frac{a_n-1}{a_n{K}_n}+{\beta }_n\text{,}\text{(38)}& c_n=h_n^I-b_n{x}^{\ast }-{\beta }_n\text{,}\end{array}$$for any ${{\beta }_n}_0^{\infty }\subset \mathbf{R}$ chosen such that $h_n^I-b_n{x}^{\ast }\ge {\beta }_n\ge 1-a_n/a_n{K}_n$; $\forall n\in {\mathbf{Z}}_{0+}$ which guarantees that ${\omega }_n\ge 0$ and $c_n\ge 0$; $\forall n\in {\mathbf{Z}}_{0+}$. Note that, $b_n\in 0,1$; $\forall n\in {\mathbf{Z}}_{0+}$. Sufficient conditions for convergence ${{\tilde{x}}_n}_0^{\infty }\subset {\mathbf{R}}_{0+}⟶0⟺{x_n}_0^{\infty }⟶x^{\ast }$ are $\widetilde{{\mathrm{x}}_0}\ge 0$, and$$\begin{array}{}\text{(39)}& {\sum }_{n=0}^{\infty }{b}_n={\sum }_{n=0}^{\infty }\frac{a_n-1}{a_n}=\infty \text{,}\text{(40)}& {\omega }_n=o{b}_n\iff \underset{n⟶\infty }{\lim }\frac{{\beta }_n}{a_n-1}+\frac{1}{a_n{K}_n}=0\text{,}\text{(41)}& {\sum }_{n=0}^{\infty }{c}_n<\infty \text{.}\end{array}$$

The condition (40) is equivalent to ${\lim }_{n⟶\infty }{\beta }_n-1-a_n/a_n{K}_n=0$ and the condition (41) is guaranteed if$$\begin{array}{}\text{(42)}& {\sum }_{n=0}^{\infty }{h}_n^I+\frac{1-a_n}{a_n}{x}^{\ast }-K_n^{-1}<\infty \text{,}\end{array}$$since, from (38) and the condition $h_n^I-b_n{x}^{\ast }\ge {\beta }_n\ge 1-a_n/a_n{K}_n$; $\forall n\in {\mathbf{Z}}_{0+}$, one gets:$$\begin{array}{c}\text{(43)}& {\displaystyle \sum }_{n=0}^{\infty }{c}_n={\displaystyle \sum }_{n=0}^{\infty }{h}_n^I-b_n{x}^{\ast }-{\beta }_n\le {\displaystyle \sum }_{n=0}^{\infty }{h}_n^I-\frac{a_n-1}{a_n}{x}^{\ast }+\frac{a_n-1}{a_n{K}_n}={\displaystyle \sum }_{n=0}^{\infty }{h}_n^I+\frac{1-a_n}{a_n}{x}^{\ast }-K_n^{-1}<\infty \text{,}\end{array}$$which induces also the further necessary condition $\underset{n⟶\infty }{\lim }{h}_n^I+1-a_n/a_n{x}^{\ast }-K_n^{-1}=0$, since $c_n⟶0$ as $n⟶\infty$, guaranteed in turn if

  ${h_n^I}_{n=0}^{\infty }⟶0$ (equivalently, if ${h_n}_{n=0}^{\infty }⟶0$) and ${K_n^{-1}}_{n=0}^{\infty }⟶x^{\ast }$, equivalently, if ${K_n}_{n=0}^{\infty }⟶y^{\ast }$, or if

  ${h_n^I}_{n=0}^{\infty }⟶0\text{\hspace{0.17em}and\hspace{0.17em}}{a_n}_{n=0}^{\infty }⟶1$

3.1. Harvesting Control Law Based on d´ Alembert Convergence Criterion

a) ${\sum }_{n=0}^{\infty }{c}_n<\infty$ with $c_n\ge 0$; $\forall n\in {\mathbf{Z}}_{0+}$ is guaranteed under d´ Alembert convergence criterion in order to ${x_n}_{n=0}^{\infty }⟶x^{\ast }=1/y^{\ast }>0$, equivalently, ${y_n}_{n=0}^{\infty }⟶y^{\ast }>0$, if$$\begin{array}{}\text{(44)}& \frac{h_{n+1}^I-b_{n+1}{x}^{\ast }-{\beta }_{n+1}}{h_n^I-b_n{x}^{\ast }-{\beta }_n}={\gamma }_n\le \gamma <1\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$if $h_n^I\ne b_n{x}^{\ast }+{\beta }_n$ and $h_n^I\ge b_n{x}^{\ast }+{\beta }_n$; $\forall n\in {\mathbf{Z}}_{0+}$. Then, since $b_n=a_n-1/a_n$; $\forall n\in {\mathbf{Z}}_{0+}$, equation (44) leads to$$\begin{array}{c}\text{(45)}& b_{n+1}{x}^{\ast }+{\beta }_{n+1}=\frac{a_{n+1}-1}{a_{n+1}}{x}^{\ast }+{\beta }_{n+1}\le h_{n+1}^I=b_{n+1}{x}^{\ast }+{\beta }_{n+1}+{\gamma }_n{h}_n^I-b_n{x}^{\ast }-{\beta }_n\\ =\frac{a_{n+1}-1}{a_{n+1}}-{\gamma }_n\frac{a_n-1}{a_n}{x}^{\ast }+{\beta }_{n+1}+{\gamma }_n{h}_n^I-{\beta }_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$with ${\lim }_{n⟶\infty }{\beta }_n-1-a_n/a_n{K}_n=0$. One gets from (12), by using $x_n=y_n^{-1}$; $\forall n\in {\mathbf{Z}}_{0+}$, that$$\begin{array}{}\text{(46)}& y_{n+1}=x_{n+1}^{-1}=\frac{a_n{K}_n{x}_n^{-1}-h_n{K}_n+a_n-1x_n^{-1}}{K_n+a_n-1x_n^{-1}}=\frac{a_n{K}_n-h_n{K}_n{x}_n+a_n-1}{K_n{x}_n+a_n-1}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$so that, equivalently,$$\begin{array}{}\text{(47)}& x_{n+1}=\frac{K_n{x}_n+a_n-1}{a_n{K}_n-h_n{K}_n{x}_n+a_n-1}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$which equalized to (7) in Theorem 2 (i) leads to:$$\begin{array}{}\text{(48)}& h_n^I=h_n^I{h}_n=\frac{K_n{x}_n+a_n-1}{a_n{K}_n-h_n{K}_n{x}_n+a_n-1}+\frac{1}{a_n}\frac{1-a_n}{K_n}-x_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

Then, combining (45) and (48)$$\begin{array}{}\text{(49)}& h_{n+1}^I=\frac{a_{n+1}-1}{a_{n+1}}-{\gamma }_n\frac{a_n-1}{a_n}{x}^{\ast }+{\gamma }_n\frac{K_n{x}_n+a_n-1}{a_n{K}_n-h_n{K}_n{x}_n+a_n-1}+\frac{1}{a_n}\frac{1-a_n}{K_n}-x_n+{\beta }_{n+1}-{\gamma }_n{\beta }_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

Equivalently,$$\begin{array}{}\text{(50)}& h_{n+1}^I=\frac{a_{n+1}-1}{a_{n+1}}{x}^{\ast }+{\gamma }_n\frac{K_n{x}_n+a_n-1}{a_n{K}_n-h_n{K}_n{x}_n+a_n-1}+\frac{1}{a_n}1-a_n{K}_n^{-1}+x^{\ast }-x_n+{\beta }_{n+1}-{\gamma }_n{\beta }_n;\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$for some ${{\gamma }_n}_0^{\infty }\subset 0,\gamma \subset 0,1$ with ${\beta }_n\ge 1-a_n/a_n{K}_n$; $\forall n\in {\mathbf{Z}}_{0+}$. If the necessary limit condition ${\lim }_{n⟶\infty }{\beta }_n-1-a_n/a_n{K}_n=0$ for ${{\beta }_n}_{n=0}^{\infty }$ is forced for all values of the sequence, one has in the above equation that$$\begin{array}{}\text{(51)}& h_{n+1}^I=\frac{a_{n+1}-1}{a_{n+1}}{x}^{\ast }+{\gamma }_n\frac{K_n{x}_n+a_n-1}{a_n{K}_n-h_n{K}_n{x}_n+a_n-1}+\frac{1}{a_n}1-a_n{K}_n^{-1}+x^{\ast }-x_n+\frac{1-a_{n+1}}{a_{n+1}{K}_{n+1}}-{\gamma }_n\frac{1-a_n}{a_n{K}_n}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$so that for each $n\in {\mathbf{Z}}_{0+}$:

  Step 1: for a given scalar $\gamma \in 0,1$ and a prefixed equilibrium point $x^{\ast }>0$ such that ${x_n}_{n=0}^{\infty }⟶x^{\ast }$, so that, equivalently, ${y_n}_{n=0}^{\infty }⟶y^{\ast }=1/x^{\ast }$ and given at the $n$-th sampled time $K_n$, $K_{n+1}$, $a_n$, $a_{n+1}$, $x_n$, and $h_n$, one calculates $h_{n+1}^I$ for some gain ${\gamma }_n\le \gamma$ according to (51).

  Step 2: one calculates $y_{n+2}^{-1}=x_{n+2}=a_{n+1}^{-1}{x}_{n+1}+a_{n+1}-1/a_{n+1}{K}_{n+1}+h_{n+1}^I$ and then$$\begin{array}{c}\text{(52)}& h_{n+1}=\frac{a_{n+1}{K}_{n+1}{y}_{n+1}}{K_{n+1}+a_{n+1}-1y_{n+1}}-y_{n+2}=\frac{a_{n+1}{K}_{n+1}}{K_{n+1}{x}_{n+1}+a_{n+1}-1}-x_{n+2}^{-1}=\frac{a_{n+1}{K}_{n+1}}{K_{n+1}{x}_{n+1}+a_{n+1}-1}-\frac{a_{n+1}{K}_{n+1}}{K_{n+1}{x}_{n+1}+a_{n+1}-1-a_{n+1}{K}_{n+1}{h}_{n+1}^I}\text{.}\end{array}$$

Note that, the condition $h_n^I\ge b_n{x}^{\ast }+{\beta }_n$ with $h_n^I<0$ (and then $h_n<0$ implying hunting/fishing action) if $K_n^{-1}>x^{\ast }$ so that $K_n<y^{\ast }$. Conversely, $h_n^I<b_n{x}^{\ast }+{\beta }_n$ with $h_n^I\ge 0$ (and then $h_n\ge 0$ implying repopulation) if $K_n^{-1}\le x^{\ast }$ so that $K_n\ge y^{\ast }$.

3.2. Harvesting Control Law Based on Cauchy Root Test

Note that for a given $\gamma \in 0,1$ the sequence generated by$$\begin{array}{}\text{(53)}& h_n^I=b_n{x}^{\ast }+{\beta }_n+{\gamma }^n{\beta }_n\ge \frac{1-a_n}{a_n{K}_n}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$guarantees that ${\sum }_{n=0}^{\infty }{c}_n<\infty$ since $0\le c_n=h_n^I-b_n{x}^{\ast }-{\beta }_n={\gamma }_n^n\le {\gamma }^n<1$; $\forall n\in {\mathbf{Z}}_{0+}$. Furthermore, ${\omega }_n=o{b}_n$ is equivalent to ${\lim }_{n⟶\infty }{\beta }_n-1-a_n/a_n{K}_n=0$. It can be taken ${\beta }_n=1-a_n/a_n{K}_n+{\epsilon }_n$ with ${{\epsilon }_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}⟶0$. The control law is derived by Step 1 and Step 2 of Section 3.1 under the replacement of (51) by the subsequent equation which implements (51) by respecting the constraint on the sequence:$$\begin{array}{}\text{(54)}& h_n^I=b_n{x}^{\ast }+\frac{1-a_n}{a_n{K}_n}+{\epsilon }_n+{\gamma }^n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

3.3. Harvesting Control Law Based on Raabe Convergence Criterion

The condition ${\sum }_{n=0}^{\infty }{c}_n<\infty$ with $c_n=h_n^I-b_n{x}^{\ast }-{\beta }_n$; $\forall n\in {\mathbf{Z}}_{0+}$ is achieved from Raabe criterion if

$L_n=n1-c_{n+1}/c_n$; $\forall n\in {\mathbf{Z}}_{0+}$ with ${L_n}_{n=0}^{\infty }⟶L>1$. Furthermore, the condition ${\omega }_n=o{b}_n$ is achieved with ${\beta }_n=1-a_n/a_n{K}_n+{\epsilon }_n$; $\forall n\in {\mathbf{Z}}_{0+}$ with ${{\epsilon }_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}⟶0$. Then, the joint conditions $c_{n+1}=c_{n}1-L_n/n$ and ${\beta }_n=1-a_n/a_n{K}_n+{\epsilon }_n$; $\forall n\in {\mathbf{Z}}_{0+}$ with ${L_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}⟶L>1$ and ${{\epsilon }_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}⟶0$ are achieved for a given harvesting control $h_n$ at the $n$- sample if$$\begin{array}{c}\text{(55)}& h_{n+1}^I=b_{n+1}{x}^{\ast }+{\beta }_{n+1}+1-\frac{L_n}{n}{h}_n^I-b_n{x}^{\ast }-{\beta }_n=\frac{a_{n+1}-1}{a_{n+1}}-1-\frac{L_n}{n}\frac{a_n-1}{a_n}{x}^{\ast }+\frac{1-a_{n+1}}{a_{n+1}{K}_{n+1}}+{\epsilon }_{n+1}+1-\frac{L_n}{n}{h}_n^I-\frac{1-a_n}{a_n{K}_n}-{\epsilon }_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

Then, Step 1 and Step 2 of Section 3.1 are executed to get the harvesting control $h_{n+1}$ at the next $n+1$-sample.

3.4. Control Law Based on a Reference Model

It is assumed that ${y_n^{\ast }}_{n=1}^{\infty }$ is a suitable reference implicit model for the population. Note that, from (1) that $y_{n+1}=y_{n+1}^{\ast }$; $\forall n\in {\mathbf{Z}}_{0+}$ if $h_n=y_{n+1}^0-y_{n+1}^{\ast }$; $\forall n\in {\mathbf{Z}}_{0+}$, where the “a priori” population (that is, the harvesting free one) at the $n+1$-the sample is given by$$\begin{array}{}\text{(56)}& y_{n+1}^0=\frac{a_n{K}_n{y}_n}{K_n+a_n-1y_n}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

While the “a posteriori” one is $y_{n+1}$ which equalizes the targeted value by the reference model. If $y_{n+1}^0>y_{n+1}^{\ast }$ then the harvesting action $h_n$ implies a fishing or a hunting action while if $y_{n+1}^0<y_{n+1}^{\ast }$ then the harvesting consists in repopulation and if $y_{n+1}^0=y_{n+1}^{\ast }$ then no harvesting is performed. The following result is obvious from the above equations:

Proposition 3.

Assume that ${y_n^{\ast }}_{n=0}^{\infty }⟶y^{\ast }>0$. Then, the following properties hold:

(i) Then, ${y_{n+1}^0-h_n}_{n=0}^{\infty }⟶y^{\ast }$ and ${y_n}_{n=0}^{\infty }⟶y^{\ast }$.

(ii) Define ${{\lambda }_n}_{n=0}^{\infty }$ by ${\lambda }_n=a_n/1+a_n-1/K_n{y}_n$ so that $y_{n+1}^0={\lambda }_n{y}_n$ and define ${{\sigma }_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}$, with ${{\sigma }_n}_{n=0}^{\infty }⟶1$, such that $y_{n+1}=y_{n+1}^{\ast }={\sigma }_n{y}^{\ast }$. Then, $h_n=y_{n+1}^0-y_{n+1}^{\ast }={\lambda }_n{y}_n-{\sigma }_n{y}^{\ast }$ and ${h_n+1-{\lambda }_n{y}^{\ast }}_{n=0}^{\infty }⟶0$.

(iii) If ${y_n-K_n}_{n=0}^{\infty }⟶0$, equivalently, if ${K_n}_{n=0}^{\infty }⟶K^{\ast }=y^{\ast }$ (since ${y_n^{\ast }}_{n=0}^{\infty }⟶y^{\ast }$) then ${{\lambda }_n}_{n=0}^{\infty }⟶1$, ${y_n^0}_{n=0}^{\infty }⟶y^{\ast }$ and ${h_n}_{n=0}^{\infty }⟶0$.

The above-given result establishes that the population sequence convergence to a limit neither requires the convergence of its parameterizing sequences defining the intrinsic growth rates and the carrying capacities nor implies the convergence to zero of the harvesting control sequence (Proposition 1 (ii)).

4. Model Extensions by considering Allee and Overpopulation Effects

Note that, $y_n=h_n=0$ in (1) for some $n\in {\mathbf{Z}}_{0+}$ implies that $y_{n+1}=0$, that is, extinction in finite time in the absence of positive harvesting implies that the extinction remains afterwards for all time. However, the asymptotic extinction is not possible according to this model if $y_0>0$ in the presence of negative or null harvesting as it is obvious from the instability of the zero equilibrium point. Note also that ${y_n}_{n=0}^{\infty }$ cannot be strictly decreasing as addressed as follows:

Proposition 4.

Assume that ${h_n}_{n=0}^{\infty }$ is nonzero. If ${\lim \sup }_{n⟶\infty }{h}_n-a_n-1/K_n+a_n-1y_n{K}_n-y_n{y}_n<0$ then the asymptotic extinction is not possible. If ${h_n}_{n=0}^{\infty }\equiv 0$, then the asymptotic extinction is not possible under positive initial conditions.

Proof.

Note that, if for any $n\in {\mathbf{Z}}_{0+}$, $y_{n+1}<y_n$ then$$\begin{array}{}\text{(57)}& y_{n+1}-y_n=y_n\frac{a_n{K}_n}{K_n+a_n-1y_n}-1-h_n=y_n\frac{a_n-1}{K_n+a_n-1y_n}{K}_n-y_n-h_n<0\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

Thus, if ${\lim \sup }_{n⟶\infty }{h}_n-a_n-1/K_n+a_n-1y_n{K}_n-y_n{y}_n<0$ then the asymptotic extinction is not possible. Next, assume that $h_n=0$; $\forall n\in {\mathbf{Z}}_{0+}$. Then, $y_{n+1}<y_n$ for some $n\in {\mathbf{Z}}_{0+}$ implies that $y_n>K_n$. Take a positive real constant $\epsilon <{\inf }_{n\in {\mathbf{Z}}_{0+}}{K}_n$. Assume asymptotic extinction so that ${y_n}_0^{\infty }⟶0$ and assume also that $\max y_{n+1},0<y_n\le \epsilon$ for some $n\in {\mathbf{Z}}_{0+}$.Then, $y_n>\max y_{n+1},K_n>\max y_{n+1},\epsilon \ge \epsilon$ which is a contradiction to $y_n\le \epsilon$. As a result, if $y_n\le \epsilon$ then $y_{n+1}\ge y_n$ so that ${y_n}_0^{\infty }⟶0$ is not possible.

4.1. Consider Allee Effect for Small Number of Individuals

The well-known Allee effect establishes that under small population numbers, extinction is possible because of the difficulties for the individuals to meet members of the cohort. Now, we modify (1) with a density-dependent function $f:{\mathbf{R}}_{0+}⟶{\mathbf{R}}_{+}$ which penalizes the presence of small numbers of individuals as follows. It is assumed in the sequel that $h_n\equiv 0$ in (1) which is coherent with forbidding harvesting under very small amounts of individuals in the environment: It is still assumed that the intrinsic growth rate exceeds unity for all time.$$\begin{array}{}\text{(58)}& y_{n+1}=\frac{f{y}_n{a}_n{K}_n{y}_n}{K_n+a_n-1y_n}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

Theorem 3.

The following properties hold:

(i) If $f{y}_n<a_n^{-1}$ for all $n\in {\mathbf{Z}}_{0+}$ such that $y_n\ne 0$ which is guaranteed under the stronger sufficient condition ${\sup }_{n\in {\mathbf{Z}}_{0+}}f{y}_n<1/{\sup }_{n\in {\mathbf{Z}}_{0+}}{a}_n$ then ${y_n}_{n=0}^{\infty }$ is strictly decreasing so that it converges to zero.

(ii) Assume that $f{y}_n={\alpha }_n{y}_n^p+{\beta }_n$ where $p\in {\mathbf{R}}_{+}$, ${{\alpha }_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}$, ${{\beta }_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}$. If ${\beta }_n=0$ for $y_n=0$; ${\alpha }_n{\beta }_n=0$, ${\beta }_n<a_n^{-1}$; $\forall n\in {\mathbf{Z}}_{0+}$, and$$\begin{array}{}\text{(59)}& {\alpha }_n<\frac{a_n-1}{a_n{K}_n}{y}_n^{1-p}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$

  then ${y_n}_{n=0}^{\infty }$ is strictly decreasing so that it converges to zero. If (59) is replaced with ${\lim \sup }_{n⟶\infty }{\alpha }_n-a_n-1/a_n{K}_n{y}_n^{1-p}<0$ then ${y_n}_{n=0}^{\infty }$ is not necessarily strictly decreasing while still ${y_n}_{n=0}^{\infty }⟶0$.

Proof.

Note from (58) that $y_{n+1}<y_n$ if and only if $K_n+a_n-1y_n>f{y}_n{a}_n{K}_n$, and equivalently$$\begin{array}{}\text{(60)}& f{y}_n<\frac{1}{a_n}+\frac{a_n-1}{a_n{K}_n}{y}_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$which is guaranteed for all $n\in {\mathbf{Z}}_{0+}$ if $f{y}_n<a_n^{-1}$ such that $y_n\ne 0$; $\forall n\in {\mathbf{Z}}_{0+}$, if ${\sup }_{n\in {\mathbf{Z}}_{0+}}f{y}_n<1/{\sup }_{n\in {\mathbf{Z}}_{0+}}{a}_n$. Property (i) has been proved. On the other hand, the condition $y_{n+1}<y_n$; $\forall n\in {\mathbf{Z}}_{0+}$ implies that$$\begin{array}{}\text{(61)}& f{y}_n={\alpha }_n{y}_n^p+{\beta }_n<\frac{K_n+a_n-1y_n}{a_n{K}_n}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$which is guaranteed if$$\begin{array}{}\text{(62)}& {\alpha }_n<\frac{a_n-1}{a_n{K}_n}{y}_n^{1-p},{\beta }_n<\frac{1}{a_n}+\frac{a_n-1}{a_n{K}_n}-{\alpha }_n{y}_n^{p-1}{y}_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$the second above-given constraint being already guaranteed by the constraint ${\beta }_n<a_n^{-1}$; $\forall n\in {\mathbf{Z}}_{0+}$, implying also that $a_n-1/a_n{K}_n-{\alpha }_n{y}_n^{p-1}\ge 0$ for the given constraint on ${\alpha }_n$ in (62) and then ${y_n}_{n=0}^{\infty }$ is strictly decreasing so that it converges to zero. Note also that if ${\lim \sup }_{n⟶\infty }{\alpha }_n-a_n-1/a_n{K}_n{y}_n^{1-p}<0$ then ${y_n}_{n=0}^{\infty }⟶0$. Property (ii) has been proved.

Remark 3.

Note that, if $p=1$ with ${\beta }_n$ being zero or close to zero, then Theorem 3(iii) follows if ${\beta }_n<a_n^{-1}$ and ${\alpha }_n<a_n-1/a_n{K}_n$; $\forall n\in {\mathbf{Z}}_{0+}$. Note also that a coherent $f{y}_n$ under the above-given theorem to describe the Allee effect is $f{y}_n={\alpha }_n{y}_n^p+{\beta }_n$ with $0<p\le 1$ reflecting the difficulties for small numbers to find partners. Exponents $p>1$ might lead to large numbers of individuals, in mathematical terms, to the unbounded growing of the population. A parallel description for $p<0$ in $f{y}_n$ might be appropriate to describe competence within the cohort, for instance for food or refuge seeking, since it penalizes the presence of huge numbers of individuals in the cohort along the evolution process.

The penalty for high levels of populations mentioned in the last remark can also be described by introducing quadratic or higher terms in the denominator of the evolution equation as it is now addressed.

4.2. Incorporation of a Penalty Term for High Number of Individuals

In the following, we introduce a penalty term to deal with the presence of a high number of individuals. It is still assumed that the intrinsic growth rate sequence exceeds unity.

Proposition 5.

Consider the modified Beverton–Holt equation:$$\begin{array}{}\text{(63)}& y_{n+1}=\frac{a_n{K}_n{y}_n}{K_n+a_n-1y_n+c_n{y}_n^2}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$with ${c_n}_0^{\infty }\subset {\mathbf{R}}_{+}$. The following properties hold:

(i) Neither asymptotic extinction nor asymptotic population unboundedness are possible for any given finite initial condition $y_0>0$.

(ii) Assume that ${a_n}_{n=0}^{\infty }⟶a>1$, ${K_n}_{n=0}^{\infty }⟶K>0$ and ${c_n}_{n=0}^{\infty }⟶c>0$. Then, there is a unique positive equilibrium point $\overline{y}=a-1/2c\sqrt{1+4Kc/a-1}-1<\min K,\sqrt{Ka-1/c}$ for $c>0$, that is, it is smaller than the equilibrium point $\overline{y}=K$ for $c=0$ and it is a globally asymptotically stable attractor. The extinction equilibrium point $\overline{y}=0$ is unstable for $c\ge 0$.

Proof.

Rewrite the Beverton–Holt equation as$$\begin{array}{}\text{(64)}& y_{n+1}=\frac{a_n{K}_n{y}_n}{K_n+a_n-1+c_n{y}_n{y}_n}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

Thus, the population inverse is described by the following equation:$$\begin{array}{}\text{(65)}& x_{n+1}=\frac{K_n+c_n{y}_n^2}{a_n{K}_n}{x}_n+\frac{a_n-1}{a_n{K}_n}=\frac{K_n{x}_n+c_n{y}_n}{a_n{K}_n}+\frac{a_n-1}{a_n{K}_n}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

Since $y_n{x}_n=1$ and $y_n^2{x}_n=y_n$; $\forall n\in {\mathbf{Z}}_{0+}$. Then, if $y_n⟶\infty$ as $n⟶\infty$, $x_n⟶0$ as $n⟶\infty$ and, from (65), $x_{n+1}⟶\infty$ and $y_{n+1}⟶0$, and $y_{n+k}⟶0$; $\forall k\in {\mathbf{Z}}_{0+}$, a contradiction to the unboundedness of ${y_n}_{n=0}^{\infty }$. Thus, ${y_n}_{n=0}^{\infty }$ is bounded for any given finite $y_0>0$.On the other hand, if ${y_n}_{n=0}^{\infty }⟶0$, for any arbitrarily small $\epsilon \in {\mathbf{R}}_{+}$ there is some arbitrarily large $m=m\epsilon \in {\mathbf{Z}}_{0+}$ such that $y_{m+1}<y_m\le \epsilon$, equivalently, $a_m-1+c_m{y}_m{y}_m>a_m-1K_m$. But then $a_m-1+c\epsilon \epsilon >a_m-1K_m$ so that $\epsilon$ cannot be arbitrarily small. Therefore, asymptotic extinction is not possible and for any finite initial condition $y_0>0$, one has 0 < ${\lim \inf }_{n⟶\infty }{y}_n\le {\lim \sup }_{n⟶\infty }{y}_n<+\infty$. Property (i) has been proved.

To prove Property (ii), note that the limiting equation at an eventual equilibrium point is given by$$\begin{array}{}\text{(66)}& \overline{y}=\frac{aK\overline{y}}{K+a-1+c\overline{y}\overline{y}}\text{,}\end{array}$$which is satisfied by the extinction equilibrium point ${\overline{y}}_0=0$ and by $\overline{y}>0$ which satisfies the constraint $a-1K-\overline{y}=c{\overline{y}}^2$ from (66). If $c=0$ then $\overline{y}=K$, since $a>1$, and if $c>0$ then $\overline{y}<K$. Now, we prove that $\overline{y}$ is globally asymptotically stable. First, we prove that it is locally asymptotically stable and in a second step that it is globally asymptotically stable. Consider the linearized evolution of the inverse $x_n=\overline{x}+\delta x_n$ of $y_n$, where $\overline{x}={\overline{y}}^{-1}$ under a perturbation around the equilibrium point leading to $x_{n+1}=\overline{x}+\delta x_{n+1}$. One gets:$$\begin{array}{c}\text{(67)}& \overline{x}+\delta x_{n+1}=\frac{1}{a}\overline{x}+\delta x_n+\frac{a-1}{aK}\frac{\overline{y}+\delta y_n}{\overline{y}+\delta y_n}+\frac{c}{aK}\overline{y}+\delta y_n=\frac{1}{a}x+\delta x_n+\frac{a-1}{aK}+\frac{c}{aK}\frac{1}{\overline{x}}-{\overline{y}}^2\delta x_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$since $\delta y_n=\delta 1/x_n=-{1/x_n^2}_{x_n=\overline{x}}\delta x_n=-1/{\overline{x}}^2\delta x_n=-{\overline{y}}^2\delta x_n$; $\forall n\in {\mathbf{Z}}_{0+}$. From the above-given identity, the one-step ahead evolution of the linearized perturbation of the inverse of the solution is given by the following equation:$$\begin{array}{}\text{(68)}& \delta x_{n+1}=\frac{1}{a}1-\frac{c{\overline{y}}^2}{K}\delta x_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

Now, we prove that so that $c{\overline{y}}^2/K<1$. It has been proved before that the equilibrium point $y$ satisfies the constraints:$$\begin{array}{}\text{(69)}& a-1K-\overline{y}=c{\overline{y}}^2=aK+\overline{y}-a\overline{y}-K=aK-a-1\overline{y}-K<aK\text{.}\end{array}$$

since $a>1$ so that $0<c{\overline{y}}^2/K<1$. Furthermore, $0<1/a1-c{\overline{y}}^2/K<1$ since $a>1-c{\overline{y}}^2/K\iff a+c{\overline{y}}^2/K>1$ which holds since $a>1$. Since $1/a1-c{\overline{y}}^2/K<1$ then $\delta x_{n+1}/\delta x_n<1$ and $\delta y_{n+1}/\delta y_n={\overline{y}}^2\delta x_{n+1}/{\overline{y}}^2\delta x_n=\delta x_{n+1}/\delta x_n<1$; $\forall n\in {\mathbf{Z}}_{0+}$ since $y\ne 0$ and any eventual perturbation of the equilibrium point at any sample becomes reduced in size at the next sample. Therefore, the equilibrium point $\overline{y}$ is locally asymptotically stable.

We now prove that the extinction equilibrium point ${\overline{y}}_0=0$ is unstable. Take a small perturbation $\delta y=\delta y_n<\epsilon$, at any sample $n$, being less than $\epsilon >0$, of ${\overline{y}}_0=0$ which needs to be positive since a negative perturbed equilibrium point is unfeasible since it implies a negative population. Then, for $c>0$$$\begin{array}{}\text{(70)}& \frac{\delta y_{n+1}}{\delta y_n}>\frac{a}{1+a-1/K+c/K\epsilon \epsilon }>1\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$and the zero equilibrium point is unstable provided that$$\begin{array}{}\text{(71)}& \frac{c}{K}{\epsilon }^2+\frac{a-1}{K}\epsilon -a-1<0\text{,}\end{array}$$which holds if $\epsilon \in 0,a-1/2c\sqrt{1+4cK/a-1}-1$. If $c=0$ then, $\delta y_{n+1}/\delta y_n>1$ if $\epsilon \in 0,\sqrt{2}-1K$. Thus, the zero equilibrium point ${\overline{y}}_0=0$ associated with extinction is unstable. Since the nonzero equilibrium point $\overline{y}$ has already proved to be locally asymptotically stable then it is also globally asymptotically stable. The proof is complete.

5. Extinction Issues and Associated Stability Results in the Absence of Haversting

It is now investigated the extinction of the solution sequence ${y_n}_{n=0}^{\infty }$; $\forall n\in {\mathbf{Z}}_{0+}$ of the modified Beverton–Holt equation (63) if ${a_n}_{n=0}^{\infty }⟶a<1$, ${K_n}_{n=0}^{\infty }⟶K>0$ and ${c_n}_{n=0}^{\infty }⟶c\ge 0$. Basically, it is found that:

(i) If $c=0$ then the only locally stable equilibrium point is $\overline{y}=0$ (extinction). There is another equilibrium point ${\overline{y}}_1=K$ which is unstable. As a result, the solution sequence is globally stable (that is, bounded for all time and any given finite initial condition) but the population extinguishes asymptotically. On the other hand, remember from the former sections that if $a>1$, the extinction equilibrium point $\overline{y}=0$ is unstable while ${\overline{y}}_1=K$ is stable. As a result, the solution sequence is also globally stable and the population never extinguishes asymptotically if Allee´s effect is not considered.

(ii) If $0<c\le 1-a/4K$, there are, at least two, equilibrium points ${\overline{y}}_0=0$ and ${\overline{y}}_1>K$. If $c=1-a/4K$ there are no more equilibrium points than the extinction point and ${\overline{y}}_1$. If $c<1-a/4K$ then there is another equilibrium point ${\overline{y}}_2>{\overline{y}}_1$, The extinction equilibrium point ${\overline{y}}_0=0$ is locally asymptotically stable, ${\overline{y}}_1$ is unstable and ${\overline{y}}_2$ is locally stable or instable depending on the combined values of the parameterization triple $K>0,a\in 0,1,c\in 0,1-a/4K$. The solution sequence is still globally stable as before.

The equilibrium points satisfy $y_{n+1}=y_n=\overline{y}$; $\forall n\in {\mathbf{Z}}_{0+}$ so that they satisfy either $\overline{y}=0$ or $\overline{y}>0$. Note that if $\overline{y}>0$ is an equilibrium point then $0<c{\overline{y}}^2=1-a\overline{y}-K$ with $a<1$ and $\overline{y}>K$. The subsequent result addresses the existence and local asymptotic stability of the extinction and nonextinction equilibrium points and the global stability of the evolution sequence for small intrinsic growth rates such that $a\in 0,1$.

Theorem 4.

Assume that $a<1$. Then, the following properties hold:

(i) ${\overline{y}}_0=0$ is an equilibrium points. If $a<1$ and $c=1-a/4K$ then ${\overline{y}}_0=0$ and $\overline{y}=1-a/2c$ are the equilibrium points.

(ii) If $a<1$ then the nonzero feasibly equilibrium points are the positive real zeros of$$\begin{array}{}\text{(72)}& c{\overline{y}}^2-1-a\overline{y}+1-aK=0\text{,}\end{array}$$

  that is,$$\begin{array}{}\text{(73)}& {\overline{y}}_{\mathrm{1,2}}=\frac{1-a\pm \sqrt{{1-a}^2-4cK1-a}}{2c}\text{,}\end{array}$$

  so that, if $c>1-a/4K$, then ${\overline{y}}_0=0$ is an equilibrium point and there is no nonzero equilibrium point. If $0<c<1-a/4K$ then $\overline{y}=0$ and also there are two distinct positive equilibrium points, namely,$$\begin{array}{c}\text{(74)}& {\overline{y}}_1=\frac{1-a-\sqrt{{1-a}^2-4cK1-a}}{2c}>K\text{;}{\overline{y}}_2=\frac{1-a+\sqrt{{1-a}^2-4cK1-a}}{2c}<\frac{1-a}{c}\text{,}\end{array}$$

  which become identical if $c=1-a/4K$.

(iii) The inverse sequence ${x_n}_{n=0}^{\infty }$ of the evolution sequence ${y_n}_{n=0}^{\infty }$ of the limiting Beverton-Holt equation is given by the following equation:$$\begin{array}{}\text{(75)}& x_{n+1}=\frac{1}{a}{x}_n-\frac{1-a}{Ka}+\frac{c}{Ka{x}_n}={\frac{1}{a}}^{n+1}{x}_0+{\displaystyle \sum }_{i=0}^n{\frac{1}{a}}^i\frac{c}{Ka{x}_{n-i}}-\frac{1-a}{aK}\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

  The sequence ${y_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}$ is ultimately bounded for any finite initial condition $y_0>0$ and it can converge to some limit $0\le \overline{y}<1-a/c$ (already proved in (ii) since $0<{\overline{y}}_1\le {\overline{y}}_2<1-a/c$).

(iv) The null equilibrium point is locally asymptotically stable. In particular,

(1) ${y_n}_{n=0}^{\infty }\subset 0,y_0⟶0$ and it is strictly decreasing if $c=0$ and $0\le y_0<K$

(2) If $c>1-a/4K$ then ${y_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}⟶0$ is bounded and strictly decreasing for any finite $y_0\ge 0$ and ${y_n}_{n=0}^{\infty }0,y_0⟶0$ and it is strictly decreasing if $c\le 1-a/4K$ and $y_0\in 0,{\overline{y}}_1$, since ${\overline{y}}_1=1-a-\sqrt{{1-a}^2-4cK1-a}/2c$

(v) If$0<c\le 1-a/4K$ then the evolution sequence is globally stable for any given finite initial condition. If $c>1-a/4K$ then the evolution sequence is globally asymptotically stable for any given finite initial condition and the evolution sequence converges asymptotically towards the species extinction.

If $0\le c\le 1-a/4K$, then the extinction point is locally stable and the equilibrium point ${\overline{y}}_1$ is unstable. The equilibrium point ${\overline{y}}_2$ is unstable if ${\overline{y}}_2\in {\overline{y}}_1,2K\cup 2K/1-a,+\infty$ and it is locally asymptotically stable if ${\overline{y}}_2\in 2-a/1-aK,2K/1-a$.

Define the parameters ${\lambda }_c=c/a>0$ and ${\lambda }_K=K/a>0$ for $c>0$. Then, the above-given local asymptotic stability condition of ${\overline{y}}_2$ is that for given ${\lambda }_c$ and ${\lambda }_K$, $a\in 0,1$ satisfies the subsequent constraint:$$\begin{array}{}\text{(76)}& {2{\lambda }_K{\lambda }_c{a}^2+a-1}^2+\frac{{\lambda }_c{\lambda }_K}{1-a}{a}^2<\frac{1}{4}<{\frac{2{\lambda }_K{\lambda }_c{a}^2}{1-a}+a-1}^2+\frac{{\lambda }_c{\lambda }_K}{1-a}{a}^2\text{.}\end{array}$$

Proof.

Properties (i)-(ii) are direct by calculating the real positive zeros $y_{\mathrm{1,2}}$, equation (74), which satisfy the subsequent equation of the stationary modified Beverton–Holt equation (63), with $a<1$ and $c>0$, such that ${y_n}_0^{\infty }\subset {\mathbf{R}}_{0+}$ with $y_{n+1}/y_n<1$ is satisfied:$$\begin{array}{}\text{(77)}& c{\overline{y}}^2-1-a\overline{y}+1-aK=0\text{,}\end{array}$$requiring that $c\le 1-a/4K$ and implying also, if $\overline{y}\ne 0$, that $0<c{\overline{y}}^2=1-a\overline{y}-K$ so that $1-a/c\ge {\overline{y}}_2\ge {\overline{y}}_1>K$ if $c\le 1-a/4K$ (${\overline{y}}_2>{\overline{y}}_1>K$ if and only if $c<1-a/4K$). If $c>1-a/4K$ then the zeros are complex conjugate with nonzero imaginary parts so that ${\overline{y}}_0=0$ (extinction) is the unique equilibrium point. Properties (i) and (ii) have been proved.

The evolution of the inverse sequence ${x_n=y_n^{-1}}_{n=0}^{\infty }$ of the population evolution sequence ${y_n}_{n=0}^{\infty }$ of the stationary modified Beverton–Holt equation is given by the following equation:$$\begin{array}{}\text{(78)}& x_{n+1}=\frac{1}{a}{x}_n-\frac{1-a}{Ka}+\frac{c}{Ka}{y}_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$so that one gets the solution sequence (75) which is got via direct recursive calculation from initial conditions. The subsequent expression is equivalent to (75):$$\begin{array}{}\text{(79)}& x_{n+1}+\frac{1-a}{aK}{{\sum }_{i=0}^n\frac{1}{a}}^i={\frac{1}{a}}^{n+1}{x}_0+{{\sum }_{i=0}^n\frac{1}{a}}^i\frac{c}{Ka{x}_{n-i}}\text{;}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

Assume that ${x_n}_{n=0}^{\infty }⟶\overline{x}\in 0,+\infty$. Then, since $a<1$, one gets in the above equation for any given finite initial condition $x_0=1/y_0>0$:$$\begin{array}{}\text{(80)}& \underset{n⟶\infty }{\lim }\frac{1-a}{aK}{{\sum }_{i=0}^n\frac{1}{a}}^i=x_0\underset{n⟶\infty }{\lim }{\frac{1}{a}}^{n+1}=\underset{n⟶\infty }{\lim }{{\sum }_{i=0}^n\frac{1}{a}}^i\frac{c}{Ka{x}_{n-i}}=+\infty \text{,}\end{array}$$so that$$\begin{array}{c}\text{(81)}& \overline{x}=\underset{n,m⟶\infty }{\lim }{x}_{n+m+1}=\underset{n,m⟶\infty }{\lim }{x}_m\underset{n⟶\infty }{\lim }{\frac{1}{a}}^{n+1}+\underset{m,n⟶\infty }{\lim }{{\sum }_{i=m}^{m+n}\frac{1}{a}}^i\frac{c}{Ka{x}_{n+m-i}}-\frac{1-a}{aK}=x_0\underset{n+m⟶\infty }{\lim }{\frac{1}{a}}^{n+m+1}+\underset{m,n⟶\infty }{\lim }{{\sum }_{i=m}^{m+n}\frac{1}{a}}^i\frac{c}{Ka\overline{x}}-\frac{1-a}{aK}\\ =+\infty +\underset{m,n⟶\infty }{\lim }{{\sum }_{i=m}^{m+n}\frac{1}{a}}^i\frac{c}{Ka\overline{x}}-\frac{1-a}{aK}\text{.}\end{array}$$

Consider the subsequent claims:

(a) $\overline{x}=0$ then $\overline{y}=1/\overline{x}=+\infty$. This claim is impossible since then $+\infty =+\infty +\infty$ is a contradiction. Then, the sequence ${y_n}_{n=0}^{\infty }$ is ultimately bounded, that is, ${\lim \sup }_{n⟶\infty }{y}_n<+\infty$, for any given finite $y_0>0$.

(b) $\overline{x}\in 0,+\infty$ then $y=1/\overline{x}\in 0,+\infty$. This claim is possible with $\overline{x}>c/1-a$, so that $\overline{y}<1-a/c$, (including the case $\overline{x}=+\infty$, $\overline{y}=1/\overline{x}=0$) since then $0=+\infty -\infty$ is not a contradiction. It is not possible with $\overline{x}=c/1-a$, so that $\overline{y}=1-a/c$ since then $\overline{x}=c/1-a=+\infty$ is a contradiction, and it is not possible either with $\overline{x}<c/1-a$, so that $\overline{y}>1-a/c$, since then $0=+\infty +\infty$ is also a contradiction. Property (iii) has been proved.

To prove Property (iv), first take $c=0$ with $a<1$. Then,$$\begin{array}{}\text{(82)}& \frac{y_{n+1}}{y_n}=\frac{a}{1-1-a/K{y}_n}<1\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$which implies that $y_n<K$ with the extra constraint that $y_n<K/1-a$ in order that $0\le y_{n+1}<+\infty$. Both constraints together reduce to the first one since $y_n<\min K,K/1-a=K$. Thus, if $y_0<K$ then ${y_n}_{n=0}^{\infty }\subset 0,K⟶0$ and it is strictly decreasing. Thus, the extinction equilibrium point is locally asymptotically stable, as claimed, if $c=0$ and $a<1$. Now, assume that $c>0$. Then,$$\begin{array}{}\text{(83)}& \frac{y_{n+1}}{y_n}=\frac{Ka}{K-1-a{y}_n+c{y}_n^2}<1\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$with a positive denominator to ensure that ${y_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}⟶0$ is bounded and strictly decreasing. Both constraints jointly hold if $K-1-a{y}_n+c{y}_n^2>Ka$, equivalently,$$\begin{array}{}\text{(84)}& 1-aK-y_n+c{y}_n^2>0\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{,}\end{array}$$which holds unconditionally for any finite $y_0\ge 0$ if $c>1-a/4K$, since in that case ${\overline{y}}_1$ and ${\overline{y}}_2$ are unfeasible since they are complex conjugate with nonzero imaginary parts so that the unique equilibrium point is $\overline{y}=0$; and for $y_0\in 0,1-a-\sqrt{{1-a}^2-4cK1-a}/2c$ if $c\le 1-a/4K$. Property (iv) has been proved.

To prove Property (v), first note that that ${\lim y_1}_{y_0⟶\infty }=0$ if $c>0$. Thus, for any large initial condition $y_0$, it follows that $y_1$ is arbitrarily close to the extinction equilibrium point ${\overline{y}}_0=0$, which is locally asymptotically stable so that ${y_n}_{n=0}^{\infty }⟶0$ for sufficiently large initial condition $y_0$ and ${y_n}_{n=0}^{\infty }$ is bounded for any given finite initial condition. Therefore, any evolution sequence is globally stable and it converges asymptotically to zero if $c>1-a/4K$ since then the unique locally asymptotically stable equilibrium point is the extinction so that it is also a global attractor. Now, the local stability/instability properties of ${\overline{y}}_1$ and ${\overline{y}}_2$ of Property (v) are proved. Note from Property (iv) that:

(a) If $c>1-a/4K$ then ${y_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}⟶0$ is bounded and strictly decreasing for any finite initial condition $y_0\ge 0$. Also, both ${\overline{y}}_1$ and ${\overline{y}}_2$ are unfeasible;

(b) If $0<c\le 1-a/4K$ then ${y_n}_{n=0}^{\infty }\subset {\mathbf{R}}_{0+}⟶0$ if $y_0\in 0,{\overline{y}}_1$ so that $\overline{y}=0$ is locally asymptotically stable and ${\overline{y}}_1$ is unstable since any solution from any initial condition $y_0<y_1$, being arbitrarily close to $y_1$ converges asymptotically to $\overline{y}=0$.

It remains to prove that ${\overline{y}}_2$ is conditionally stable. Note that, for a first-order perturbation of the equilibrium point ${\overline{y}}_2$ at any $n$-th sample, one has that $x_n={\overline{x}}_2+\delta x_n$, $x_{n+1}={\overline{x}}_2+\delta x_{n+1}$, $\delta y_n=\delta 1/x_n=-1/{\overline{x}}_2^2\delta x_n=-{\overline{y}}_2^2\delta x_n$ with ${\overline{x}}_2=1/{\overline{y}}_2$ so that$$\begin{array}{}\text{(85)}& \delta x_{n+1}=\frac{1}{a}\delta x_n-\frac{c}{Ka{\overline{x}}_2^2}\delta x_n=\frac{1}{a}1-\frac{c}{K}{\overline{y}}_2^2\delta x_n\text{;\hspace{0.17em}}\forall n\in {\mathbf{Z}}_{0+}\text{.}\end{array}$$

Then, ${\overline{y}}_2$ is unstable if ${\delta x_{n+1}/\delta x_n}_{y={\overline{y}}_2}=1/a1-c/K{\overline{y}}_2^2\ge 1$; $\forall n\in {\mathbf{Z}}_{0+}$, equivalently if,$$\begin{array}{}\text{(86)}& 1-\frac{c}{K}{\overline{y}}_2^2=2-a-\frac{1-a}{K}{\overline{y}}_2\ge a\text{,}\end{array}$$which is satisfied under the two subsequent conditions:

(a) $2-a-1-a/K{\overline{y}}_2\ge a$, which is equivalent to ${\overline{y}}_2\le 2K$, which together with the former constraints ${\overline{y}}_2>K$ and $K<1-a/c$, yields $K<{\overline{y}}_1\le {\overline{y}}_2\le 2K<21-a/c$;

(b) $1-a/K{\overline{y}}_2-2-a\ge a$ which is equivalent to ${\overline{y}}_2\ge 2K/1-a$.

   As a result, ${\overline{\mathrm{y}}}_2$ is instable if ${\overline{\mathrm{y}}}_2\in {\overline{\mathrm{y}}}_1,2\mathrm{K}\cup 2\mathrm{K}/1-\mathrm{a},+\infty$.

   In the same way, in order for ${\overline{y}}_2$ to be locally asymptotically stable, $2-a-1-a/K{\overline{y}}_2<a$ which is fulfilled under two conditions, namely,

(c) $0\le 2-a-1-a/K{\overline{y}}_2<a$ which holds if and only if $\overline{{\mathrm{y}}_2}\in 2\mathrm{K},2-\mathrm{a}/1-\mathrm{a}\mathrm{K}$;

(d) $0\le 1-a/K{\overline{y}}_2-2-a<a$ which holds if and only if ${\overline{y}}_2\in 2-aK/1-a,2K/1-a$

which are equivalently combined into the local stability condition ${\overline{y}}_2\in 2K,2K/1-a$. Taking into account that ${\overline{y}}_2=1-a/c1+\sqrt{1-4cK/1-a}$, the condition ${\overline{y}}_2\in 2-aK/1-a,2K/1-a$ becomes equivalent to$$\begin{array}{}\text{(87)}& {\frac{2c2-a}{{1-a}^2}-1}^2\le 1-\frac{4cK}{1-a}<{\frac{4Kc}{{1-a}^2}-1}^2\text{,}\end{array}$$subject to $0\le c<1-a/4K$ and define parameters ${\lambda }_c=c/a>0$ (note that ${\lambda }_c=0$ does not need to be considered since then ${\overline{y}}_2={\overline{y}}_1=K$ is stable) and ${\lambda }_K=K/a>0$. Then, the constraint $c<1-a/4K$ takes the form $4{\lambda }_c{\lambda }_K{a}^2+a-1<0$ and the above local stability constraint (87) of ${\overline{y}}_2$ takes the form (76).

Remark 4.

Note from Theorem 4 (ii) that any nonzero equilibrium point is less than the carrying capacity if $c>0$ contrarily to the case when $c=0$ where $\overline{y}=K$ is an equilibrium point.

Remark 5.

Note that, contrarily to Theorem 4(iv), if $a\ge 1$ and $c=0$, then, in the absence of harvesting, ${y_n}_{n=0}^{\infty }$ cannot be strictly decreasing converging to zero since, for the limit solution sequence to be strictly decreasing, it is necessary that

$y_{n+1}/y_n=a/1+a-1/K{y}_n<1$ what implies that $y_n>K$ so that ${y_n}_{n=0}^{\infty }⟶0$ is impossible for any given positive finite initial condition if $a_n\ge 1$; $\forall n\in {\mathbf{Z}}_{0+}$. Thus, the extinction equilibrium point is unstable if $a_n\ge 1$; $\forall n\in {\mathbf{Z}}_{0+}$.

6. Numerical Simulations

6.1. Example 1

This example illustrates the results of Theorem 1. The sequences ${a_n}_0^{\infty }$, ${K_n}_0^{\infty }$ and ${z_n}_0^{\infty }$ are respectively generated by means of the following difference equationsc:$$\begin{array}{c}\text{(88)}& a_{n+1}={\epsilon }_1{a}_n+{\rho }_1\text{,}{K}_{n+1}={\epsilon }_2{K}_n+{\rho }_2\text{,}\\ z_{n+1}={\epsilon }_3{z}_n+{\rho }_3\text{,}\end{array}$$with the following values for the parameters:$$\begin{array}{}\text{(89)}& {\epsilon }_1=0.9\text{,\hspace{0.17em}}{\rho }_1=0.4\text{,\hspace{0.17em}}{\epsilon }_2=0.8\text{,\hspace{0.17em}}{\rho }_2=200\text{,\hspace{0.17em}}{\epsilon }_3=0.75\text{\hspace{0.17em}and\hspace{0.17em}}{\rho }_3=100\text{,}\end{array}$$and the following initial conditions:$$\begin{array}{}\text{(90)}& a_0=1.5\text{,\hspace{0.17em}}{K}_0=500\text{\hspace{0.17em}and\hspace{0.17em}}{z}_0=10\text{,}\end{array}$$

In this way the conditions of Theorem 1 about the sequences ${a_n}_0^{\infty }$ and ${K_n}_0^{\infty }$ are fulfilled since $a={\lim }_{n⟶\infty }{a}_n={\rho }_1/1-{\epsilon }_1=4$ and $K={\lim }_{n⟶\infty }{K}_n={\rho }_2/1-{\epsilon }_2=1000$.

In a first simulation, the harvesting sequence ${h_n}_0^{\infty }$ with $h_n=5\ast {0.7}^n$ is considered so that $h={\lim }_{n⟶\infty }{h}_n=0$ and then the conditions of Theorem 1(i) are satisfied. Figure 1 shows the evolution of the species population $y_n$ and that of the environment carrying capacity $K_n$ if the population is initially $y_0=200$. One can see that $\overline{y}={\lim }_{n⟶\infty }{y}_n={\lim }_{n⟶\infty }{K}_n=K=1000$ as Theorem 1(i) points out.

[figure(s) omitted; refer to PDF]

On the other hand, Figure 2 displays the evolution of the species population and that of the environment carrying capacity if the population is initially $y_0=10$. Again, one can see that $\overline{y}={\lim }_{n⟶\infty }{y}_n={\lim }_{n⟶\infty }{K}_n=K=1000$ although the initial population is close to the equilibrium point $\overline{y}=0$. The results displayed in Figures 1 and 2 illustrate the fact that the equilibrium point $\overline{y}=0$ is globally asymptotically unstable while the equilibrium point $\overline{y}=K$ is globally asymptotically stable.

[figure(s) omitted; refer to PDF]

In a second simulation the same values of (89) and (90) are maintained but the sequence ${h_n}_0^{\infty }$ converges to a nonzero value $h$. Concretely, the time evolution of ${h_n}_0^{\infty }$ is displayed in Figure 3 while that of $y_n$ and $K_n$ if the population is initially $y_0=145$ is shown in Figure 4.

[figure(s) omitted; refer to PDF]

In this example the fact that $0\ne h={\lim }_{n⟶\infty }{h}_n<{\lim }_{n⟶\infty }{K}_n=K$ is observed so that the conditions of Theorem 1(ii) and Theorem 1(iii) are satisfied. In fact, the Beverton–Holt equation (1) possesses two equilibrium points given by (1). Concretely, such equilibrium points are ${\overline{y}}_1=142.855$ and ${\overline{y}}_2=600$. One can see that ${\lim }_{n⟶\infty }{y}_n={\overline{y}}_2=600$ although the initial condition $y_0=145$ is close to the equilibrium point ${\overline{y}}_1$. Such a fact illustrates that ${\overline{y}}_1$ is unstable while ${\overline{y}}_2$ is globally asymptotically stable in the sense that all solutions generated by finite initial conditions converge to such an equilibrium.

6.2. Example 2

This example illustrates the results of Theorem 3 related with the Allee effect for small number of individuals in the species. The intrinsic growth rate and carrying capacity sequences, ${a_n}_0^{\infty }$ and ${K_n}_0^{\infty }$, are given by (88) with the same values for the parameters ${\epsilon }_1$, ${\epsilon }_2$, ${\rho }_1$, and ${\rho }_2$ than those pointed out in (89). Moreover, the harvesting sequence ${h_n}_0^{\infty }$ is zero for all $n\epsilon {\mathbb{Z}}_{0+}$ and the function $f{y}_n$ appearing in the modified Beverton-Holt equation (58) is given by the following equation:$$\begin{array}{}\text{(91)}& f{y}_n=\frac{y_n+{\alpha }_1}{a_n{y}_n+{\alpha }_2}\text{,}\end{array}$$with ${\alpha }_1=0.1$ and ${\alpha }_2=1$. In this way, $f{y}_n<a_n^{-1}$ for all $n\epsilon {\mathbb{Z}}_{0+}$ and the conditions of Theorem 3(i) are fulfilled. Figure 5 shows the evolution of the species population $y_n$ if the population is initially $y_0=15$. Figure 6 displays the evolution of the function $f{y}_n$ and the inverse of the sequence ${a_n}_0^{\infty }$. In Figure 5, one can see that the species population converges to the extinction as Theorem 3(i) establishes since $f{y}_n<a_n^{-1}$ for all $n\epsilon {\mathbb{Z}}_{0+}$ as it is shown in Figure 6.

[figure(s) omitted; refer to PDF]

Now, the function $f{y}_n$ appearing in the modified Beverton–Holt equation (58) is given by the following equation:$$\begin{array}{}\text{(92)}& f{y}_n={\alpha }_n{y}_n^p+{\beta }_n\text{,}\end{array}$$with ${\alpha }_n=0.7a_n-1/a_n{K}_n{y}_n^{1-p}$ and ${\beta }_n=0.8a_n^{-1}$ for all $n\epsilon {\mathbb{Z}}_{0+}$. In this way, the conditions of Theorem 3(ii) are fulfilled. Figure 7 shows the evolution of the species population $y_n$ if the population is initially $y_0=15$. One can see that the species population converges to the extinction as Theorem 3(ii) establishes.

[figure(s) omitted; refer to PDF]

6.3. Example 3

The following two examples illustrate the results of Proposition 5 about the modified Beverton–Holt equation (63).

(i) The sequences ${a_n}_0^{\infty }$, ${K_n}_0^{\infty }$ and ${c_n}_0^{\infty }$ are respectively defined as$$\begin{array}{}\text{(93)}& a_n=a_{0}1+0.5\sin 2\pi f_{a}n\text{,}\text{(94)}& K_n=K_{0}1+0.02\sin 2\pi f_{K}n\text{,}\text{(95)}& c_n=c_{0}1+0.3\sin 2\pi f_{c}n\text{,}\end{array}$$

  to illustrated the result (i) of such a proposition with $a_0=1.8$, $f_a=0.01$, $K_0=500$, $f_K=0.02$, $c_0=0.03$ and $f_c=0.03$. Figure 8 shows the evolution of the species population $y_n$ if the population is initially $y_0=100$. One can see that the species population neither extinguishes nor increases in an unboundedness way as Proposition 5 (i) establishes.

(ii) The sequences ${a_n}_0^{\infty }$, ${K_n}_0^{\infty }$, and ${c_n}_0^{\infty }$ are, respectively, generated by means of the following difference equations:$$\begin{array}{c}\text{(96)}& a_{n+1}={\epsilon }_1{a}_n+{\rho }_1\text{,}{K}_{n+1}={\epsilon }_2{K}_n+{\rho }_2\text{,}\\ c_{n+1}={\epsilon }_3{c}_n+{\rho }_3\text{,}\end{array}$$

  with the following values for the parameters:$$\begin{array}{}\text{(97)}& {\epsilon }_1=09\text{,\hspace{0.17em}}{\rho }_1=0.4\text{,\hspace{0.17em}}{\epsilon }_2=0.8\text{,\hspace{0.17em}}{\rho }_2=200\text{,\hspace{0.17em}}{\epsilon }_3=0.75\text{\hspace{0.17em}and\hspace{0.17em}}{\rho }_3=0.005\text{,}\end{array}$$

  and the following initial conditions:$$\begin{array}{}\text{(98)}& a_0=1.5\text{,\hspace{0.17em}}{K}_0=500\text{\hspace{0.17em}and\hspace{0.17em}}{c}_0=0.01\text{.}\end{array}$$

In this way the conditions of Proposition 5(ii) about the sequences ${a_n}_0^{\infty }$, ${K_n}_0^{\infty }$ and ${c_n}_0^{\infty }$ are fulfilled since $a={\lim }_{n⟶\infty }{a}_n={\rho }_1/1-{\epsilon }_1=4$, $K={\lim }_{n⟶\infty }{K}_n={\rho }_2/1-{\epsilon }_2=1000$ and $c={\lim }_{n⟶\infty }{c}_n={\rho }_3/1-{\epsilon }_3=0.02$. Figure 9 shows the evolution of the species population $y_n$ and that of the environment carrying capacity $K_n$ if the population is initially $y_0=200$. One can see that $\overline{y}={\lim }_{n⟶\infty }{y}_n=a-1/2c\sqrt{1+4Kc/a-1}-1\approx 319<\min K,\sqrt{Ka-1/c}=\min 1000,\sqrt{Ka-1/c}\approx 387$ as Proposition 1(ii) points out.

[figure(s) omitted; refer to PDF]

On the other hand, Figure 10 displays the evolution of the species population and that of the environment carrying capacity if the population is initially $y_0=10$. Again, one can see that $\overline{y}={\lim }_{n⟶\infty }{y}_n\approx 319$ although the initial population is close to the equilibrium point $\overline{y}=0$. The results displayed in Figures 9 and 10 illustrate the fact that the equilibrium point $\overline{y}=0$ is globally asymptotically unstable while the equilibrium point $\overline{y}=K$ is globally asymptotically stable for $c>0$.

[figure(s) omitted; refer to PDF]

6.4. Example 4

The following examples illustrate the results of Section 5 about the modified Beverton–Holt equation (63). The sequences ${a_n}_0^{\infty }$ and ${K_n}_0^{\infty }$ are, respectively, generated by means of the following difference equations:$$\begin{array}{c}\text{(99)}& a_{n+1}={\epsilon }_1{a}_n+{\rho }_1\text{,}{K}_{n+1}={\epsilon }_2{K}_n+{\rho }_2\text{,}\end{array}$$with the following values for the parameters:$$\begin{array}{}\text{(100)}& {\epsilon }_1=0.8\text{,\hspace{0.17em}}{\rho }_1=0.1\text{,\hspace{0.17em}}{\epsilon }_2=0.9\text{\hspace{0.17em}and\hspace{0.17em}}{\rho }_2=50\text{,}\end{array}$$and the following initial conditions:$$\begin{array}{}\text{(101)}& a_0=1.5\text{\hspace{0.17em}and\hspace{0.17em}}{K}_0=300\text{.}\end{array}$$

In this way the conditions of Section 5 about the sequences ${a_n}_0^{\infty }$ and ${K_n}_0^{\infty }$ are fulfilled since $a={\lim }_{n⟶\infty }{a}_n={\rho }_1/1-{\epsilon }_1=0.5$ and $K={\lim }_{n⟶\infty }{K}_n={\rho }_2/1-{\epsilon }_2=500$. Several choices for the sequence ${c_n}_0^{\infty }$ are considered to illustrate the results of Section 5:

(i) In the first case ${c_n}_0^{\infty }$ is given by the difference equation:$$\begin{array}{}\text{(102)}& c_{n+1}={\epsilon }_3{c}_n+{\rho }_3\text{,}\end{array}$$

  with the values for the parameters ${\epsilon }_3=0.75$ and ${\rho }_3=0$ and the initial condition $c_0=0.01$. In this way, $c={\lim }_{n⟶\infty }{c}_n={\rho }_3/1-{\epsilon }_3=0$ and then, the conditions that $a<1$, $K>0$ and $c=0$ are fulfilled so that the modified Beverton–Holt equation has two equilibrium points, namely, the stable point ${\overline{y}}_1=0$ (extinction) and the unstable one ${\overline{y}}_2=K$.

  Figure 11 shows the evolution of the species population $y_n$ and that of the environment carrying capacity $K_n$ if the population is initially $y_0=250$. The results displayed in Figure 11 illustrates the fact that the equilibrium point ${\overline{y}}_1=0$ is globally asymptotically stable while the equilibrium point ${\overline{y}}_2=500$ is unstable for $c=0$ as Section 5 points out.

(ii) In the second case ${c_n}_0^{\infty }$ is given by (102) with the values for the parameters ${\epsilon }_3=0.75$ and ${\rho }_3=1-{\epsilon }_{3}1-a/4K=6.25x{10}^{-5}$ and the initial condition $c_0=0.01$. In this way, $c={\lim }_{n⟶\infty }{c}_n={\rho }_3/1-{\epsilon }_3=2.5x{10}^{-4}$ and then, the conditions that $a<1$, $K>0$ and $c=1-a/4K$ are fulfilled so that the modified Beverton–Holt equation has two equilibrium points, namely, the stable point ${\overline{y}}_1=0$ (extinction) and the unstable one ${\overline{y}}_2=1-a/2c=1000$. Figure 12 shows the evolution of the species population $y_n$ and that of the environment carrying capacity $K_n$ if the population is initially $y_0=1050$. The results displayed in Figure 12 illustrate the fact that the equilibrium point ${\overline{y}}_1=0$ is globally asymptotically stable while the equilibrium point ${\overline{y}}_2=1000$ is unstable as Theorem 4 (i) establishes.

(iii) In the third case ${c_n}_0^{\infty }$ is given by (102) with the values for the parameters ${\epsilon }_3=0.75$ and ${\rho }_3=0.001>1-{\epsilon }_{3}1-a/4K=6.25x{10}^{-5}$ and the initial condition $c_0=0.01$. In this way, $c={\lim }_{n⟶\infty }{c}_n={\rho }_3/1-{\epsilon }_3=0.004$ and then, the conditions that $a<1$, $K>0$ and $c>1-a/4K=2.5x{10}^{-4}$ are fulfilled so that the modified Beverton–Holt equation has only one equilibrium point, namely, the stable point ${\overline{y}}_1=0$ (extinction). Figure 13 shows the evolution of the species population $y_n$ and that of the environment carrying capacity $K_n$ if the population is initially $y_0=250$. The results displayed in Figure 13 illustrate the fact that the unique equilibrium point ${\overline{y}}_1=0$ is globally asymptotically stable, that is asymptotically stable for any finite initial condition, as Theorem 4 (ii) establishes.

(iv) In the fourth case ${c_n}_0^{\infty }$ is given by (102) with the values for the parameters ${\epsilon }_3=0.75$ and ${\rho }_3=3x{10}^{-5}<1-{\epsilon }_{3}1-a/4K=6.25x{10}^{-5}$ and the initial condition $c_0=0.01$. In this way, $c={\lim }_{n⟶\infty }{c}_n={\rho }_3/1-{\epsilon }_3=1.2x{10}^{-4}$ and then, the conditions that $a<1$, $K>0$ and $0<c<1-a/4K=2.5x{10}^{-4}$ are fulfilled so that the modified Beverton–Holt equation has three equilibrium points, namely, the stable point ${\overline{y}}_1=0$ (extinction) and the unstable ones ${\overline{y}}_2=1-a-\sqrt{{1-a}^2-4ck1-a}/2c\approx 581>K=500$ and ${\overline{y}}_3=1-a-\sqrt{{1-a}^2-4ck1-a}/2c\approx 3586<1-a/c\approx 4167$. Figure 14 shows the evolution of the species population $y_n$ and that of the environment carrying capacity $K_n$ if the population is initially $y_0=700$. The results displayed in Figure 14 illustrate the fact that the equilibrium point ${\overline{y}}_1=0$ is globally asymptotically stable while the other ones are unstable as Theorem 4 (ii) establishes.

[figure(s) omitted; refer to PDF]

7. Conclusions

This paper has discussed a generalized time-varying Beverton–Holt equation which considers the presence of positive or negative harvesting and, eventually, a quadratic-type penalty for the population excess. Such a term takes account for the potential internal competence between the cohort individuals for food, refuge, etc. The harvesting action (describing hunting/fishing actions) is considered jointly with eventually present independent consumption (describing migrations from outside of the habitat to inside or vice-versa). It is seen that the presence of the penalty term can translate into the presence of two other positive equilibrium points. Some particular stability results have been also derived for the stationary equation, which arises when its parameterizing sequences converge, for the case of small levels of population by introducing a term taking account for the Allee effect. The paper has also designed some species evolution control laws by monitoring the harvesting action and has discussed the influence in the stability results of considering a modelling function of Allee effect which makes difficult growing or even can cause extinction for small numbers of reproductive individuals.

The equilibrium points of the stationary solution in the presence and absence of harvesting action have been characterized and their local asymptotic stability properties have been investigated in the case of intrinsic growth rate exceeding unity and eventual execution of harvesting actions and in the case of the intrinsic growth rate being less than unity. Some numerical examples have been also discussed.

Acknowledgments

The authors are grateful to the Basque Government for its support through Grant no. IT1555-22 and to MCIN/AEI 269.10.13039/501100011033 for Grant no. PID2021-1235430B-C21/C22. The authors are also grateful to the referees by their useful comments.

References

[1] F. Zimmermann, K. Enberg, M. Mangel, "Density-independent mortality at early life stages increases the probability of overlooking and underlying stock-recruitment relationship," ICES Journal of Marine Science, vol. 78 no. 6, pp. 2193-2203, 2021.

[2] S. Stevic, "A short proof of the Cushing-Henson conjecture," Discrete Dynamics in Nature and Society, vol. 2006, 2006.

[3] M. Bohner, S. H. Streipert, "The second Cushing-Henson conjecture for the Beverton-Holt q-difference equation," Opuscula Mathematica, vol. 37 no. 6, pp. 795-819, DOI: 10.7494/opmath.2017.37.6.795, 2017.

[4] M. De la Sen, S. Alonso-Quesada, "A Control Theory point of view on Beverton–Holt equation in population dynamics and some of its generalizations," Applied Mathematics and Computation, vol. 199 no. 2, pp. 464-481, DOI: 10.1016/j.amc.2007.10.021, 2008.

[5] M. De la Sen, S. Alonso-Quesada, "Model-matching-based control of the Beverton-Holt equation in ecology," Discrete Dynamics in Nature and Society, vol. 2008,DOI: 10.1155/2008/793512, 2008.

[6] M. De la Sen, "About the properties of a modified generalized Beverton-Holt equation in ecology models," Discrete Dynamics in Nature and Society, vol. 2008,DOI: 10.1155/2008/592950, 2008.

[7] M. De la Sen, "The environment carrying capacity is not independent of the intrinsic growth rate for subcritical spawning stock biomass in the Beverton-Holt equation," Ecological Modelling, vol. 204 no. 1-2, pp. 271-273, DOI: 10.1016/j.ecolmodel.2006.11.029, 2007.

[8] M. De la Sen, A. Ibeas, S. Alonso-Quesada, A. J. Garrido, I. Garrido, "On the properties of a class of impulsive competition Beverton-Holt equations," Applied Sciences, vol. 11 no. 19,DOI: 10.3390/app11199020, 2021.

[9] T. Quinn, "Population Dynamics," Encyclopedia Of Environmetrics, 2012.

[10] Z. Alsharawi, M. B. H. Rhouma, "The Beverton-Holt model with periodic and conditional harvesting," Journal of Biological Dynamics, vol. 3 no. 5, pp. 463-478, DOI: 10.1080/17513750802601041, 2009.

[11] Y. Li, J. Li, "Discrete-time models for releases of sterile mosquitoes with Beverton-Holt type of survivability," Ricerche di Matematica, vol. 67 no. 1, pp. 141-162, DOI: 10.1007/s11587-018-0361-4, 2018.

[12] Y. Li, J. Li, "Stage-structured discrete-time models for interacting wild and sterile mosquitoes with beverton-holt survivability," Mathematical Biosciences and Engineering, vol. 16 no. 2, pp. 572-602, DOI: 10.3934/mbe.2019028, 2019.

[13] S. Al-Nasir, A. H. Lafta, "Analysis of a harvested discrete-time biological models," International Journal of Nonlinear Analysis and Applications, vol. 12 no. 2, pp. 2235-2246, 2021.

[14] G. R. J. Gaut, K. Goldring, F. Grogan, C. Haskell, R. J. Sacker, "Difference Equations with the Allee effect and the periodic sigmoid Beverton-Holt equation revisited," Journal of Biological Dynamics, vol. 6 no. 2, pp. 1019-1033, DOI: 10.1080/17513758.2012.719039, 2012.

[15] A. J. Harry, C. M. Kent, V. L. Kocic, "Global behavior of solutions of a periodically forced Sigmoid Beverton–Holt model," Journal of Biological Dynamics, vol. 6 no. 2, pp. 212-234, DOI: 10.1080/17513758.2011.552738, 2012.

[16] E. J. Bertrand, M. R. S. Kulenovic, "Global dynamics of higher-order transcendental-type generalized Beverton-Holt equations," International Journal of Differential Equations, vol. 13 no. 2, pp. 71-84, 2018.

[17] T. Khyat, M. R. S. Kulenovic, "The invariant curve caused by Neimark-Sacker bifurcation of a perturbed Beverton-Holt difference equation," International Journal of Differential Equations, vol. 12 no. 2, pp. 267-280, 2017.

[18] M. R. S. Kulenovic, S. Moranjkic, Z. Nurkanovic, "Global dynamics and bifurcation of a perturbed Sigmoid Beverton-Holt difference equation: m. R. S. KULENOVIĆ, S. Moranjkić and Z. nurkanovic," Mathematical Methods in the Applied Sciences, vol. 39 no. 10, pp. 2696-2715, DOI: 10.1002/mma.3722, 2016.

[19] C. M. Kent, "Attenuance and resonance in a periodically forced sigmoid Beveron-Holt model," International Journal of Difference Equations, vol. 7 no. 1, pp. 35-50, 2012.

[20] M. Bohner, J. Mesquita, S. Streipert, "The Beverton–Hold model on isolated time scales," Mathematical Biosciences and Engineering, vol. 19 no. 11, pp. 11693-11716, DOI: 10.3934/mbe.2022544, 2022.

[21] V. L. Kocic, "A note on the nonautonomous delay Beverton-Holt model," Journal of Biological Dynamics, vol. 4 no. 2, pp. 1131-1139, DOI: 10.1080/17513750902803588, 2010.

[22] H. Sparholt, "Fish species interaction in the Baltic sea," Dana, vol. 10, 1994.

[23] M. P. Sissenwine, N. Daan, "Multispecies models relevant to management of living resources," Actes ICES mar Science, vol. 193, 1989.

[24] V. Berinde, "Iterative approximation of fixed points," Lecture Notes in Mathematics, vol. 1912, 2006.