The General Analytic Expression of a Harvested Logistic Model with Slowly Varying Coefﬁcients

: The harvested logistic model with a slow variation in coefﬁcients has been considered. Two cases, which depend on the harvest rate, were identiﬁed. The ﬁrst one is when the harvest is subcritical, where the population evolves to an equilibrium. The other is supercritical harvesting, where the population decreases to zero at ﬁnite times. The single analytic approximate expression, which is capable of describing both harvesting cases, is readily and explicitly obtained using the multi-time scaling method together with the perturbation approach. This solution ﬁts for a wide range of coefﬁcient values. In addition, such an expression is validated by utilizing numerical computations, which are obtained by using the fourth-order Runge–Kutta technique. Finally, the comparison shows a very good agreement between the two methods.


Introduction
Utilizing the harvested model to characterize the behaviour of a single species population has involved extensive investigations over the years. Such a model is a simple logistic model containing a harvesting contribution that is proportional to the population [1]. The harvested model has been implemented in numerous biological applications, such as the fishery industry [2], modeling the evolution of the Sandhill crane population [3], and other applications, shown in [3,4].
The exact solution of the harvested model with constant parameters is well-known in the literature-see, for example, [3,5]. In reality, the regular fluctuation of the environment has a significant influence on the evolution of the population. Accordingly, there are many studies that been conducted by considering variable environmental parameters. For example, Rosenblat [6] and Legović et al. [7] investigated population models in view of a periodic fluctuation in the environment. In cases such as those in [6,7] and among others [8,9], the exact solution of non-autonomous systems is often impossible, and numerical approximate methods must be utilized. However, numerical results provide limited information for the population's evolution by responding to particular input data-for example, [10]. Thus, Shepherd et al. [11] considered obtaining an analytic expression for the logistic model with variable carrying capacities, while the growth rate is held fixed. In their calculations, it has been assumed that there is a slow variation in the carrying capacity compared to the total change in the population. Such assumptions allow applying a multi-time scaling method to obtain the analytic expression of the system. Similarly, Grozdanvski et al. [12] modified the technique of those in [11] to deal with the logistic model with two slowly varying coefficients. The analytic approximation is successfully obtained and compared well with the numerical result. In addition, Alharbi [13] investigated a similar model but in view of the asymmetric variations in coefficients. Thus, the application of the multi-time scaling method has successfully been developed to obtain the analytic expression of the system. In this analysis, such an expression has been compared with numerical results and has shown very good agreement. For the logistic model with a harvesting coefficient that is proportional to the population, there are two situations of the population's evolution that may be recognised, which involve subcritical harvesting where the population evolves to a surviving state and supercritical harvesting where the population decreases to zero in finite time. Grozdanovski et al. [14] investigated such models and obtained the analytic expressions of subcritical and supercritical harvesting. However, in their calculations, they obtained the expression of the two situations separately, which leads to obtaining two definitions for the time variation: one for subcritical harvesting and the other for supercritical harvesting. Applying an analogous analysis as those in [14], the single-species logistic model with Holling type II functional responses was investigated by Idlango et al. [15]. The subcritical and supercritical harvesting expansions have been obtained analytically but only by considering a slow variation in the harvesting rate, whereas the other coefficients are held fixed. Furthermore, analytic and theoretical approximations along with the dynamical system analysis of the predator-prey model with Holling type III functional responses have been obtained by Tassaddiq et al. [16]. Their approximations are validated additionally with numerical simulations, and the comparison shows very good agreement.
Therefore, we aim in this paper to investigate the harvested logistic model with variable parameters to display both subcritical and supercritical harvesting behaviours in a single comprehensive expression. The analytic approximation will be obtained by using the multi-time scaling method together with the perturbation approach. Due to the absence of the exact solution of such a system, we consider a numerical computation using a fourth-order Runge-Kutta technique in order to validate the analytic approximation. This comparison shows a favorable agreement between the two methods.

Analytic Analysis
The harvesting model, which arises from the modification on the Logistic model by including density-dependant harvesting, may be written in the following form: whereN(T) represents the population density on T > 0, while P(T), K(T), and H(T) are a positive-valued functions on T > 0 of the growth rate, the carrying capacity, and the harvesting rate, respectively.N 0 is the initial population when T = 0.
In what follows, we first non-dimensionalize the harvesting model (1) by assuming P(T), K(T), and H(T) varying on time scales T P , T K , and T H , respectively. Thus, functions P(T), K(T), and P(T) may be expressed in dimensionless forms as follows: where ρ, κ, and η are the dimensionless growth rate, carrying capacity, and harvesting rate, respectively, while P 0 , K 0 , and H 0 denote characteristic constants.
Here, we introduce the dimensionless population and time byn =N/K 0 and t = P 0 T, respectively. Accordingly, dimensionless harvesting model (1) can be expressed as follows: where positive dimensionless constant σ is a ratio of the harvesting to the growth characteristic parameter, while µ is the dimensionless initial population corresponding toN 0 , and it is defined by the following.
If we assume P, K, and H vary on the same time scales, i.e., T P = T K = T H = T * , then function (3) might be reintroduced in the following form: where = (T * P 0 ) −1 .
As model (4) stands, the exact solution of populationn(t) is impossible; hence, numerical computations must be employed to obtain the numerical approximation. However, when we consider that the time scale for P, K, and H largely varies than that of population N, we aim to have a small . Accordingly, the coefficients ρ( t), κ( t), and η( t) vary slowly relative to populationn(t). Then, two time scales may be observed in (4), which are the intrinsic time scale, t, and slow time scale, t, in view that the analytic approximate expression may be obtained by using the multi-time scaling method together with the perturbation approach based on the small .

Multi-Time Scaling
Following the justification above, we construct a multi-time scaling method by presenting the generalized ordinary time variable, t o , corresponding to t, and the slow time variable, t s , corresponding to t, as follows: where subscripts o and s refer to the intrinsic time scale and slow time scale, respectively. Note that w(t s ) is expected to be a positive, continuous, and differentiable function on With respect to t s = t, differentiating (5) yields the following: and note, from (6), that the relationship between t o and t is governed by function w (t s ). Thus, we assume w (t s ) > 0 to ensure such relationships are monotonic. Considering the justification above, population functionn(t, ) may be expressed as a function, n(t o , t s , ), of these variables. Thus, utilizing the chain rule, system (4) becomes a partial differential equation for the unknown n(t o , t s , ) as follows: where D o and D s represent the partial derivatives of t o and t s , respectively. Therefore, multi-time scaling is successfully applied and such an equation is analogous to that of (4). Note that, since appears explicitly in (7) , this allows us to construct the perturbation approach as → 0, which leads to explicitly obtaining the analytic approximate expression for the population of (4) that is valid for all t > 0.

Perturbation Analysis
Here, applying the perturbation approach based on a small to (7) yields the following: expressing n(t o , t s , ) in the power of .
Substituting (8) into (7) provides the following upon equating the similar power of : and this determines the leading and higher order terms that are independent of . Solving (9) for n 0 (t o , t s ) obtains the following: where A(t s ) is an arbitrary function of t s that is determined later, while θ(t s ) is given by the following.
In a similar manner, substituting (11) in (10), we obtain the particular solution for n 1 (t o , t s ) upon solving the following: where the primes in A (t s ) and κ (t s ) represent the derivative with respect to t s , while ν(t o , t s ) is given by the following.
Observe that by fixing t s while t o → ∞, the rate of convergence of (11) is exponential, while it is not in (13) due to the presence of the quadratic function of t o . Consequently, we consider aligning the rate of convergence of both n 0 (t o , t s ) and n 1 (t o , t s ) to reach their limits at the same time. In order to perform this, an analogous line of reasoning of those in [13,14,17] may be applied to eliminate the quadratic function of t o in (13) by separately setting the coefficients of t o and t 2 o to zero. This yields the following.
In order to satisfy both equations of (15), we choose A (t s ) = 0 for the first and, hence, A(t s ) = c, while we select θ (t s ) = 0 in the second; hence, θ(t s ) is constant, which will be chosen later.
Rearranging the expression of θ(t s ), given in (12), yields the following: which shows that w (t s ) is defined as a difference between the growth and harvesting rate function. However, it is important to remember that w (t s ) must always be positive as mentioned earlier. To maintain this condition, we carefully have to choose the value of θ(t s ).
Observe that σ is the positive characteristic parameter and define σ = H 0 P 0 as a ratio between the harvesting and growth rate parameters. Accordingly, when ρ(t s ) > η(t s ) on interval 0 ≤ t s < ∞, σ < 1 and the population survives. In contrast, for η(t s ) > ρ(t s ) on interval 0 ≤ t s < ∞, σ > 1 and the population dies out. Consequently, choosing θ(t s ) = 1 − σ always maintains condition w (t s ) > 0 no matter what the situation of the population is. From this observation, (16) is rewritten as follows.
By considering A(t s ) = c and θ(t s ) = 1 − σ, (11) and (13) become the following: where the prime represents the derivative with respect to t s . Applying (17) to (5) yields the expression of time scale t o as follows: defining the variable t o completely in terms of growth and harvesting rate functions.
To summarize our analysis, applying (18) and (19) into the expansion (8) yields the analytic approximate expression as follows: where the prime denotes to the derivative with respect to t s , and γ(t s ) and Ψ(t o , ts) are given by the following: while t o is as given in (20), and c is a constant. As expansion (8) stands, it is constructed with O(1) and O( ). Accordingly, we assume that c in (21) takes a similar form, i.e., c = c 0 + c 1 . Thus, applying c expansion to (21) and expanding the powers of produce the general analytic approximate expression of the harvesting model, which is capable of describing both subcritical and supercritical harvesting as follows: where γ is as given in (22) and ξ(t o , t s ) is given by the following: while c 0 and c 1 are obtained by applying initial condition (4) into (24) and then equating the coefficients in powers of . This produces the following when solved. ,

Results and Discussion
In this paper, the multi-time scaling method was successfully applied to the harvesting model, and expansion (24) represents the result. Such an expansion (24) is capable of describing the population behaviour in the two situations, which are subcritical and supercritical harvesting. Note that the perturbation expansion (8) has been limited to the first two terms due to their remarkable accuracy without a swell of expressions. In addition, when → 0 is considered with constant coefficients, expression (24) reduces to the exact solution.
For σ < 1, the population survives to the equilibrium state, which can be obtained from (24) as t → ∞. This is given as follows: varying slowly as t s → ∞. Note from (26) that for a large class of time, the variation in the carrying capacity, the growth rate, and the harvesting rate becomes constant; hence, their slopes will be zero, i.e., κ (t s ) = ρ (t s ) = η (t s ) = 0. Accordingly, population n(t o , t s , ) eventually tends to κ(t s ) − σ η(t s ) ρ(t s ) , varying close to the carrying capacity. In contrast, when σ > 1, ρ(t s ) < η(t s ); hence, the population is extinct at a finite time.
Fixing the growth and the harvesting rate, say ρ(t s ) = η(t s ) ≡ 1, the time-variable relationship (20) yields t o ≡ 1 t s and t s = t, which is realistic for describing the behaviour of the population by considering a slow variation in a biotic environment.
In order to validate our analytic approximation, we consider solving our problem numerically by using the fourth-order Runge-Kutta technique with appropriate choices of carrying capacity, growth rate, and harvesting rate functions. Consequently, Figures 1a and 2a, were generated using expansion (24) together with the numerical result, which shows very good agreement between the analytic (solid curve) and the numerical (dotted curve) solutions, and this agreement only falls off at higher values of . Graphically, Figure 3a confirms that the accuracy between the two methods fails at = 0.5 where the absolute difference between the numerical and analytic data reaches 2.5% according to the curve peak value in Figure 3b. Figure 1a displays the population evolution from initial value µ = 0.1 to the limiting state, given in (26), for slow and periodic functions of κ( t), ρ( t) and η( t), where = 0.2. Figure 1b was generated by computing the absolute difference between the analytic and numerical solutions as t → ∞, and it displays favorable accuracies by showing a maximum discrepancy that is less than 6 × 10 −3 . Figure 2a shows the extinct population where the curve decreases from initial value µ = 0.1 to zero at a finite time for the periodic functions of κ( t), ρ( t), and η( t), where = 0.2. Figure 2b shows that the agreement between the analytic and numerical results as good t → ∞, and the absolute difference between the two methods reaches 1.2 × 10 −6 .

Conclusions
In this study, we investigated the harvested logistic model with a slow variation in the coefficients. Two cases, which depend on the harvest rate, have been determined. The first one is when the harvest is subcritical where the population evolves to the equilibrium. The other is supercritical harvesting where the population decreases to zero in a finite time period. The single analytic approximate expression, which is capable of describing both harvesting cases, is readily and explicitly obtained using the multi-time scaling method together with the perturbation approach. This solution fits for a wide range of coefficients values. In addition, such an expression is validated by utilizing numerical computations, which are obtained by using the fourth-order Runge-Kutta technique. Thus, the comparison between the two methods shows very good accuracy, as → 0.