Global Dynamics of an Exploited Prey-Predator Model with Constant Prey Refuge

Abstract

This paper describes a prey-predator model with Holling type II functional response incorporating constant prey refuge and harvesting to both prey and predator species. We have analyzed the boundedness of the system and existence of all possible feasible equilibria and discussed local as well as global stabilities at interior equilibrium of the system. The occurrence of Hopf bifurcation of the system is examined, and it was observed that the bifurcation is either supercritical or subcritical. Influences of prey refuge and harvesting efforts are also discussed. Some numerical simulations are carried out for the validity of theoretical results.

1. Introduction

Prey-predator models are of great interest to researchers in mathematics and ecology because they deal with environmental problems such as community’s morbidity and how to control it and optimal harvest policy to sustain a community. In the physical sciences, generic models can be constructed to explain a variety of phenomena. However, in the life sciences, a model only describes a particular situation. So a variety of models are needed due to the complexity of the ecosystem. Theoretical and numerical studies of these models are able to give us an understanding of the interactions that are taking place.

While investigating biological phenomena, there are many factors which affect dynamical properties of biological and mathematical models. One of the familiar nonlinear factors is the functional response. In population dynamics, a functional response of the predator to the prey density refers to the change in the density of prey attached per unit time per predator as the prey density changes [1]. Holling [2] suggested three different kinds of functional response for different kinds of species to model the phenomena of predation, which made the classical Lotka-Volterra system more realistic. The other factor which affects dynamical properties is harvesting. The subject of harvesting in predator-prey systems has been of interest to economists, ecologists, and natural resource managers for some time now. There are basically three types of harvesting reported in the literature: (i) constant harvesting where a constant number of individuals are harvested per unit time, (ii) proportional harvesting where the number of individuals harvested per unit time is proportional to the current population, and (iii) nonlinear harvesting. Most research has focused attention on optimal exploitation, guided entirely by profits from harvesting. First of all, depending on the nature of the applied harvesting strategy, the long-run stationary density of population may be significantly smaller than the long-run stationary density of a population in the absence of harvesting. Therefore, while a population can, in the absence of harvesting, be free of extinction risk, harvesting can lead to the incorporation of a positive extinction probability and therefore, to potential extinction in finite time. Secondly, if a population is subjected to a positive extinction rate, then harvesting can drive the population density to a dangerously low level at which extinction becomes a sure thing no matter how the harvester affects the population afterwards.

Effects of harvesting on various types of prey-predator models have been considered by many researchers. Mathematical modeling with harvesting renewable resources started with the studies of Clark [3, 4]. Xiao et al. [5] have investigated the dynamics of a system with constant harvesting on the predator, whereas Leard et al. [6] have studied the dynamics of ratio-dependent models that include nonconstant harvesting on the prey. Xia et al. [7] studied the combined effects of harvesting and time delay. Recently, Pahari and Kar [8] analyzed a prey-predator fishery model with harvesting where they applied a tax to regulate the fishery.

The existence of prey refuges can clearly have important effects on the coexistence of predators and their prey. According to Taylor [9], the different kinds of refuges can be arranged into three types: (a) those which provide permanent spatial protection for a small subset of the prey population, (b) those which provide temporary spatial protection, and (c) those which provide a temporal refuge in numbers, that is, decrease the risk of predation by increasing the abundance of vulnerable prey. It is also observed that refuge has a stabilizing effect on the equilibrium for a simple Lotka-Volterra model. The problem of predator-prey interactions under a prey refuge has been studied by some authors. McNair [10] showed that several kinds of refuges could exert a locally destabilizing effort and create stable, large-amplitude oscillations which would damp out if no refuge was present. McNair [11] obtained that a prey refuge with legitimate entry existing dynamics was quite capable of amplifying rather than damping predator-prey oscillation. Kar [12] proposed a predator-prey model incorporating a prey refuge and independent harvesting on either species. He showed that, using the harvesting efforts as control, it is possible to break the cyclic behavior of the system. Wang et al. [13] considered a prey-predator system where individuals from a prey fish population could hide in holes where predators were unable to enter. Ji and Wu [14] considered a predator-prey model incorporating a constant prey refuge and a constant-rate prey harvesting and showed that the system is controlled by using constant harvesting or prey refuge. Prepredator interactions with prey refuge also may be found in the work of Huang et al. [15], Wang and Pan [16], Kar et al. [17], and so forth.

The objective of this paper is to study the combined effects of harvesting and constant prey refuge on the dynamics of predator-prey model. We also examined the existence of a Hopf bifurcation and determined the conditions for which the bifurcation is either supercritical or subcritical.

González-Olivares and Ramos-Jiliberto [18] studied the following prey-predator model with constant number of prey using refuges as shown in the Figure 1: 𝑑𝑥𝑥𝑑𝑡=𝛼𝑥1−𝑘−𝛽(𝑥−𝑚)𝑦(),1+𝑎𝑥−𝑚𝑑𝑦()𝑦𝑑𝑡=−𝑑𝑦+𝑐𝛽𝑥−𝑚(),1+𝑎𝑥−𝑚 (1) with 𝑥(0)>𝑚, 𝑦(0)>0, where 𝑥 and 𝑦 denote the prey and predator populations, respectively, at any time 𝑡 and 𝛼 and 𝑘, respectively, represent the intrinsic growth rate and carrying capacity of the prey. 𝑚 is a constant, which represents number of prey which seeks refuge from predation. 𝑑 is the death rate of the predator and 𝑐 is the conversion factor. The term 𝛽𝑥∣(1+𝑎𝑥) denotes the functional response of the predator, which is termed Holling type II response function (see [2]). They examined the local stability of equilibria and existence of limit cycle. Chen et al. [19] considered the same model and examined also global stability and uniqueness of limit cycle. Also Ji and Wu [14] considered the prey-predator model with constant prey refuge and a constant-rate prey harvesting ℎ>0, 𝑑𝑥𝑥𝑑𝑡=𝛼𝑥1−𝑘−𝛽(𝑥−𝑚)𝑦1+𝑎(𝑥−𝑚)−ℎ,𝑑𝑦()𝑦𝑑𝑡=−𝑑𝑦+𝑐𝛽𝑥−𝑚(),1+𝑎𝑥−𝑚 (2) and studied instability and global stability of the equilibria and uniqueness of limit cycle and showed also the influence of constant prey refuge and constant-rate prey harvesting.

Figure 1 

Conceptual diagram of the Rosenzweig-MacArthur model with prey refuge. Small quadrilateral and oval shapes represent the prey and predator, respectively. The vertical dashed line represents the boundary of a refuge, where prey species are invulnerable to predators. The refuge could protect either a constant number or a constant fraction of the total prey population.

Motivated by the paper of González-Olivares and Ramos-Jiliberto [18], we considered the following prey-predator system: 𝑑𝑥𝑥𝑑𝑡=𝛼𝑥1−𝑘−𝛽(𝑥−𝑚)𝑦1+𝑎(𝑥−𝑚)−𝑞₁𝐸₁𝑥,𝑑𝑦()𝑦𝑑𝑡=−𝑑𝑦+𝑐𝛽𝑥−𝑚1+𝑎(𝑥−𝑚)−𝑞₂𝐸₂𝑦, (3) where 𝐸₁≥0 and 𝐸₂≥0 denote the harvesting efforts for the prey and predator, respectively. 𝑞₁𝐸₁𝑥 and 𝑞₂𝐸₂𝑦 represent the catch of the prey and predator population, where 𝑞₁ and 𝑞₂ represent the catchability coefficients, respectively. 𝛼, 𝛽, 𝑘, 𝑑, 𝑎, 𝑐, 𝑚, 𝑞₁, and 𝑞₂ are positive constants.

Our paper is organized in the following way. In Section 2, we have discussed the boundedness of the solutions of system (3). All possible equilibria of the system (8) and the stability criterion at those equilibria and uniqueness of limit cycles at interior equilibrium are discussed in Section 3. In Section 4, we have discussed Hopf bifurcation at interior equilibrium. Influence of refuge parameter 𝑚 and harvesting efforts 𝐸₁ and 𝐸₂ are discussed in Section 5. Numerical simulations are given in Section 6. A brief concluding remark is given in Section 7.

2. Boundedness of the System

Boundedness of a model guarantees its validity. The following theorem establishes the uniform boundedness of the model system (3).

Theorem 1. All solutions of the system (3) which start in 𝑅²₊ are uniformly bounded.

Proof. Let (𝑥(𝑡),𝑦(𝑡)) be any solution of the system with positive initial conditions 𝑥(0)>𝑚 and 𝑦(0)>0.
Now we define the function 𝑊=𝑥+𝑦/𝑐.
Therefore, time derivative gives 𝑑𝑊=𝑑𝑡𝑑𝑥+1𝑑𝑡𝑐𝑑𝑦𝑑𝑡=𝛼𝑥−𝛼𝑥²𝑘−𝑞₁𝐸₁𝑥−𝑑𝑦𝑐−𝑞₂𝐸₂𝑦𝑐. (4) Now for each 𝑣>0, we have 𝑑𝑊𝑑𝑡+𝑣𝑊=𝛼𝑥−𝛼𝑥²𝑘−𝑞₁𝐸₁𝑥+𝑣𝑥+𝑣𝑦𝑐−𝑑𝑦𝑐−𝑞₂𝐸₂𝑦𝑐. (5) We have 𝑑𝑊/𝑑𝑡+𝑣𝑊≤(𝑘/4𝛼)(𝛼+𝑣−𝑞₁𝐸₁)²−(1/𝑐)(𝑑+𝑞₂𝐸₂−𝑣).
Let us choose 𝑣>𝑑+𝑞₂𝐸₂; then the right-hand side is positive. As we assume that both 𝐸₁ and 𝐸₂ are bounded, the right-hand side is bounded for all (𝑥,𝑦)∈𝑅²₊.
Thus we choose a 𝑢>0 such that 𝑑𝑊/𝑑𝑡+𝑣𝑊<𝑢.
Applying the theory of differential inequality [20], we obtain <𝑢0<𝑊𝑥,𝑦𝑣1−𝑒^−𝑣𝑡+𝑊𝑥(0),𝑦(0)𝑒^𝑣𝑡,𝑢for𝑡⟶∞,0<𝑊<𝑣. (6) Thus, all solutions of the system (3) that start in 𝑅²₊ are confined to the region 𝐵, where 𝐵=𝑥,𝑦∈𝑅²₊𝑢∶𝑊=𝑣+𝜀,forany𝜀>0. (7) This completes the theorem.

Now for simplicity, let us introduce 𝑋=𝑥−𝑚; then the system (3) of equations changes to (still denote 𝑋=𝑥) 𝑑𝑥()𝑑𝑡=𝛼𝑥+𝑚1−𝑥+𝑚𝑘−𝛽𝑥𝑦1+𝑎𝑥−𝑞₁𝐸₁(𝑥+𝑚),𝑑𝑦𝑑𝑡=−𝑑𝑦+𝑐𝛽𝑥𝑦1+𝑎𝑥−𝑞₂𝐸₂𝑦. (8)

3. The Steady States and Their Stability

We now study the existence and nature of the steady states. Particularly we are interested in the interior equilibrium of the system. To begin with, we list all possible steady states of the system (8) as follows:(i)equilibrium in the absence of predator 𝑃₁(𝑥₁,0), where 𝑥₁=(𝑘𝛼−𝑚𝛼−𝑞₁𝐸₁𝑘)/𝛼 exists provided 𝑘𝛼>𝑚𝛼+𝑞₁𝐸₁𝑘,(ii)the interior equilibrium 𝑃^∗(𝑥^∗,𝑦^∗), where 𝑥^∗=𝑑+𝑞₂𝐸₂𝛽𝑐−𝑎𝑑−𝑎𝑞₂𝐸₂,𝑦^∗=1+𝑎𝑥^∗𝛽𝑘𝑥^∗𝛼𝑥^∗+𝑚𝑘−𝑚−𝑥^∗−𝑞₁𝐸₁𝑘𝑥^∗+𝑚. (9)If 𝛽𝑐>𝑎𝑑+𝑎𝑞₂𝐸₂ holds and any of the following conditions hold:(H1)𝐴>0, 𝐵>0, 𝐴²≥4𝐵, and (𝐴−𝐴²−4𝐵)/2<𝑚<(𝐴+𝐴²−4𝐵)/2,(H2)𝐴>0, 𝐵<0, and 0<𝑚<(𝐴+𝐴²−4𝐵)/2,(H3)𝐴<0, 𝐵<0, and 0<𝑚<(𝐴+𝐴²−4𝐵)/2,where 𝐴=𝑘−2𝑥^∗−𝑞₁𝐸₁𝑘/𝛼 and 𝐵=(𝑥^∗)²−𝑘𝑥^∗+𝑞₁𝐸₁𝑘𝑥^∗/𝛼, then the interior equilibrium 𝑃^∗(𝑥^∗,𝑦^∗) exists.

For the stability analysis of the equilibria, let the variational matrix of the system (8) at an arbitrary point (𝑥,𝑦) be 𝑀(𝑥,𝑦).

At 𝑃₁(𝑥₁,0), the eigenvalues of the variational matrix 𝑀(𝑥₁,0) are 𝐸₁𝑞₁−𝛼 and 𝑐𝛽𝑥₁/(1+𝑎𝑥₁)−𝑑−𝑞₂𝐸₂.

Now if (𝛼/𝑞₁){(𝑘−𝑚)/𝑘−(𝑑+𝑞₂𝐸₂)/𝑘(𝑐𝛽−𝑎𝑑−𝑎𝑞₂𝐸₂)}<𝐸₁<𝛼/𝑞₁ and 𝛽𝑐>𝑎𝑑+𝑎𝑞₂𝐸₂, then both eigenvalues are negative and hence 𝑃₁(𝑥₁,0) is locally asymptotically stable.

To determine the local stability character of the interior equilibrium 𝑃^∗(𝑥^∗,𝑦^∗), we compute the variational matrix 𝑀(𝑥^∗,𝑦^∗) at (𝑥^∗,𝑦^∗).

The characteristic equation of the variational matrix 𝑀(𝑥^∗,𝑦^∗) is given by 𝛾²+𝛾2𝑚𝛼𝑘−𝑀₁+𝑀₂=0, (10) where 𝑀₁=𝛼−𝑞₁𝐸₁−𝛽𝑦^∗(1+𝑎𝑥)²−2𝑥^∗𝛼𝑘,𝑀₂=𝑐𝛽²𝑥^∗𝑦^∗(1+𝑎𝑥^∗)³>0. (11) Routh-Hurwitz criterion states that all roots of the characteristic equation (10) have negative real parts if (2𝑚𝛼/𝑘−𝑀₁)>0.

Therefore, interior equilibrium 𝑃^∗(𝑥^∗,𝑦^∗) will be asymptotically stable if 𝑚>𝑀₁𝑘/2𝛼, and unstable if 𝑚<𝑀₁𝑘/2𝛼.

Theorem 2. If 𝑃^∗(𝑥^∗,𝑦^∗) exists with 𝑚>𝑀₁𝑘/2𝛼 then 𝑃^∗(𝑥^∗,𝑦^∗) is locally asymptotically stable.

3.1. Global Stability for the Interior Equilibrium

Theorem 3. If 𝑅(𝑥)>0, then the system (8) will be globally asymptotically stable around the interior equilibrium 𝑃^∗(𝑥^∗,𝑦^∗).

Proof. To show the global stability of system (8), we define a Lyapunov function as follows: 𝑉𝑥,𝑦=𝐴₁𝑥−𝑥^∗−𝑥^∗𝑥ln𝑥^∗+𝐵₁𝑦−𝑦^∗−𝑦^∗𝑦ln𝑦^∗, (12) where 𝐴₁ and 𝐵₁ are positive constants to be determined in the subsequent steps.
The time derivative along the trajectories of (8) is 𝑑𝑉𝑑𝑡=𝐴₁𝑥−𝑥^∗𝑥𝑑𝑥𝑑𝑡+𝐵₁𝑦−𝑦^∗𝑦𝑑𝑦𝑑𝑡=𝐴₁𝑥−𝑥^∗𝛼𝑚1+𝑥1−𝑥+𝑚𝑘−𝛽𝑦1+𝑎𝑥−𝑞₁𝐸₁−𝑞₁𝐸₁𝑚𝑥+𝐵₁𝑦−𝑦^∗−𝑑−𝑞₂𝐸₂+𝑐𝛽𝑥1+𝑎𝑥=−𝐴₁𝑥−𝑥^∗2𝑥^∗𝛼𝑥^∗𝑘+𝛼𝑚𝑥𝑚1−𝑘−𝑞₁𝐸₁𝑚𝑥−𝑎𝛽𝑥^∗𝑦^∗(1+𝑎𝑥)(1+𝑎𝑥^∗)+𝛽𝑐𝐵₁−𝐴₁1+𝑎𝑥^∗𝑥−𝑥^∗𝑦−𝑦^∗(1+𝑎𝑥)(1+𝑎𝑥^∗)=−𝐴₁𝑥−𝑥^∗2𝑥^∗𝑅(𝑥)+𝛽𝑐𝐵₁−𝐴₁1+𝑎𝑥^∗𝑥−𝑥^∗𝑦−𝑦^∗(1+𝑎𝑥)(1+𝑎𝑥^∗), (13) where 𝑅(𝑥)=[𝛼𝑥^∗/𝑘+(𝛼𝑚/𝑥)(1−𝑚/𝑘)−𝑞₁𝐸₁𝑚/𝑥−𝑎𝛽𝑥^∗𝑦^∗/(1+𝑎𝑥)(1+𝑎𝑥^∗)].
Choosing 𝐴₁=1 and 𝐵₁=(1+𝑎𝑥^∗)/𝑐, we have 𝑑𝑉𝑑𝑡=−𝑥−𝑥^∗2𝑥^∗𝑅(𝑥). (14) Thus if 𝑅(𝑥)>0, then 𝑑𝑉/𝑑𝑡<0. This completes the proof.

In Figure 2, we show that 𝑅(𝑥) may be positive for some positive value of 𝑥 by using the numerical value of the parameters 𝛼=8, 𝑎=0.5, 𝛽=4, 𝑐=0.8, 𝑑=2, 𝑘=50, 𝑞₁=0.6, 𝑞₂=0.2, 𝐸₁=0.4, 𝐸₂=0.2, and 𝑚=0.75.

Figure 2 

Graph of 𝑅(𝑥).

3.2. Uniqueness of Limit Cycles

It is known that for prey-predator systems, existence and stability of a limit cycle are related to the existence and stability of a positive equilibrium. If the limit cycles do not exist, in this case, the equilibrium is globally asymptotically stable. On the other hand, if the positive equilibrium exists and is unstable, there must occur at least one limit cycle.

Let us consider system (8) in the form 𝑑𝑥(𝑥)(𝑥)(0)𝑑𝑡=𝑥𝑔−𝑦ℎ,𝑥>0,𝑑𝑦𝑑𝑡=𝑦−𝑑−𝑞₂𝐸₂+𝑞(𝑥),𝑦(0)>0, (15) where 𝑔(𝑥)=𝛼(1−𝑥/𝑘)−2𝛼𝑚/𝑘−𝑞₁𝐸₁+(1/𝑥)(𝛼𝑚−𝛼𝑚²/𝑘−𝑞₁𝐸₁𝑚),ℎ(𝑥)=𝛽𝑥/(1+𝑎𝑥), and 𝑞(𝑥)=𝑐𝛽𝑥/(1+𝑎𝑥).

Now we consider the following theorem (see [21]) regarding uniqueness of limit cycles of the previous system.

Theorem 4. Suppose for system (15) that 𝑑𝑑𝑥𝑥𝑔^′(𝑥)+𝑔(𝑥)−𝑥𝑔(𝑥)ℎ^′(𝑥)/ℎ(𝑥)−𝑑−𝑞₂𝐸₂+𝑞(𝑥)≤0 (16) in 0≤𝑥<𝑥^∗ and 𝑥^∗<𝑥≤𝑘. Then the previous system has exactly one limit cycle which is globally asymptotically stable with respect to the set 𝑥,𝑦∶𝑥>0,𝑦>0⧵𝑃^∗𝑥^∗,𝑦^∗. (17)

Following Theorem 4, we may state that when 𝑚<𝑀₁𝑘/2𝛼, system (8) has unique globally stable limit cycle. Thus we see that when the system is unstable, there exists a unique globally stable limit cycle.

4. Hopf Bifurcation

To discuss Hopf bifurcation of the system (8) we take the help of the paper [22].

Let us now consider system (8) in the form 𝑑𝑥(𝑥)𝑓(𝑥),𝑑𝑡=ℎ−𝑦𝑑𝑦ℎ(𝑥)𝑑𝑡=𝑐𝑦−𝑑^′, (18) where 𝑓(𝑥)=(𝛼/𝑘𝛽)(𝑘−𝑥)(1+𝑎𝑥)−(2𝛼𝑚/𝑘𝛽)(1+𝑎𝑥)−(𝑞₁𝐸₁/𝛽)(1+𝑎𝑥)+(1/𝑥𝛽)(𝛼𝑚−𝛼𝑚²/𝑘−𝑞₁𝐸₁𝑚)(1+𝑎𝑥),ℎ(𝑥)=𝛽𝑥/(1+𝑎𝑥), and 𝑑^′=(𝑑+𝑞₂𝐸₂)/𝑐.

Let (𝛿,𝑦_𝛿) be the interior equilibrium of the system (18), where 𝛿=𝑑+𝑞₂𝐸₂𝛽𝑐−𝑎𝑑−𝑎𝑞₂𝐸₂,𝑦_𝛿=1+𝑎𝛿𝛼()𝛽𝑘𝛿𝛿+𝑚)(𝑘−𝑚−𝛿−𝑞₁𝐸₁𝑘(𝛿+𝑚). (19) In this section, we examined the Hopf bifurcation occurring at (𝛿,𝑦_𝛿) taking 𝛿 as a bifurcation parameter.

Variational matrix of the system (18) at (𝛿,𝑦_𝛿) is 𝑊(𝛿)𝑌(𝛿)𝑍(𝛿)0𝐽=, (20) where 𝑊(𝛿)=ℎ(𝛿)𝑓^′(𝛿), 𝑌(𝛿)=−ℎ(𝛿), and 𝑍(𝛿)=𝑐𝑓(𝛿)ℎ^′(𝛿).

The characteristic equation of the variational matrix 𝐽 is given by 𝜇²−𝜇𝑇+𝐷=0, where 𝑇=tr𝐽=𝑊(𝛿) and 𝐷=det𝐽=−𝑌(𝛿)𝑍(𝛿).

Now if 𝑇<0 and 𝐷>0, then (𝛿,𝑦_𝛿) is locally asymptotically stable. On the other hand, if 𝑇>0 and 𝐷>0 then (𝛿,𝑦_𝛿), is unstable. Also 𝑇(𝛿)=0 and the Jacobian matrix 𝐽(𝛿) has a pair of imaginary eigenvalues 𝜇=±𝑖−𝑌(𝛿)𝑍(𝛿). Let 𝜇=𝛽(𝛿)±𝑖𝜔(𝛿) be the roots of 𝜇²−𝜇𝑇+𝐷=0 when 𝛿 is near 𝛿; then 𝛽(𝛿)=𝐴(𝛿)2,𝜔(𝛿)=−4𝑌(𝛿)𝑍(𝛿)−{𝑊(𝛿)}²2,𝛽^′(𝛿)=ℎ(𝛿)𝑓^′(𝛿)^′2=ℎ^′(𝛿)𝑓^′(𝛿)+ℎ(𝛿)𝑓^′′(𝛿)2. (21)𝛿 is a root of 𝑓^′(𝛿)=0  (as ℎ(𝛿)=𝑑^′≠0); that is 𝛿=𝛿, satisfies the equation 2𝛼𝑎𝛿³+2𝛼𝑎𝑚+𝑞₁𝐸₁𝑎𝑘+𝛼−𝑎𝛼𝑘𝛿²+𝛼𝑚𝑘−𝛼𝑚²−𝑞₁𝐸₁𝑘𝑚=0, (22) and also 𝑓^′′(𝛿)≠0.

Therefore, 𝛽^′(𝛿)|_𝛿=𝛿=ℎ(𝛿)𝑓^′′(𝛿)/2≠0 as ℎ(𝛿)≠0.

The system (18) undergoes a Hopf bifurcation at (𝛿,𝑦_𝛿) when 𝛿=𝛿. By further analysis, we determined that the bifurcation is either supercritical or subcritical by the first Lyapunov coefficient (see [23–25]) 𝑎𝛿=116ℎ^′𝛿𝑓^′′′𝛿ℎ𝛿ℎ^′𝛿+2ℎ^′𝛿²𝑓^′′𝛿−𝑓^′′𝛿ℎ^′′𝛿ℎ𝛿. (23) The computation of 𝑎(𝛿) is technical, and the detailed calculations are given in the appendix and only the results are stated in the following theorem.

Theorem 5. The system (18) undergoes a Hopf bifurcation at (𝛿,𝑦_𝛿); the Hopf bifurcation is supercritical and backward (subcritical and forward, resp.) if 𝑎(𝛿)<0(𝑎(𝛿)>0), where 𝑎(𝛿) is defined in (23).

5. Influence of the Parameter 𝑚 and Harvesting Efforts 𝐸₁ and 𝐸₂

Interior equilibrium of the system (8) is 𝑃^∗(𝑥^∗,𝑦^∗), where 𝑥^∗=𝑑+𝑞₂𝐸₂𝛽𝑐−𝑎𝑑−𝑎𝑞₂𝐸₂,𝑦^∗=1+𝑎𝑥^∗𝛽𝑘𝑥^∗𝛼𝑥^∗+𝑚𝑘−𝑚−𝑥^∗−𝑞₁𝐸₁𝑘𝑥^∗+𝑚. (24) Notice that the scaling from system (3) to system (8) is 𝑋=𝑥−𝑚. Let the interior equilibrium of the system (3) be 𝑃^′′(𝑥^′′,𝑦^′′), where 𝑥^′′=𝑑+𝑞₂𝐸₂𝛽𝑐−𝑎𝑑−𝑎𝑞₂𝐸₂𝑦+𝑚,_′′=1+𝑎𝑥^′′−𝑚𝛽𝑘𝑥^′′−𝑚𝛼𝑥^′′𝑘−𝑥^′′−𝑞₁𝐸₁𝑘𝑥^′′. (25) The model (3) is a one-prey-one-predator model with constant prey refuge. Now we discuss the influence of the prey refuge on the model dynamics.

Now 𝑑𝑥^′′/𝑑𝑚=1>0.

Therefore, 𝑥^′′ is a strictly increasing function of 𝑚.

That is, increasing the amount of prey refuge can increase prey population.

To check the effect of 𝑚 on predator, differentiating 𝑦^′′ with respect to 𝑚, we get 𝑑𝑦^′′/𝑑𝑚=−(2𝛼(1+𝑎𝑥^∗)/𝛽𝑥^∗𝑘)(𝑚+𝑥^∗−𝑘/2+𝑞₁𝐸₁𝑘/2𝛼), where 𝑥^∗=(𝑑+𝑞₂𝐸₂)/(𝛽𝑐−𝑎𝑑−𝑎𝑞₂𝐸₂).

Case 1. If 𝑚>𝑘/2−𝑥^∗−𝑞₁𝐸₁𝑘/2𝛼, then 𝑦^′′ is a monotonically decreasing function of 𝑚.
That is, if the value of prey refuge 𝑚 gradually increases above threshold value 𝑘/2−𝑥^∗−𝑞₁𝐸₁𝑘/2𝛼, then the predator population gradually decreases.

Case 2. If 𝑚<𝑘/2−𝑥^∗−𝑞₁𝐸₁𝑘/2𝛼, then 𝑦^′′ is a monotonically increasing function of 𝑚.
That is, if the value of prey refuge 𝑚 gradually decreases below threshold value 𝑘/2−𝑥^∗−𝑞₁𝐸₁𝑘/2𝛼, then the predator population gradually increases.

Case 3. If =𝑘/2−𝑥^∗−𝑞₁𝐸₁𝑘/2𝛼, then𝑦^′′ has a maximum value.

To construct Table 1 and Figure 3, we take 𝛼=8, 𝑎=0.5, 𝛽=4, 𝑐=0.8, 𝑑=2, 𝑘=50, 𝑞₁=0.6, 𝑞₂=0.2, 𝐸₁=0.4, and 𝐸₂=0.2 in appropriate units. From the analysis shown we see that increasing the amount of prey refuge can increase prey population and that increasing the amount of prey refuge can increase the density of predator species, and this happened due to predator species still having enough food for predation with 𝑚 being small, but if the prey refuge is larger than a threshold, that is, as the prey refuge becomes large enough, then the increasing amount of prey refuge can decrease predator species, and this happened due to the loss of food for predator species.

Table 1 

Figure 3 

Change of 𝑥′′ and 𝑦′′ with the change of the refuge parameter 𝑚.

Influence of the Harvesting Efforts E₁ and E₂. Interior equilibrium of the system (8) is 𝑃^∗(𝑥^∗,𝑦^∗), where 𝑥^∗=𝑑+𝑞₂𝐸₂𝛽𝑐−𝑎𝑑−𝑎𝑞₂𝐸₂,𝑦^∗=1+𝑎𝑥^∗𝛽𝑘𝑥^∗𝛼𝑥^∗+𝑚𝑘−𝑚−𝑥^∗−𝑞₁𝐸₁𝑘𝑥^∗+𝑚. (26) At first let 𝑚 and 𝐸₂ be fixed.

It is observed that 𝑥^∗ is independent of the parameter 𝐸₁, whereas the value of 𝑦^∗ depends on 𝐸₁. Therefore 𝐸₁ has no effect on the interior equilibrium level on 𝑥.

To check the effect of 𝐸₁ on predator, differentiating 𝑦^∗ with respect to 𝐸₁, we get 𝑑𝑦^∗𝑑𝐸₁𝑞=−¹1+𝑎𝑥^∗𝑚+𝑥^∗𝛽𝑥^∗<0. (27)

This shows that 𝑦^∗ is a strictly decreasing function of 𝐸₁.

To construct Table 2 and Figure 4, we take 𝛼=8, 𝑎=0.5, 𝛽=4, 𝑐=0.8, 𝑑=2, 𝑘=50, 𝑞₁=0.6, 𝑞₂=0.2, 𝑚=0.35, and 𝐸₂=0.2 in appropriate units. From the analysis shown we see that as the harvesting effort 𝐸₁ increases, the prey species remain unchanged but predator species decrease and this happens due to loss of food for predator species and goes to extinction when 𝐸₁ is large.

Table 2 

Figure 4 

Change of 𝑦∗ with the change of the harvesting effort 𝐸1.

Now let 𝑚 and 𝐸₁ be fixed; then we have 𝑑𝑥^∗𝑑𝐸₂=𝛽𝑐𝑞₂𝛽𝑐−𝑎𝑑−𝑎𝐸₂𝑞₂²>0. (28) Therefore 𝑥^∗ is a strictly increasing function of 𝐸₂.

Also 𝑑𝑦^∗𝑑𝐸₂=𝛼𝑚²(𝑥^∗)²+𝑞₁𝐸₁𝑘𝑚(𝑥^∗)²+2𝑎𝑚𝛼+𝑎𝑘𝛼−𝛼−2𝛼𝑎𝑥^∗−𝑞₁𝐸₁𝑘𝑎−𝑘𝛼𝑚(𝑥^∗)²𝑑𝑥^∗𝑑𝐸₂. (29) Since 𝑑𝑥^∗/𝑑𝐸₂>0, therefore 𝑑𝑦^∗𝑑𝐸₂>0if𝛼𝑚²(𝑥^∗)²+𝑞₁𝐸₁𝑘𝑚(𝑥^∗)²+2𝑎𝑚𝛼+𝑎𝑘𝛼>𝛼+2𝑎𝛼𝑥^∗+𝑞₁𝐸₁𝑘𝑎+𝑘𝛼𝑚(𝑥^∗)²,𝑑𝑦^∗𝑑𝐸₂<0if𝛼𝑚²(𝑥^∗)²+𝑞₁𝐸₁𝑘𝑚(𝑥^∗)²+2𝑎𝑚𝛼+𝑎𝑘𝛼<𝛼+2𝑎𝛼𝑥^∗+𝑞₁𝐸₁𝑘𝑎+𝑘𝛼𝑚(𝑥^∗)². (30) Thus as the harvesting effort 𝐸₂ increases, the predator population increases when 𝐸₂ is smaller than the threshold value. But if the harvesting effort 𝐸₂ gradually increases above the threshold value, that is, as the harvesting effort 𝐸₂ becomes large enough, then the increasing amount of the harvesting effort 𝐸₂ can decrease predator population.

To construct Table 3 and Figure 5, we take 𝛼=8, 𝑎=0.5, 𝛽=4, 𝑐=0.8, 𝑑=2, 𝑘=50, 𝑞₁=0.6, 𝑞₂=0.2, 𝐸₁=0.5, and 𝑚=0.35 in appropriate units. From the analysis shown we see that as the harvesting effort 𝐸₂ increases, the prey and predator species increase. But if the harvesting effort 𝐸₂ increases larger than a threshold value, that is, as the harvesting effort 𝐸₂ becomes large enough, then increasing amount of the harvesting effort 𝐸₂ can increase the prey species but will decrease predator species and go to extinction of the predator species when 𝐸₂ is large.

Table 3 

Figure 5 

Variation of 𝑥∗ and 𝑦∗ with the change of the harvesting effort 𝐸2.

6. Numerical Simulation

As the problem is not a case study, the real-world data are not available for this model. We, therefore, take here some hypothetical data with the sole purpose of illustrating the results that we have established in the previous sections. Let us consider the parameters of the system as 𝛼=8, 𝑎=0.5, 𝛽=4, 𝑐=0.8, 𝑑=2, 𝑘=50, 𝑞₁=0.6, 𝑞₂=0.2, 𝐸₁=0.4, and 𝐸₂=0.2 in appropriate units. For these value of parameters, we get the critical value of 𝑚 as 𝑚^∗=0.398969. Thus it is easy to verify that for this set of parameters, the system (8) is locally asymptotically stable around its interior equilibrium 𝑃^∗(𝑥^∗,𝑦^∗) for 𝑚>𝑚^∗ and is unstable for 𝑚<𝑚^∗. Thus for 𝑚=𝑚^∗=0.398969 the system (8) undergoes a Hopf bifurcation. Now for 𝑚=0.75, we have interior equilibrium (0.93578,4.95175) which is asymptotically stable (see Figure 6), but 𝑚=0.2, and the interior equilibrium (0.93578,3.37539) is unstable (see Figure 7). Thus taking 𝑚 as a control parameter, it is possible to drive the prey-predator system to require equilibrium and to prevent the cycle behaviour of the system. From Figures 8, 9, 10, and 11, we see that 𝐸₁ and 𝐸₂ may also be used as controls for the system (8); Hopf bifurcation occurs when 𝐸₁=𝐸^∗₁=7.62878 (here 𝐸₂=0.2, 𝑚=0.35) and 𝐸₂=𝐸^∗₂=19.4068 (here 𝐸₁=0.5, 𝑚=0.35). For the previous values of parameters and 𝑚=0.398969, we obtained one value of 𝛿 say 𝛿=0.93578 and 𝑎(𝛿)=−0.332712<0. Thus we may conclude that the Hopf bifurcation around the interior equilibrium is supercritical and backward.

Figure 6 

Phase space trajectories corresponding to different initial levels, which shows that (0.93578,4.95175) is a global attractor.

Figure 7 

There is a stable limit cycle surrounding (0.93578,3.3754) with 𝑚=0.2.

Figure 8 

Phase space trajectories corresponding to different initial levels. Here 𝐸1=9.6, 𝐸2=0.2, and 𝑚=0.35.

Figure 9 

There is a stable limit cycle surrounding (0.93578,2.11487) with 𝐸1=6, 𝐸2=0.2, and 𝑚=0.35.

Figure 10 

Phase space trajectories corresponding to different initial levels with 𝐸1=0.5, 𝐸2=19.6, and 𝑚=0.35.

Figure 11 

There exist Hopf-bifurcating small amplitude periodic solutions with 𝐸1=0.5, 𝐸2=17, and 𝑚=0.35.

7. Concluding Remarks

This paper deals with a prey-predator model with Holling type II functional response incorporating a constant prey refuge and independent harvesting in either species. Oscillatory behavior and existence of limit cycles in harvested prey-predator system are common in nature. It is noted that constant prey refuge plays an important role in the dynamics of the proposed model system. It is also observed from the obtained results that constant prey refuge can cause an unstable equilibrium to become stable, and even a simple Hopf bifurcation occurred when the parameter 𝑚 passes through its critical value. There exists a threshold value of 𝑚, such that for the prey refuge smaller than this threshold, increasing the amount of prey refuge can increase the predator population, and if the prey refuge is larger than the threshold, increasing the amount of prey refuge can decrease the predator population. We have proved that exactly one stable limit cycle occurs when the positive equilibrium is unstable. We also determined the critical value of 𝛿 at which Hopf bifurcation occurs and observed that the bifurcation is supercritical and backward. It was also found that it is possible to control the system in such a way that the system approaches a required state, using the efforts 𝐸₁ and 𝐸₂ as controls.

Our analytical results and numerical simulation also indicate that dynamic behavior of the model not only depends on the prey refuge parameter 𝑚 but also depends on harvesting efforts 𝐸₁ and 𝐸₂. Hence it is possible to control the system in such a way that the system approaches a required state, using the harvesting efforts 𝐸₁ and 𝐸₂ or prey refuge 𝑚 as controls.

In our model, we have considered the catch-rate function based on catch-per-unit-effort hypothesis; that is, ℎ=𝑞𝐸𝑢 (𝑢 and 𝐸 denote the prey or predator population and harvesting effort, resp.). But this type of catch-rate function embodies some defects in that (i) it assumes random search for fish, (ii) it assumes equal likelihood of being captured for every fish, (iii) there is unbounded linear increase of ℎ with respect to 𝐸 for fixed 𝑢, and (iv) there is unbounded linear increase ℎ with respect to 𝑢 for a fixed 𝐸. These unrealistic features can largely be removed by adopting the alternative functional form ℎ=𝑞𝐸𝑢/(𝑛₁𝐸+𝑛₂𝑢), where 𝑛₁ and 𝑛₂ are positive constants, but we leave it for our future research work. The entire study of the paper is mainly based on the deterministic framework. On the other hand, it will be more realistic if it is possible to consider the model system in the stochastic environment due to some ecological fluctuations and other factors. Thus, a future research problem would be considered in stochastic environment.

Appendix

Detailed Calculation of Formula (23)

First we translate the equilibrium (𝛿,𝑦_𝛿) to the origin by translation 𝑥=𝑥−𝛿, 𝑦=𝑦−𝑦_𝛿. (Still denote 𝑥 and 𝑦 by 𝑥 and 𝑦, resp.). Thus the system (18) becomes 𝑑𝑥()𝑓()−𝑑𝑡=ℎ𝑥+𝛿𝑥+𝛿𝑦+𝑦_𝛿,𝑑𝑦𝑑𝑡=𝑐𝑦+𝑦_𝛿ℎ(𝑥+𝛿)−𝑑^′. (A.1) We write the system (A.1) as follows: 𝑑𝑥()𝑓()(𝛿)𝑓(𝛿)−ℎ()𝑑𝑡=ℎ𝑥+𝛿𝑥+𝛿−ℎ𝑥+𝛿𝑦+𝑦_𝛿−ℎ(𝛿)𝑦_𝛿,𝑑𝑦𝑑𝑡=−𝑐𝑑^′𝑦+𝑐ℎ(𝑥+𝛿)𝑦+𝑦_𝛿−𝑐ℎ(𝛿)𝑦_𝛿, (A.2) where 𝑓(𝛿)=𝑦_𝛿 and ℎ(𝛿)=𝑑^′.

Now compute the Taylor expansion of related functions 𝑦+𝑦_𝛿ℎ(𝑥+𝛿)−ℎ(𝛿)𝑦_𝛿=𝑎₁₀𝑥+𝑎₀₁𝑦+𝑎₂₀𝑥²+𝑎₁₁𝑥𝑦+𝑎₃₀𝑥³+𝑎₂₁𝑥²𝑦+𝑂|𝑥|⁴,|𝑥|³𝑦,ℎ(𝑥+𝛿)𝑓(𝑥+𝛿)−ℎ(𝛿)𝑓(𝛿)=𝑏₁₀𝑥+𝑏₂₀𝑥²+𝑏₃₀𝑥³+𝑂|𝑥|⁴, (A.3) where 𝑎₁₀=ℎ^′(𝛿)𝑦_𝛿,𝑎₂₀=12ℎ^′′(𝛿)𝑦_𝛿,𝑎₃₀=16ℎ^′′′(𝛿)𝑦_𝛿,𝑎₀₁=ℎ(𝛿),𝑎₁₁=ℎ^′(𝛿),𝑎₂₁=12ℎ^′′(𝛿),𝑏₁₀=𝑓ℎ^′(𝛿),𝑏₂₀=12𝑓ℎ^′′(𝛿),𝑏₃₀=16𝑓ℎ^′′′(𝛿). (A.4) Then the system (A.1) becomes 𝑑𝑥𝑑𝑡𝑑𝑦𝑥𝑦+𝐹𝑑𝑡=𝐽¹𝐹𝑥,𝑦,𝛿₂𝑥,𝑦,𝛿, (A.5) where 𝑊(𝛿)𝑌(𝛿)𝑍(𝛿)0,𝐹𝐽=₁𝑥,𝑦,𝛿=𝑏₂₀−𝑎₂₀𝑥²−𝑎₁₁𝑥𝑦+𝑏₃₀−𝑎₃₀𝑥³−𝑎₂₁𝑥²𝑦+𝑂|𝑥|⁴,|𝑥|³𝑦,𝐹₂𝑥,𝑦,𝛿=𝑐𝑎₂₀𝑥²+𝑎₁₁𝑥𝑦+𝑎₃₀𝑥³+𝑎₂₁𝑥²𝑦+𝑂|𝑥|⁴,|𝑥|³𝑦. (A.6) Now we define a matrix𝑃=10𝑁𝑀, (A.7) where 𝑁=−𝑊(𝛿)/2𝑌(𝛿) and 𝑀=−−4𝑌(𝛿)𝑍(𝛿)−𝑊²(𝛿)/2𝑌(𝛿); then 𝑃⁻¹=−𝑁10𝑀1𝑀. (A.8) By using the linear transformation 𝑥𝑦𝑢𝑣=𝑃, (A.9) we have 𝑥𝑦=𝑢𝑁𝑢+𝑀𝑣,𝑃⁻¹𝑢𝑣(𝛿)𝑢𝑣𝐽𝑃=𝐽, (A.10) where 𝐽(𝛿)=𝛽(𝛿)(𝛿)𝜔(𝛿)𝛽(𝛿)−𝜔. (A.11) Then system (A.5) becomes 𝑑𝑢𝑑𝑡𝑑𝑣(𝛿)𝑢𝑣+𝐹𝑑𝑡=𝐽¹()𝐹𝑢,𝑣,𝛿₂(𝑢,𝑣,𝛿), (A.12) where 𝛽(𝛿) and 𝜔(𝛿) are defined in (21), and 𝐹¹(𝑢,𝑣,𝛿)=𝐹₁(𝑢,𝑁𝑢+𝑀𝑣,𝛿)=𝐴₂₀𝑢²+𝐴₁₁𝑢𝑣+𝐴₃₀𝑢³+𝐴₂₁𝑢²𝑣+𝑂|𝑢|⁴,|𝑢|³|𝑣|,𝐹²()𝑁𝑢,𝑣,𝛿=−𝑀𝐹₁()+1𝑢,𝑁𝑢+𝑀𝑣,𝛿𝑀𝐹₂(𝑢,𝑁𝑢+𝑀𝑣,𝛿)=𝐵₂₀𝑢²+𝐵₁₁𝑢𝑣+𝐵₃₀𝑢³+𝐵₂₁𝑢²𝑣+𝑂|𝑢|⁴,|𝑢|³|𝑣|, (A.13) where 𝐴₂₀=𝑏₂₀−𝑎₂₀−𝑎₁₁𝑁,𝐴₁₁=−𝑎₁₁𝐴𝑀,₃₀=𝑏₃₀−𝑎₃₀−𝑎₂₁𝑁,𝐴₂₁=−𝑎₂₁𝐵𝑀,₂₀=𝑐𝑀𝑎₂₀+𝑎₁₁𝑁−𝑁𝑀𝑏₂₀−𝑎₂₀−𝑎₁₁𝑁,𝐵₁₁=𝑐𝑎₁₁+𝑎₁₁𝐵𝑁,₃₀=𝑐𝑀𝑎₃₀+𝑎₂₁𝑁−𝑁𝑀𝑏₃₀−𝑎₃₀−𝑎₂₁𝑁,𝐵₂₁=𝑐𝑎₂₁+𝑎₂₁𝑁. (A.14) Rewrite the system (A.12) in a polar coordinate form as: 𝑟̇=𝛽(𝛿)𝑟+𝑎(𝛿)𝑟³̇(𝛿)(𝛿)𝑟+⋯,𝜃=𝜔+𝑐²+⋯. (A.15) Then the Taylor expansion of (A.15) at 𝛿=𝛿 yields 𝑟̇=𝛽^′𝛿𝛿−𝛿𝑟+𝑎𝛿𝑟³+𝑂𝑟𝛿−𝛿²,𝑟³𝛿−𝛿,𝑟⁵,𝜃̇=𝜔𝛿+𝜔^′𝛿𝛿−𝛿+𝑐𝛿𝑟²+𝑂𝛿−𝛿²,𝑟²𝛿−𝛿,𝑟⁴. (A.16) In order to determine the stability of the periodic solution, we need to calculate the sign of the coefficient 𝑎(𝛿), which is given by 𝑎𝛿=1𝐹16^{1𝑢𝑢𝑢}+𝐹¹_𝑢𝑣𝑣+𝐹²_𝑢𝑢𝑣+𝐹²_𝑣𝑣𝑣+116𝜔𝛿𝐹¹_𝑢𝑣𝐹¹_𝑢𝑢+𝐹¹_𝑣𝑣−𝐹²_𝑢𝑣𝐹²_𝑢𝑢+𝐹²_𝑣𝑣−𝐹¹_𝑢𝑢𝐹²_𝑢𝑢+𝐹¹_𝑣𝑣𝐹²_𝑣𝑣, (A.17) where all partial derivatives are evaluated at the bifurcation point; that is, (𝑢,𝑣,𝛿)=(0,0,𝛿).

Since 𝑣 is linear in both 𝐹¹(𝑢,𝑣,𝛿) and 𝐹²(𝑢,𝑣,𝛿), we have that 𝐹¹_𝑢𝑣𝑣=𝐹²_𝑣𝑣𝑣=𝐹¹_𝑣𝑣=𝐹²_𝑣𝑣=0, (A.18)when𝛿=𝛿, 𝑁_𝛿=𝑁|_𝛿=𝛿=0, 𝑀_𝛿=𝑀|_𝛿=𝛿=𝑐𝑓(𝛿)ℎ^′(𝛿)/ℎ(𝛿), also 𝜔(𝛿)𝑀_𝛿=𝑐𝑓(𝛿)ℎ^′(𝛿).

At (0,0,𝛿), we have by simple calculation that 𝐹¹_𝑢𝑢𝑢+𝐹²_𝑢𝑢𝑣=6𝐴₃₀+2𝐵₂₁=6𝑏₃₀−𝑎₃₀+2𝑐𝑎₂₁,𝐹¹_𝑢𝑣𝐹¹_𝑢𝑢=2𝐴₁₁𝐴₂₀=−2𝑏₂₀−𝑎₂₀𝑎₁₁𝑀_𝛿,𝐹²_𝑢𝑣𝐹²_𝑢𝑢=2𝐵₂₀𝐵₁₁=2𝑐²𝑀_𝛿𝑎₂₀𝑎₁₁,𝐹¹_𝑢𝑢𝐹²_𝑢𝑢=4𝐴₂₀𝐵₂₀=4𝑐𝑀_𝛿𝑎₂₀𝑏₂₀−𝑎₂₀. (A.19) Now 𝑎𝛿=1𝐹16^{1𝑢𝑢𝑢}+𝐹²_𝑢𝑢𝑣+116𝜔𝛿𝐹¹_𝑢𝑣𝐹¹_𝑢𝑢−𝐹²_𝑢𝑣𝐹²_𝑢𝑢−𝐹¹_𝑢𝑢𝐹²_𝑢𝑢=16𝑏16³⁰−𝑎₃₀+2𝑐𝑎₂₁+116𝜔𝛿−2𝑏₂₀−𝑎₂₀𝑎₁₁𝑀_𝛿−2𝑐²𝑀_𝛿𝑎₁₁𝑎₂₀−4𝑐𝑀_𝛿𝑎₂₀𝑏₂₀−𝑎₂₀=116𝑓ℎ^′′′𝛿−ℎ^′′′𝛿𝑓𝛿+𝑐ℎ^′′𝛿+116𝑐𝑓𝛿ℎ^′𝛿×ℎ^′′𝛿𝑓𝛿−𝑓𝑔^′′𝛿𝑐𝑓𝛿ℎ^′𝛿²ℎ𝛿−𝑐²𝑓𝛿ℎ^′𝛿ℎ^′′𝛿−𝑐𝑓𝛿ℎ^′′𝛿𝑓ℎ^′′𝛿−ℎ^′′𝛿𝑓𝛿=116𝑓ℎ^′′′𝛿−ℎ^′′′𝛿𝑓𝛿+1ℎ16^′′𝛿𝑓𝛿−𝑓𝑔^′′𝛿×ℎ^′𝛿²+ℎ𝛿ℎ^′′𝛿ℎ𝛿ℎ^′𝛿=116ℎ^′𝛿𝑓^′′′𝛿ℎ𝛿ℎ^′𝛿+2ℎ^′𝛿²𝑓^′′𝛿−𝑓^′′𝛿ℎ^′′𝛿ℎ𝛿 (A.20) since 𝑓^′(𝛿)=0.