Solving Ratio-Dependent Predator-Prey System with Constant Effort Harvesting Using Homotopy Perturbation Method

Abstract

Due to wide range of interest in use of bioeconomic models to gain insight into the scientific management of renewable resources like fisheries and forestry, homotopy perturbation method is employed to approximate the solution of the ratio-dependent predator-prey system with constant effort prey harvesting. The results are compared with the results obtained by Adomian decomposition method. The results show that, in new model, there are less computations needed in comparison to Adomian decomposition method.

1. Introduction

Partial differential equations which arise in real-world physical problems are often too complicated to be solved exactly, and even if an exact solution is obtainable, the required calculations may be practically too complicated, or it might be difficult to interpret the outcome. Very recently, some promising approximate analytical solutions are proposed such as Exp-function method, Adomian decomposition method (ADM), variational iteration method (VIM), and homotopy perturbation method (HPM).

HPM is the most effective and convenient method for both linear and nonlinear equations. This method does not depend on a small parameter. Using homotopy technique in topology, a homotopy is constructed with an embedding parameter 𝑝∈[0,1], which is considered as a “small parameter.” HPM has been shown to effectively, easily, and accurately solve a large class of linear and nonlinear problems with components converging to accurate solutions. HPM was first proposed by He [1–7] and was successfully applied to various engineering problems.

The motivation of this paper is to extend the homotopy perturbation method (HPM) [8–17] to solve the ratio-dependent predator-prey system. The results of HPM are compared with those obtained by the ADM [18]. Different from ADM, where specific algorithms are usually used to determine the Adomian polynomials, HPM handles linear and nonlinear problems in simple manner by deforming a difficult problem into a simple one. The HPM is useful to obtain exact and approximate solutions of linear and nonlinear differential equations.

In this paper, we assume that the predator in model is not of commercial importance. The prey is subjected to constant effort harvesting with 𝑟, a parameter that measures the effort being spent by a harvesting agency. The harvesting activity does not affect the predator population directly. It is obvious that the harvesting activity does reduce the predator population indirectly by reducing the availability of the prey to the predator. Adopting a simple logistic growth for prey population with 𝑒>0,𝑏>0, and 𝑐>0 standing for the predator death rate, capturing rate, and conversion rate, respectively, we formulate the problem as𝑑𝑥𝑑𝑑𝑦𝑡=𝑥(1−𝑥)−𝑦𝑏+𝑥𝑦𝑥−𝑟𝑥,𝑑𝑡=𝑦𝑐+𝑥𝑦𝑥−𝑒𝑦,(1.1) where 𝑥(𝑡) and 𝑦(𝑡) represent the fractions of population densities for prey and predator at time 𝑡, respectively. Equations (1.1) are to be solved according to biologically meaningful initial conditions 𝑥(0)≥0 and 𝑦(0)≥0 [18].

2. Applications

In this section, we will apply the HPM to nonlinear differential system of ratio-dependant predator-prey,𝐻𝜈,𝑝=1−𝑝𝐿𝜈−𝐿𝑢₀+𝑝𝐴𝜈−𝑓𝑟=0,𝑝∈0,1,𝑟𝜀Ω,(2.1)where 𝐴(𝜈) is a general differential operator which can be divided into a linear part 𝐿(𝜈) and a nonlinear part 𝑁(𝜈) and 𝑓(𝑟) is a known analytical function. 𝑝∈[0,1] is an embedding parameter, while 𝑢₀ is an initial approximation of the equation which should be solved, and satisfies the boundary conditions.

According to the HPM (relation (2.1)), we can construct a homotopy of system as follows: 11−−𝑝𝑝×𝜈²𝜈𝜈̇₁₂+𝜈₁𝜈̇₁−𝑥̇₀y₀−𝑥̇₀𝑥₀+𝑝𝜈²𝜈̇₁+𝜈₁𝜈̇₁−1−𝑏−𝑟𝜈₁𝜈₂+𝜈₂𝜈²¹−1−𝑟𝜈²¹+𝜈³¹=0,𝜈̇₂+𝜈₁𝑣̇₂−ẏ₀y₀−x₀ẏ₀+𝑝𝜈²𝜈̇₂+𝜈₁𝜈̇₂+𝑒−𝑐𝜈₁𝜈₂+e𝜈²²=0,(2.2) where dot denotes differentiation with respect to 𝑡, and the initial approximations are as follows:𝑣𝑣_1,0(𝑡)=𝑥₀(𝑡)=𝑥(0),_2,0(𝑡)=𝑦₀(𝑡)=𝑦(0).(2.3) Assume that the solution of (2.2) can be written as a power series in 𝑝 as follows:𝜈𝜈₁=𝜈_1,0+𝑝𝜈_1,1+𝑝²𝜈_1,2+𝑝³𝜈_1,3+⋯,₂=𝜈_2,0+𝑝𝜈_2,1+𝑝²𝜈_2,2+𝑝³𝜈_2,3+⋯,(2.4) where 𝜈_𝑖,𝑗(𝑖,𝑗=1,2,3,…) are functions yet to be determined. Substituting (2.3) and (2.4) into (2.2), and arranging the coefficients of p powers, we have𝑣^2,0𝑣̇+^1,0+𝑣𝑣_1,031,0𝑣̇_1,0+𝑣_1,1−𝑣^21,0+𝑣_1,0𝑣̇_1,1+𝑣_2,0𝑣̇_1,1+𝑟𝑣_1,0𝑣_2,0+𝑏𝑣_1,0𝑣_2,0−𝑣_1,0𝑣_2,0+𝑣_2,0𝑣^21,0+𝑟𝑣^21,0𝑝𝑣̇_1,1+𝑣_1,0𝑣̇_1,2+𝑣_2,0𝑣̇_1,2+𝑣_2,1𝑣̇_1,1+2𝑟𝑣_1,0𝑣_1,1+𝑏𝑣_1,0𝑣_2,1+2𝑣_2,0𝑣_1,0𝑣_1,1+𝑟𝑣_1,1𝑣_2,0+𝑟𝑣_1,0𝑣_2,1+𝑏𝑣_1,1𝑣_2,0−𝑣_1,0𝑣_2,1−𝑣_1,1𝑣_2,0+𝑣_2,1𝑣^21,0−2𝑣_1,0𝑣_1,1+3𝑣^21,0𝑣_1,1𝑝²+𝑣_1,1𝑣̇_1,2+𝑣_1,2𝑣̇_1,1+𝑣_1,0𝑣̇_1,3+𝑣_2,1𝑣̇_1,2+𝑣_2,0𝑣̇_1,3+𝑣_2,2𝑣̇_1,1+𝑣_2,0𝑣^21,1−𝑣_1,0𝑣_2,2−𝑣_1,2𝑣_2,0−𝑣_1,1𝑣_2,1+𝑣_2,2𝑣^21,0+𝑟𝑣^21,1+3𝑣_1,0𝑣^21,1−𝑣^21,1+𝑏𝑣_1,1𝑣_2,1+𝑏𝑣_1,0𝑣_2,2+𝑏𝑣_1,2𝑣_2,0+𝑟𝑣_1,0𝑣_2,2+𝑟𝑣_1,1𝑣_2,1+𝑟𝑣_1,2𝑣_2,0+2𝑣_2,0𝑣_1,0𝑣_1,2+2𝑟𝑣_1,0𝑣_1,2+2𝑣_2,1𝑣_1,0𝑣_1,1+3𝑣^21,0𝑣_1,2𝑣−2𝑣_1,0𝑣_1,2𝑝³+⋯=0,_2,0𝑣̇+^2,0+𝑒𝑣𝑣_1,01,0𝑣̇_2,0+𝑣_2,1𝑣_2,0−𝑐𝑣_1,0𝑣_2,0+𝑣_2,0𝑣̇_2,1+𝑣_1,0𝑣̇_2,1+𝑒𝑣^22,0𝑝𝑣̇_2,1+𝑒𝑣_1,0𝑣_2,1−𝑐𝑣_1,0𝑣_2,1+𝑒𝑣_1,1𝑣_2,0−𝑐𝑣_1,1𝑣_2,0+2𝑒𝑣_2,0𝑣_2,1+𝑣_2,0𝑣̇_2,2+𝑣_1,1𝑣̇_2,1+𝑣_1,0𝑣̇_2,2𝑝²+𝑒𝑣^22,1+𝑣_2,1𝑣̇_2,2+𝑣_2,2𝑣̇_2,1+𝑣_2,0𝑣̇_2,3+𝑣_1,1𝑣̇_2,2+𝑣_1,2𝑣̇_2,1+𝑣_1,0𝑣̇_2,3+𝑒𝑣_1,0𝑣_2,2+𝑒𝑣_1,1𝑣_2,1−𝑐𝑣_1,0𝑣_2,2−𝑐𝑣_1,1𝑣_2,1+𝑒𝑣_1,2𝑣_2,0−𝑐𝑣_1,2𝑣_2,0+2𝑒𝑣_2,0𝑣_2,2𝑝³+⋯=0.(2.5) In order to obtain the unknown of 𝜈_𝑖,𝑗(𝑥,𝑡),𝑖,𝑗=1,2,3,…, we must construct and solve the following system which includes 6 equations, considering the initial conditions of 𝜈_𝑖,𝑗(0)=0,𝑖,𝑗=1,2,3,… :𝑣𝑣_2,0𝑣̇_1,0+𝑣_1,0𝑣̇_1,0=0,_31,0𝑣−𝑣^21,0+𝑣_1,0𝑣̇_1,1+𝑣_2,0𝑣̇_1,1+𝑣_1,0𝑣_2,0+𝑏𝑣_1,0𝑣_2,0−𝑣_1,0𝑣_2,0+𝑣_2,0𝑣^21,0+𝑟𝑣^21,0=0,_1,1𝑣̇_1,1+𝑣_1,0𝑣̇_1,2+𝑣_2,0𝑣̇_1,2+𝑣_2,1𝑣̇_1,1+2𝑟𝑣_1,0𝑣_1,1+𝑏𝑣_1,0𝑣_2,1+2𝑣_2,0𝑣_1,0𝑣_1,1+𝑟𝑣_1,1𝑣_2,0𝑣+𝑟𝑣_1,0𝑣_2,1+𝑏𝑣_1,1𝑣_2,0−𝑣_1,0𝑣_2,1−𝑣_1,1𝑣_2,0+𝑣_2,1𝑣^21,0−2𝑣_1,0𝑣_1,1+3𝑣^21,0𝑣_1,1=0,_2,0𝑒𝑣𝑣_1,0̇_2,0+𝑣_1,0𝑣̇_2,0=0,𝑣𝑣_2,0−𝑐𝑣_1,0𝑣_2,0+𝑣_2,0𝑣̇_2,1+𝑣_1,0𝑣̇_2,1+𝑒𝑣^22,0=0,_2,1𝑣̇_2,1+𝑒𝑣_1,0𝑣_2,1−𝑐𝑣_1,0𝑣_2,1+𝑒𝑣_1,1𝑣_2,0−𝑐𝑣_1,1𝑣_2,0+2𝑒𝑣_2,0𝑣_2,1+𝑣_2,0𝑣̇_2,2+𝑣_1,1𝑣̇_2,1+𝑣_1,0𝑣̇_2,2=0.(2.6) From (2.4), if the first three approximations are sufficient, then setting 𝑝=1 yields the approximate solution of (1.1) to𝑥(𝑡)=lim_𝑝→1𝑣₁(𝑡)=^𝑘=3_𝑘=0𝑦(𝑡)=lim_𝑝→1𝑣_1,𝑘(𝑡),𝑣₂(𝑡)=^𝑘=3_𝑘=0𝑣_2,𝑘(𝑡).(2.7) Therefore,𝑣v_1,0(𝑡)=𝑥₀(𝑡)=𝑥(0),(2.8)_1,1(𝑡)=−x⁰𝑥²⁰−𝑥₀−𝑦₀+𝑥₀𝑦₀+𝑟𝑦₀+𝑏𝑦₀+𝑟𝑥₀𝑡𝑥₀+𝑦₀𝑣,(2.9)_1,2(𝑡)=12𝑥₀+𝑦₀³𝑥⁰𝑡²3𝑦⁰𝑥²⁰−𝑥²⁰𝑏𝑦₀+2𝑥³⁰𝑏𝑦₀+3𝑥⁴⁰𝑟+6𝑥³⁰𝑦²⁰−3𝑦³⁰𝑥₀+𝑥³⁰𝑟²−9𝑥³⁰𝑦₀+6𝑥⁴⁰𝑦₀−9𝑥²⁰𝑦²⁰+2𝑦³⁰𝑥²⁰−2𝑥³⁰𝑟−2𝑟𝑦³⁰−2𝑏𝑦³⁰+𝑏²𝑦³⁰+𝑟²𝑦³⁰+𝑥²⁰𝑏𝑦₀𝑟+3𝑥₀𝑟𝑦²⁰𝑏+𝑦²⁰𝑥₀𝑒𝑏+𝑏𝑥²⁰𝑦₀𝑒−𝑏𝑥²⁰𝑦₀𝑐−3𝑥₀𝑏𝑦²⁰++93𝑥𝑥₀₃₀𝑦²⁰+3𝑦³⁰𝑥₀𝑟+3𝑦³⁰𝑥₀𝑏−6𝑥²⁰𝑟𝑦₀+2𝑥⁵⁰−3𝑥⁴⁰+𝑦³⁰+2𝑟𝑦³⁰𝑏𝑣𝑟𝑦₀−6𝑥₀𝑟𝑦²⁰+9𝑥²⁰𝑦²⁰𝑟+5𝑥²⁰𝑦²⁰𝑏+𝑥³⁰+3𝑥₀𝑟²𝑦²⁰+3𝑥²⁰𝑟²𝑦₀,(2.10)_2,0𝑣(𝑡)=𝑦₀(𝑡)=𝑦(0),(2.11)_2,1(𝑡)=y⁰−𝑒𝑥⁰+𝑐𝑥₀−𝑒𝑦₀𝑡𝑦₀+𝑥₀𝑣,(2.12)_2,2(𝑡)=−12𝑦₀+𝑥₀³+𝑦⁰𝑐𝑡𝑥²₀3𝑦⁰𝑒𝑥²⁰𝑐+𝑦²⁰𝑐𝑥₀𝑒+2𝑒𝑥³⁰𝑐−𝑐𝑥²⁰𝑦₀−𝑐𝑥₀𝑦²⁰−𝑐²𝑥³⁰+𝑐𝑥³⁰𝑦₀+𝑐𝑥²⁰𝑦₀𝑟𝑦²⁰𝑏+𝑐𝑥²⁰𝑦²⁰+𝑐𝑥₀𝑦²⁰𝑟−𝑒²𝑥³⁰−3𝑦₀𝑒²𝑥²⁰−3𝑦²⁰𝑒²𝑥₀−𝑦³⁰𝑒².(2.13) We also obtained 𝑣_1,3 and 𝑣_2,3, but because they were too long to maintain, we skip them and only use them in the final numerical results. In this manner, the other components can be easily obtained by substituting (2.8) through (2.13) into (2.7) as follows:𝑥(𝑡)=𝑥(0)−x⁰𝑥²⁰−𝑥₀−𝑦₀+𝑥₀𝑦₀+𝑟𝑦₀+𝑏𝑦₀+𝑟𝑥₀𝑡+1𝑥₀+𝑦₀2𝑥₀+𝑦₀³𝑥₀𝑡²(3𝑦₀𝑥²⁰−𝑥²⁰𝑏𝑦₀+2𝑥³⁰𝑏𝑦₀+3𝑥⁴⁰𝑟+6𝑥³⁰𝑦²⁰−3𝑦³⁰𝑥₀+𝑥³⁰𝑟²−9𝑥³⁰𝑦₀+6𝑥⁴⁰𝑦₀−9𝑥²⁰𝑦²⁰+2𝑦³⁰𝑥²⁰−2𝑥³⁰𝑟−2𝑟𝑦³⁰−2𝑏𝑦³⁰+𝑏²𝑦³⁰+𝑟²𝑦³⁰+𝑥²⁰𝑏𝑦₀𝑟+3𝑥₀𝑟𝑦²⁰𝑏+𝑦²⁰𝑥₀𝑒𝑏+𝑏𝑥²⁰𝑦₀𝑒−𝑏𝑥²⁰𝑦₀𝑐−3𝑥₀𝑏𝑦²⁰+3𝑥₀𝑦²⁰+3𝑦³⁰𝑥₀𝑟+3𝑦³⁰𝑥₀𝑏−6𝑥²⁰𝑟𝑦₀+2𝑥⁵⁰−3𝑥⁴⁰+𝑦³⁰+2𝑟𝑦³⁰𝑏+9𝑥³⁰𝑟𝑦₀−6𝑥₀𝑟𝑦²⁰𝑦(𝑡)=𝑦(0)+y₀+9𝑥²⁰𝑦²⁰𝑟+5𝑥²⁰𝑦²⁰𝑏+𝑥³⁰+3𝑥₀𝑟²𝑦²⁰+3𝑥²⁰𝑟²𝑦₀+𝑣^1,3⋯,−𝑒𝑥⁰+𝑐𝑥₀−𝑒𝑦₀𝑡𝑦₀+𝑥₀−12𝑦₀+𝑥₀³×𝑦₀𝑡²3+𝑦⁰𝑐𝑒𝑥𝑥²⁰₀𝑐+𝑦²⁰𝑐𝑥₀𝑒+2𝑒𝑥³⁰𝑐−𝑐𝑥²⁰𝑦₀−𝑐𝑥₀𝑦²⁰−𝑐²𝑥³⁰+𝑐𝑥³⁰𝑦₀+𝑐𝑥²⁰𝑦₀𝑟𝑦²⁰𝑏+𝑐𝑥²⁰𝑦²⁰+𝑐𝑥₀𝑦²⁰𝑟−𝑒²𝑥³⁰−3𝑦₀𝑒²𝑥²⁰−3𝑦²⁰𝑒²𝑥₀−𝑦³⁰𝑒²+𝑣^2,3⋯.(2.14)

3. Numerical Results and Comparison with ADM

For comparison with the results obtained by ADM [18], the parameter values in four cases are considered in Table 1.

Table 1 

Parameter values used for illustration purposes.

Results of four terms approximation for 𝑥(𝑡),𝑦(𝑡) obtained by using HPM and ADM [18] are presented in (3.1), respectively: Case1∶𝑥𝑦≈≈00..35−−00..13152𝑡5+𝑡+0.109.0417868𝑡²08−𝑡₂0.107288𝑡³,Case2∶𝑥≈0.5+0.05𝑡+0.012265𝑡₂−0.0011284𝑡³,𝑦≈0.3−0.1125𝑡+0.024433−𝑡₂0.0016032𝑡³,Case3∶𝑥≈0.3+0.0799t+0.00533t₂−0.00398199𝑡³,𝑦≈0.6−0.08𝑡+0.01866𝑡₂−0.00115𝑡³,Case4∶𝑥≈0.5+0.07857𝑡−0.0160−200𝑡.₂00231𝑡³,𝑦≈0.2+0.051428𝑡+0.0055918𝑡−₂0.00119873𝑡³,Case1∶𝑥≈0.5−0.35000𝑡+0.19476𝑡₂+0.00002245𝑡³,𝑦≈0.3−0.11250𝑡+0.018809𝑡₂−0.10728𝑡³,Case2∶𝑥≈0.5+0.05000𝑡+0.012266𝑡₂−0.0011286𝑡³,𝑦≈0.3−0.11250𝑡+0.024434𝑡₂−0.0016034𝑡³,Case3∶𝑥≈0.3+0.08000t+0.005333t₂−0.0039821𝑡³,𝑦≈0.6−0.08000𝑡+0.018667𝑡₂−0.0011555𝑡³,Case4∶𝑥≈0.5+0.07857𝑡−0.016021𝑡₂−0.0023112𝑡³,𝑦≈0.2+0.051430𝑡+0.0055920𝑡−₂0.0011984𝑡³,+0.00002246𝑡³.(3.1) Figures 1–4 show the relations between prey and predator populations versus time.

Figure 1 

Population fraction versus time for Case 1: 𝑟=0.9: (—) prey population fraction; (◦ ◦ ◦) predator population fraction.

Figure 2 

Population fraction versus time for Case 2: 𝑟=0.1: (—) prey population fraction; (◦ ◦ ◦) predator population fraction.

Figure 3 

Population fraction versus time for Case 3: 𝑟=0.1: (—) prey population fraction; (◦ ◦ ◦) predator population fraction.

Figure 4 

Population fraction versus time for Case 4: 𝑟=0.2: (—) prey population fraction; (◦ ◦ ◦) predator population fraction.

A noteworthy observation from Figure 1 is that prey and predator species can become extinct simultaneously for some values of parameters, regardless of the initial values. Thus, overexploitation of the prey population by constant effort harvesting process together with high predator capturing rate may lead to mutual extinction as a possible outcome of predator-pray interaction. In Figure 2, only the predator population gradually decreases and becomes extinct despite the availability of increasing prey population. This can be attributed to the effect of the predator death rate, being greater than the conversion rate and low constant prey harvesting as shown in Case 2 (see Table 1). Figures 3 and 4 illustrate the possibility of predator and prey long-term coexistence. Depending on the initial values, both prey and predator populations increase or reduce in order to allow long-term coexistence [18].

4. Conclusion

Homotopy perturbation method was employed to approximate the solution of the ratio-dependent predator-prey system with constant effort prey harvesting. The results obtained here were compared with results of Adomian decomposition method. The results show that there is less computations needed in comparison to ADM.