EFFECT OF HARVESTING QUOTA AND PROTECTION ZONE IN A EFFECT OF HARVESTING QUOTA AND PROTECTION ZONE IN A REACTION-DIFFUSION MODEL ARISING FROM FISHERY REACTION-DIFFUSION MODEL ARISING FROM FISHERY MANAGEMENT MANAGEMENT

. A reaction-diﬀusion logistic population model with spatially non-homogeneous harvesting is considered. It is shown that when the intrinsic growth rate is larger than the principal eigenvalue of the protection zone, then the population is always sustainable; while in the opposite case, there exists a maximum allowable catch to avoid the population extinction. The existence of steady state solutions is also studied for both cases. The existence of an optimal harvesting pattern is also shown, and theoretical results are complemented by some numerical simulations for one-dimensional domains.


1.
Introduction.With the exponential growth of modern industry and human consumption, the stock of many renewable natural resource such as marine fishes, forest timbers has been in sharp decline in the last century because of overharvesting [18,28].In the last 50 years, marine ecology and fishery industry community has debated about the merit of setting harvesting quota and no-harvesting zones [3,17,19].Even when people agreed that setting harvesting restriction and noharvesting zones may help the resource stock to rebound, a more subtle question is the strategy and implementation of the harvesting quota and no-harvesting zones [1,2,16].
To investigate the last question, we consider the following model equation where Ω is a bounded region in R N for N ≥ 1 with C 2,α smooth boundary ∂Ω, and u(x, t) is the species population density at the location x and time t.A no-flux boundary condition is imposed so the system is closed to the exterior.The steady state equation associated with (1) is Without the harvesting, the population has a logistic growth rate with maximal intrinsic growth rate a and normalized carrying capacity K = 1.The harvesting effort is described by the term ch(x)p(u), where c > 0 is the harvesting rate, the function p(u) is a harvesting response function which satisfies Holling type II growth condition (p) p ∈ C 1 (R + ), p(0) = 0, p (u) > 0 for u ∈ [0, ∞), and lim u→∞ p(u) = 1, and the harvesting distribution function h(x) satisfies (h) h ∈ L ∞ (Ω), 0 ≤ h(x) ≤ M for x ∈ Ω and some M > 0, and Here M is the maximum harvesting density at any location x.In this harvesting model (1), the function c • h(x) resembles the harvesting quota which is often set by government agencies to prevent the extinction of renewable resource.Typically a maximum total harvesting amount per unit time c is preset, and this total amount is distributed in several geographical subregions which is achieved here by defining the distribution function h(x) such that the total amount is kept at Ideally in a unit time, the amount c is harvested from the population, but the effective harvesting rate depends on the availability of the population.Thus a harvesting response function p(u) is used.When the population is abundant at a location x, then u(x) is large and p(u(x)) is close to 1 so that the effective harvesting rate is close to ch(x); but if the value of u(x) is not high, then the effective harvesting rate is only a fraction of ch(x).In the following, for simplicity of presentation, we assume that The model (1) without diffusion was also used in studying the grazing effect to grassland [23,24], in which u(t) is the biomass and c • p(u) is the grazing by cattle.
Here we try to address the following questions which are important for the fishery or other industries relying on renewable resource: 1.For which kind of harvesting distribution function h(x), there is a maximal allowable catch c h > 0 so that the population will become extinct if c > c h regardless of initial population u 0 ?2. If a maximal allowable catch c h > 0 exists for all h satisfying (h), then which h(x) can make c h the largest?The answers to these questions are closely related to the existence and profile of a protection zone (or no-harvesting zone) Ω 0 , which is a sub-region of Ω such that h(x) = 0 for x ∈ Ω 0 .For simplicity, we assume that Ω 0 has smooth boundary.Then the effective harvesting region is Ω * = Ω/Ω 0 .Since h(x) satisfies the condition (h), then we have where |Ω| is the Lebesgue measure (area if in R 2 ) of a region Ω.Thus Ω 0 satisfies In this paper we always assume that M > |Ω| −1 .Consider an eigenvalue problem Let λ M 1 (Ω 0 ) be the principal eigenvalue of (5).Then for harvesting pattern h satisfying (h) and associated protection zone Ω 0 , our main results can be summarized as follows: 1.If the intrinsic growth rate a > λ M 1 (Ω 0 ), then for any c > 0, the system (1) possesses at least one positive steady state solution, and for any initial population u 0 (x) ≥ ( ≡)0, the population is persistent for t ∈ (0, ∞).That is, there is no maximal allowable catch.2. If the intrinsic growth rate a < λ M 1 (Ω 0 ), then a maximal allowable catch c * (h) > 0 exists so that if c > c * (h) then the population described by (1) becomes extinct as t → ∞.Moreover an optimal harvesting strategy h 0 (x) exists to maximize c * (h) among all function h(x) satisfying (h), and h 0 (x) = M for x ∈ Ω 0 while h 0 (x) = 0 for x ∈ Ω 0 .We also use the bifurcation method (see [10,11,14,15]) to study the set of positive steady state solutions of (1).The role of a protection zone in reaction-diffusion predator-prey and competition models has been considered in [7,12,13,20,25,26].In this paper we focus on a single species model while these earlier work considered protection zone in an interaction relation between species.In the results above, the optimal harvesting strategy h 0 (x) is achieved by a "bang-bang" type function which is zero in the protection zone Ω 0 , and takes maximum value outside of Ω 0 .The location of Ω 0 relative to Ω in general cannot be determined, but in one-dimensional spatial domain case, it is known that the harvesting region Ω * is connected and near the boundary [22].Similar results on optimizing eigenvalues or protection zones have also been shown in [4,21,22].In Section 5, we show some numerical simulations of the maximum allowable catch c * (h) for different choices of harvesting patterns h(x) on a one-dimensional domain Ω = (0, L), which confirms the optimality of c h when Ω * is adjacent to the boundary and also provides quantitative comparison of c * (h) for different h(x).
Our results have the following biological conservation interpretation: given a habitat Ω for the harvested species, if there is a sufficiently large protection zone Ω 0 so that a > λ M 1 (Ω 0 ), then the protection zone Ω 0 becomes a "source" for the entire habitat; even when the harvesting outside of Ω 0 is at a maximum level, the population would not become extinct.If this can be achieved, then the policy of protection zone alone can make the conservation successful and there is no need to impose a maximum allowable catch.However with the condition (h) so the condition (4) holds, then the size of the protection zone is restricted.Hence it is possible that for any region Ω 0 satisfying (4), we always have a < λ M 1 (Ω 0 ).Our result implies that in this case, having a protection zone is not sufficient for the conservation, as when the catching level c is large, the species will become extinct.But one can still design the protection zone Ω 0 in a way so that the maximal allowable catch c * (h) is the largest so the harvesting can be maximized.
We structure the remaining part of the paper this way: in Section 2, we prove the existence and nonexistence of positive steady state solutions for the cases of a > λ M 1 (Ω 0 ) and 0 < a < λ M 1 (Ω 0 ), when the harvesting strategy h(x) is fixed; and in Section 3, we use bifurcation method to consider the set of positive steady state solutions with c as parameter.The optimal spatial harvesting strategy is studied in Section 4, and numerical simulations of the maximum allowable catches are shown in Section 5.In this paper, we denote by λ 2. Positive steady state solutions.In this section we consider the existence and nonexistence of positive steady state solutions of (1) for a given h(x) satisfying (h) and the protection zone Ω 0 is fixed.Apparently u = 0 is a trivial steady state solution of (1).From the maximum principle, any nonnegative solution of (2) is either 0 or a positive solution.In the following we always consider the positive solutions of (2).For p > N , let X = {u ∈ W 2,p (Ω) : ∂u ∂n = 0 on ∂Ω}.Then the conditions on p(u) and h(x) imply that any positive solution belongs to the Sobolev space X.First when a > λ M 1 (Ω 0 ), we have the following existence result for any positive c > 0.
1.For any fixed a > λ M 1 (Ω 0 ) and any c > 0, (2) has a minimal positive solution u c,m (x) and a maximal positive solution u c,M (x) such that u c,M (x) ≥ u c,m (x). 2. If 0 ≤ c < ab 2 /M , then the positive solution of (2) is unique, which we denote by u c (x).Furthermore, {(c, u c ) : 0 < c < ab 2 /M } is a smooth curve in the space R + × X, u c is strictly decreasing in c, and u c is globally asymptotically stable for any initial condition u 0 ≥ ( ≡)0.
Proof.We prove the existence of the positive solutions by the upper and lower solution method (see for example [5,9]).When a > λ M 1 (Ω 0 ), the following auxiliary scalar equation on Ω 0 has a unique positive solution θ a (x) (see for example, [27, Theorem 2.5]).From the maximum principle, 0 < θ a (x) < 1 for x ∈ Ω 0 .Note that previous results on the existence and uniqueness of positive solution of ( 6) are usually for Dirichlet or Neumann boundary conditions, but the same proof works for the mixed boundary condition in (6).Define u(x) = 1 and We now show that u and u are upper and lower solutions and 1) is negative on Ω, and ∂u/∂n = 0 on ∂Ω, it follows that u = 1 is an upper solution of (2).We note that the equation ( 2) can be rewritten as When x ∈ Ω * , it is easy to see that 0 is a lower solution and 0 < 1 = u.When x ∈ Ω 0 , then u is a lower solution of (8) in Ω 0 , and we also have u < 1 = u.By using the "glued" local lower solution in Ω 0 ∪Ω * (see [5,Theorem 1.25]), we see that u defined in ( 7) is a lower solution of (2).Then by the upper-lower solution theory (see [5,Theorem 1.24]), there exist a minimal and a maximal positive solutions of (2) between the upper and lower solution defined above.To show the minimal and maximal solutions obtained this way are indeed minimal and maximal solutions among all solutions of (2), we observe that if (2) has a positive solution u(x), then we must have u(x) ≤ u(x) ≤ u(x) from the maximum principle.
Thus ( 2) is a logistic type equation, and the uniqueness and linear stability of the positive solution of ( 2) is well-known in that case (see [5,27]).Moreover by using implicit function theorem and the linear stability of u c , we know that {(c, u c ) : 0 < c < ab 2 /M } is a smooth curve in the space R + × X (see [27,30]).
Next we prove that u c is increasing with respect to c. Since u c is differentiable with respect to c, then du c dc satisfies Since u c is linear stable, then the inverse of the linear operator in (9) maps nonnegative functions to positive functions from the maximum principle.That implies that du c dc < 0. Finally the global stability of u c follows from the uniqueness of u c and all solutions of (1) with nonnegative initial value converge to a positive steady state (see [5,Chapter 3]).
We remark that in general, the uniqueness result of positive solution for 0 ≤ c < ab 2 /M does not hold for all c values where the existence of positive solution is known.In Section 3, we will show a bifurcation result in which the positive solutions of (2) bifurcate from the trivial solution u = 0. We will show that the bifurcation can be subcritical which implies that at least near the bifurcation point, there are multiple positive solutions of (2).
Applying the maximum principle for parabolic equations, we also have the following corollary which shows the persistence of the population when a > λ M 1 (Ω 0 ).Corollary 2.2.Assume that a > λ M 1 (Ω 0 ), then for any positive solution u(x, t) of the system (1) with the initial value u 0 (x) ≥ ( ≡)0, we have where u(x) is defined in (7).Corollary 2.2 implies that the population cannot be extinct when the maximal intrinsic growth rate a is larger than the threshold λ M 1 (Ω 0 ).In this case, no matter how large the harvesting rate is, the population will not be wiped out by the harvesting.
Next we prove that when the growth rate a is smaller than the threshold λ M 1 (Ω 0 ), a large enough harvesting rate c leads to nonexistence of positive solutions of (2).
Proof.Suppose that this is not true, then for any c > 0, (2) has at least one positive solution u c .We take a sequence of positive solutions of (2) {(c n , u n )} with c n → ∞ as n → ∞.Multiplying (2) by u n and integrating over Ω, we have then there is a subsequence (which we still denote by {u n }) convergence to some û weakly in H 1 (Ω) and strongly in L p (Ω) for any p > 1. Define Multiplying ( 10) by v n and integrating over Ω, we see that This implies that v n is bounded in H 1 (Ω), then there is a subsequence (which we still denote by {v n }) converges to some v ≥ 0 weakly in H 1 (Ω) and strongly in L p (Ω) for any p > 1.In particular we have ||v|| L 2 (Ω) = 1.Again multiplying (2) by v n and integrating over Ω, we obtain that As n → ∞, the left-hand side of ( 11) is bounded, but c n → ∞, so it is necessary that Multiplying the equation ( 12) by ϕ 1 (the principal eigenfunction for λ M 1 (Ω 0 )), multiplying the equation ( 5) by û, integrating over Ω 0 and subtracting, we obtain Since 0 < a < λ M 1 (Ω 0 ), ϕ 1 > 0, û ≥ 0 and v ≥ 0, then we must have v = 0 almost everywhere in Ω 0 .Combining with the proof in the last paragraph, we conclude that v = 0 almost everywhere in Ω, which contradicts with ||v|| L 2 (Ω) = 1.Hence when 0 < a < λ M 1 (Ω 0 ), there exists a c * > 0 such that if c > c * , then (2) has no positive solution.
Similar to Corollary 2.2, we have the following result for the dynamical behavior of solutions when 0 < a ≤ λ M 1 (Ω 0 ) and large c.Corollary 2.4.Assume that 0 < a ≤ λ M 1 (Ω 0 ) and c > c * , where c * is defined as in Theorem 2.3.Then for any given initial value u 0 ≥ ( ≡)0, the corresponding solution u(x, t) of (1) converges to 0 uniformly for x ∈ Ω as t → ∞.
3. Bifurcation analysis.In this section, we prove the existence of non-constant steady state solutions of (2) using bifurcation theory.We now set up the abstract framework for our bifurcation analysis.For p > N , let and define the set of positive solutions of (2) to be and (c, u) is a solution of (2)}.
Note the u = 0 is a trivial solution of (2) for any c > 0. We apply the local and global bifurcation theorem to obtain the following general bifurcation result.
Proof.We define a mapping F : R × X → Y by It is easy to see that We start our analysis by a standard local bifurcation argument.We fix a, b > 0 and take c as the bifurcation parameter.For any c > 0, (2) has a trivial solutions (c, 0).Bifurcation could occur along the trivial branch {(c, 0) : c > 0}.By a simple calculation, the Fréchet derivatives of F at (u, v) are given by By letting u = 0, we can find that F u (c, 0)[φ] = ∆φ + aφ − ch(x)φ/b and the bifurcation point c = c 0 for positive solutions of (2) satisfies a = λ N 1 (c 0 h(x)/b, Ω).We prove that a bifurcation point c = c 0 exists if and only if 0 < a < λ M 1 (Ω 0 ).We use a similar idea as in the proof of Theorem 2.1 of [13].Define and for any function u ∈ H 1 M (Ω 0 ), we can extend it to Ω by defining u = 0 if x ∈ Ω 0 so the extended function is in H 1 (Ω).All such functions in H 1 (Ω) form a linear subspace of H 1 (Ω), which we denote it by H 1 M (Ω).By using the variational characterization of the eigenvalue λ N 1 (c 0 h(x)/b, Ω), we have Note that a strict inequality holds in (15) since the eigenfunction associated with λ N 1 (c 0 h(x)/b, Ω) is strictly positive so it does not belong to H 1 M (Ω).This shows that 0 < a < λ M 1 (Ω 0 ) is a necessary condition for the bifurcation to occur.To show the sufficiency, we define H(c) = λ N 1 (ch(x)/b, Ω).Then clearly we have H(0) = 0 and H(c) is strictly increasing.We show that lim c→∞ where {c n } is a sequence satisfying c n → ∞, and we also assume that ||φ n || ∞ = 1.From (15) and h(x) ≥ 0, we have −∆φ n ≤ λ M 1 (Ω 0 )φ n , and Hence {φ n } is bounded in H 1 (Ω), and there is a subsequence (which we still denote by {φ n }) converging to some φ weakly in H 1 (Ω) and strongly in L 2 (Ω).Since ||φ n || ∞ = 1, we can also assume that φ n → φ in L p (Ω) for any p > 1.From Lemma 2.2 of [8], we have ||φ|| ∞ = 1.On the other hand, since {λ N 1 (c n h(x)/b, Ω)} is bounded, we may assume 16) by φ n and integrating over Ω, we obtain that As n → ∞, the left hand side of ( 18) is bounded due to (17 Since h(x) > 0 in Ω * , we must have φ(x) = 0 almost everywhere in Ω * .Since ∂Ω 0 is smooth, this implies that φ| Ω0 ∈ H 1 M (Ω 0 ).Inside Ω 0 , h(x) ≡ 0, thus by letting n → ∞ in (16), and use weak formulation, we find that φ| Ω0 is a weak non-negative solution of Then λ ∞ must equal to λ M 1 (Ω 0 ) as otherwise φ = 0 in Ω 0 and hence φ = 0 in Ω, which contradicts with φ ∞ = 1.This shows that lim c→∞ H(c) = λ M 1 (Ω 0 ) thus there exists a unique c ∈ (0, ∞) such that H(c) = a since H(c) is strictly increasing.
For the proof of part 3, the existence of the connected component Σ follows from the global bifurcation theorem in [29] or [31], and it is known that Σ is either unbounded, or connects to another (c, 0), or Σ connects to another point on the boundary of S. Since ( 2) is a logistic type equation, then the projection of Σ on the u-axis is bounded.From Theorem 2.3, Σ must be bounded, then the first alternative cannot happen.From the bifurcation analysis, (c 0 , 0) is the only bifurcation point for positive solutions, hence the second option cannot happen either.Therefore Σ must connect to (0, 1).
By using the definition of c 0 , we can further compute c (0) in ( 21) to be For a general h(x) satisfying required conditions, the sign of c (0) could be positive or negative.When c (0) > 0, there exist multiple positive solutions of ( 2) for c ∈ (c 0 , c 0 + ) where > 0. This shows that the positive solution of ( 2) is not necessarily unique, and the uniqueness result in Theorem 2.1 does not hold for all c-values.To see that c (0) is possibly positive, we take h In this special case, the positive steady state solutions bifurcating from u = 0 are indeed constant ones satisfying a|Ω|(1 − u)(b + u) = c.The multiplicity of positive steady states in this case was known as early as in [23,24].

4.
Optimal spatial harvesting strategy.In this section, we continue to address the second question mentioned in the Introduction.In Sections 2-3, we have shown that for a fixed h(x), if 0 < a < λ M 1 (Ω 0 ) (which depends on h), then there exist two threshold harvesting parameter values 0 < c 0 (h) < c * (h), such that when c > c * (h), (1) has no positive steady state solution hence the population will become extinct.Here we discuss which h will optimize c 0 (h).
We adapt the techniques of [21,22].Set For h ∈ H, define and Apparently c h 0 is the principal eigenvalue of (14).Note that the definition of λ inf is essentially independent of the function h.From the restriction (4) for Ω h 0 , the definition can be made as Clearly λ inf ≥ 0 and we conjecture that In some simpler cases (for example one-dimensional Ω) we can have a better estimate of λ inf .In the following we assume that λ inf > 0. Then we have the following result.
The following theorem asserts the existence of an optimal fishing strategy.Define where α = |Ω| − M −1 .
Proof.From Lemma 4.1, we know that c sup < ∞ when a ∈ (0, λ inf ).We show that c sup is attained.By definition, there exists a sequence {h n } ∞ n=1 in H such that c hn 0 → c sup as n → ∞.Let ϕ n > 0 be the corresponding eigenfunction of c hn 0 with sup Ω ϕ n = 1.Since c sup is finite, by the standard elliptic regularity and Sobolev embedding theorems we may assume, subject to a subsequence if necessary, that ϕ n → ϕ 0 in W 2,p (Ω) ∩ C 1,γ (Ω) as n → ∞ for any p > 1 and γ ∈ (0, 1).Clearly ϕ 0 ≥ 0 and sup Ω ϕ 0 = 1.Since {h n } is uniformly bounded, we may assume that h n → h 0 weakly in L 2 (Ω).Since h n ∈ H, then we have h 0 ∈ H.This proves that c sup is attained by h 0 ∈ H.
Hence we can find a measurable subset of Ω, say E * such that M |E * | = 1 and ϕ 0 (x) ≥ η * for x ∈ E * and ϕ 0 (x) ≤ η * for x ∈ Ω \ E * .In this case we define a function h * ∈ H by From either definition, we have h * ∈ H and Ω h * (x)dx = 1, which implies that c h * 0 ≤ c sup .We show that h 0 (x) = h * (x) almost everywhere in Ω, which implies that h 0 indeed belongs to H. First we assume that h * is defined as in (31).Then we have Since Ω h 0 ϕ 2 0 dx > 0, we have Ω h * ϕ 2 0 dx > 0. It follows that Hence c h * 0 = c sup and ϕ 0 satisfies By definition, we also have It follows that Ω c sup (h * − h 0 )ϕ 2 0 dx = 0, which implies that h * = h 0 almost everywhere in Ω since ϕ 0 > 0. The proof for the case that h * is defined as in (30) is the same.Thus h 0 ∈ H α and h 0 = M χ E for some measurable set E ⊂ Ω.Note that the proof indeed shows that h * and E * are uniquely determined by h 0 , hence the second alternative does not occur.
The subset E in Theorem 4.2 is not necessarily unique.For the location of the optimal protection zone E in higher dimensional space, there is no general result.But if the spatial dimension is one, then we can use the same techniques in [22] to obtain Theorem 4.3.Suppose that N = 1 and Ω = (0, L).Then λ inf = π 2 /(4α) 2 where α = L − M −1 , and when a ∈ (0, λ inf ), c sup is only attained by a function h ∈ H α : h = M χ E .Moreover one has either E = (0, α) or E = (1 − α, 1).
Here we point out that λ inf = λ M 1 ((0, α)) which is the principal eigenvalue of ϕ (x) + λϕ(x) = 0, x ∈ (0, α), ϕ (0) = 0, ϕ(α) = 0. (32) We remark that our results in this section only optimize c 0 (h) = c h 0 which is the principal eigenvalue of a linear eigenvalue problem, and the methods here cannot be used to optimize c * (h) which can be viewed as an eigenvalue of a nonlinear eigenvalue problem.In [21], a similar optimization was achieved for a nonlinear problem via a minimization of the energy functional.But that technique does not work here as for (2), multiple positive solutions may exist and a solution is not necessarily an energy minimizer.However the numerical simulations in Section 5 indicate that at least in the one-dimensional case, the same h 0 in Theorem 4.3 also optimizes c * (h). 5. Numerical simulations.In this section we illustrate our theoretical results with some numerical simulations in Matlab.We use an implicit finite difference scheme for (1) with one-dimensional spatial domain: For a fixed harvesting pattern h(x) defined on Ω = [0, L] satisfying 0 ≤ h(x) ≤ M and L 0 h(x)dx = 1, we define the maximal allowable catch C * 1 (h) to be u(x, t) > 0 = sup{c > 0 : there exists a positive steady state solution}. (34) When c > C * 1 (h), the population becomes extinct as t → ∞, regardless of initial population u 0 .On the other hand, when c is near (but less than) C * 1 (h), the population could still be vulnerable to a large stochastic perturbation.We define another threshold value C * 2 (h) so that when c ≤ C * 2 (h), the maximal population density exceeds half of the carrying capacity as t → ∞.That is Apparently 0 < C * 2 (h) ≤ C * 1 (h) when both exist.Figure 1 shows evolution behavior of solutions of (33) when c = C * 1 (h) (left) and c = C * 2 (h) (right) for the harvesting zone being in the middle of the domain (0, L).In numerical simulations we use two types of harvesting functions which harvest in the interior of (0, L) or near the boundary of (0, L).To be more precise, we define We also compare the result with a unform harvesting function h 0 (x) = 1/L for x ∈ (0, L).Comparing the protection zone Ω h 0 for these harvesting functions, we find that Ω h0 0 = ∅, and  Figure 2 shows the plots of C * 1 (h) and C * 2 (h) versus the parameter a with harvesting functions h(x) = h * i (x) for 1 ≤ i ≤ 4. For each of i = 1, 2, 3, the curve a → C * 1 (h, a) is increasing in a and it has a = π 2 /25 ≈ 0.3947 as a vertical asymptote, which confirms the theoretical results shown in Sections 2 and 3. On the other hand, for i = 4, the curve a → C * 1 (h, a) is also increasing in a, but the vertical asymptote is a = π 2 /100 ≈ 0.0987.This confirms the result in Section 4 that the boundary protection zone is better than the interior one.Moreover Figure 2 indicates that in general, we have C * i (h * 3 ) > C * i (h * 2 ) > C * i (h * 1 ) for i = 1, 2, which implies that for the interior protection zones, the one closer to the boundary is better.Figure 3 shows the plots of C * 1 (h) and C * 2 (h) versus the parameter a with harvesting functions h(x) = h i (x) for 1 ≤ i ≤ 4. For i = 1, 2, 3, the protection zone Ω h 0 is not connected, hence the theory which we developed in Sections 2-4 does not apply to this case.One can observe that the asymptote for a → C * 1 (h 3 , a) for h 3 (x) is different from the previous theoretical values.On the other hand, we still have

D 1 (
φ, O) and λ N 1 (φ, O) the principal eigenvalues of −∆ + φ over the bounded domain O, with homogeneous Dirichlet and Neumann boundary conditions respectively.We usually omit O in the notation if the region O = Ω.If the potential function φ = 0, we simply denote them by λ D 1 (O) and λ N 1 (O).

Figure 2 .
Figure 2. Plot of C * 1 (h) (left) and C * 2 (h) (right) versus parameter a for (33) with protection zone in the interior of Ω.Here the harvesting functions are h 0 (x) and h(x) = h * i (x) for 1 ≤ i ≤ 4, L = 10 and b = 1.The values of C * 1 (h) and C * 2 (h) are numerically estimated with each increment of a by 0.025.

Figure 3 .
Figure 3. Plot of C * 1 (h) (left) and C * 2 (h) (right) versus parameter a for (33) with harvesting zone in the interior of Ω.Here the harvesting functions are h 0 (x) and h(x) = h i (x) for 1 ≤ i ≤ 4, L = 10 and b = 1.The values of C * 1 (h) and C * 2 (h) are numerically estimated with each increment of a by 0.025.