Principal eigenvalue for an elliptic problem with indefinite weight on cylindrical domains

This paper is concerned with an indefinite weight linear eigenvalue problem in cylindrical domains. We investigate the minimization of the positive principal eigenvalue under the constraint that the weight is bounded by a positive and a negative constant and the total weight is a fixed negative constant. Biologically, this minimization problem is motivated by the question of determining the optimal spatial arrangement of favorable and unfavorable regions for a species to survive. Both our analysis and numerical simulations for rectangular domains indicate that there exists a threshold value such that if the total weight is below this threshold value then the optimal favorable region is a circular type domain at one of the four corners, and a strip at the one end with shorter edge otherwise.


Introduction
Consider the following linear eigenvalue problem with indefinite weight ∆ϕ + λm(x)ϕ = 0 in Ω, ∂ϕ ∂n = 0 on ∂Ω, (1.1) where Ω is a bounded domain in R N with smooth boundary ∂Ω, n is the outward unit normal vector on ∂Ω, and the weight m is a bounded measurable function satisfying We say that λ is a principal eigenvalue of (1.1) if λ has a positive eigenfunction ϕ ∈ H 1 (Ω).It was shown by [3,19,13] that (1.1) has a positive principal eigenvalue if and only if the set Ω + = {x ∈ Ω : m(x) > 0} has positive Lebesgue measure and Moreover, λ is the only positive principal eigenvalue, and it is also the smallest positive eigenvalue of (1.1).
We are interested in the dependence of the principal eigenvalue λ = λ(m) on the weight m.The motivation comes from the diffusive logistic equation in Ω × R + , ∂u ∂n = 0 on ∂Ω × R + , u(x, 0) ≥ 0, u(x, 0) ≡ 0 in Ω, (1.4) where u(x, t) represents the density of a species at location x and time t, and the no-flux boundary condition means that no individuals cross the boundary of the habitat Ω, and ω is a positive parameter.The weight m represents the intrinsic growth rate of species: it is positive in the favorable part of habitat (Ω + ) and negative in the unfavorable one (Ω − = {x ∈ Ω : m(x) < 0}).The integral Ω m measures the total resources in a spatially heterogeneous environment.The logistic equation (1.4) plays an important role in studying the effects of dispersal and spatial heterogeneity in population dynamics, see, e.g.[4,6,7] and references therein.It is known that if ω ≤ λ(m), then u(x, t) → 0 uniformly in Ω as t → ∞ for all non-negative and non-trivial initial data, i.e., the species goes to extinction; if ω > λ(m), then u(x, t) → u * (x) uniformly in Ω as t → ∞, where u * is the unique positive steady solution of (1.4) in W 2,q (Ω) for every q > 1, i.e., the species survives.We are particularly concerned with the effects of spatial variation in the environment on species extinctions.In this connection, let m 0 < 1 be a positive constant and assume (A1) m satisfies (1.2), Ω + has positive measure, and Ω m ≤ −m 0 |Ω|.
Since the species can be maintained if and only if ω > λ(m), we see that the smaller λ(m) is, the more likely the species can survive.In this connection, the following question was raised and addressed by Cantrell and Cosner in [4,5]: Among all functions m(x) that satisfy (A1), which m(x) will yield the smallest principal eigenvalue λ(m)?From the biological point of view, finding such a minimizing function m(x) is equivalent to determining the optimal spatial arrangement of the favorable and unfavorable parts of the habitat for species to survive [4,5].This issue is important for public policy decisions on conservation of species with limited resources.We also refer to [2,7,8,12,15,20] and references therein for related work on spatial arrangement of resources and habitat fragmentation.
Given m 0 < 1 and κ > 0, we define and set The following result was established in [17]: In particular, Theorem A implies that the global minimizers of λ(m) in M must be of "bang-bang" type.When N = 1 and Ω is an interval, a complete characterization of all global minimizers of λ(m) in M is also given in [17].Theorem B implies that when Ω is an interval, then there are exactly two global minimizers of λ(m) (up to change of a set of measure zero).For the logistic model, this means that a single favorable region at one of the two ends of the whole habitat provides the best opportunity for the species to survive.
A major open problem is the characterization of the optimal set E in M for higherdimensional domains.In this paper, we focus on the case where Ω is a cylindrical domain given by Ω := (0, 1) where D is a bounded domain in R N −1 with smooth boundary ∂D.As we shall see in later discussions, even for this simple-looking case, determining the shape of the optimal set E is fairly non-trivial.
Let Ω + 0 and Ω − 0 be subsets of Ω defined by with a parameter c ∈ (0, 1), and set where x ∈ (0, 1) and y ∈ D. Note that (1.3) is equivalent to As we will see later, the principal eigenvalue is uniquely determined by √ κ tan Our interest is whether or not λ is locally minimal with respect to perturbation of Ω + 0 .Let µ be the smallest positive eigenvalue of For each c ∈ (0, c * ), we can define a positive number µ c ∈ (λκ, ∞) uniquely by (1.8) Now our main analytical result can be stated as follows.
In particular, the strip at the end with much longer edge can not be an optimal favorable region.
When µ > µ c , we expect that λ is locally minimal at least in a sufficiently wide class of perturbations, see Sect.6 for more details.Since µ is large when D is small, we conjecture that λ is locally minimal for a thin cylinder.It seems that λ is a global minimum for a sufficiently thin cylinder.

Remark.
1. Since µ c > λκ and λ → ∞ as c → 0, we have µ c → ∞ as c → 0. So another implication of Theorem 1.1 is that if the area of the favorable region is small, then a strip at any end of the rectangular domain can not be an optimal favorable region.
2. Since λ → 0 as c → c * , we have µ c → µ * as c → c * , where µ * is defined by 3. It seems that µ c is strictly decreasing in c.This means that if µ ≤ µ * , then λ is not locally minimal for any c ∈ (0, c * ).In other words, if Ω is a fat cylinder, then λ is not locally minimal regardless of the choice of c. See the remark at the end of Sect. 4.
Numerical simulations in Sect.7 on unit square indicate that the shape of the optimal Ω + depends crucially on the value of m 0 .More precisely, it seems that there exists a threshold value m * 0 ∈ (0, 1) such that if m 0 < m * 0 then Ω + is a stripe, and if m 0 > m * 0 then Ω + is a circular type domain at one of the four corners.For general rectangular domains, simulation shows that the optimal Ω + is either circular type domains located at one of the four corners or a strip at the one end with shorter edge.The numerical scheme developed in Sect.7 is based upon the projection gradient method, and it is general enough to handle both topological changes of Ω + and general domains with non-trivial topology.

One-dimensional problem
In this section, we summarize the results for the one-dimensional problem: It was shown in [17] that such m(x) is a global minimizer of the principal eigenvalue.By (2.1), we may write ϕ as From the continuity of ϕ and ϕ at x = c, we have For (A, B) = (0, 0), we obtain the characteristic equation (1.6).Since κc < 1 − c, the principal eigenvalue λ > 0 is uniquely determined from (1.6).

Formal expansion
We perturb the problem (1.1) as follows.Let g : D → R be any L 2 function satisfying D g(y)dy = 0, (3.1) and define a perturbed domain where ε is a small parameter.Then we set and consider the perturbed problem The main goal of this section is to find a formal asymptotic expansion of λ for positive small and gain some insight.To this end, we expand λ ε and ϕ ε formally as where (ϕ, λ) is an eigenpair of the one-dimensional problem (2.1).We substitute (3.3) into the weak form and compute ε 0 -, ε 1 -and ε 2 -order terms.
The first term in the left-hand side of (3.4) is written as follows For the second term, we have Here Hence we have Consequently, we obtain By using this and (3.5), we compare the ε 0 -, ε 1 -and ε 2 -order terms of (3.4).First, from ε 0 -order terms, we have This equality clearly holds by (1.1).
Next, from ε 1 -order terms, we have If we take ψ = ϕ, then Hence by (3.1), we obtain λ 1 = 0 so that This in particular means that ϕ 1 must satisfy the equation Finally, from ε 2 -order terms and λ 1 = 0, we have (3.12) Again by taking ψ = ϕ, we obtain Hence λ 2 must be expressed as The sign of λ 2 is crucial for stability.Indeed, we will show in Section 5 that the principal eigenvalue λ is NOT locally minimal if

Computation of λ 2
In this section, we compute λ 2 in the case where g(y) is an eigenfunction of (1.7) associated with a positive eigenvalue µ > 0. We note that (3.1) holds in this case.
[Case I: µ ≤ λκ] We write P as From the matching condition at x = c, we have By the first equality, we have Using this, we obtain .
[Case II: µ ≥ λκ] We write P as From the matching condition at x = c, we have By the first condition, we have Using this, we obtain

I[g]
Let us compute I[g].We note that I[g] is written as In Case I (µ ≤ λκ), we have Here, from the above computation, we have Recalling the characteristic equation (1.6), we see that the left-hand side satisfies On the other hand, since µ ≤ λκ, we have These inequalities show that (4.3) always holds in Case I. Thus we have shown that In Case II (µ ≥ λκ), we have Here

λ 2
Let Φ(s) be a function defined by Then (1.8) can be written as Noting that Φ is monotone increasing in s and Φ(s) → ∞ as s → ∞, we can define µ c uniquely by (1.8).Thus we have shown that Consequently, we obtain the following lemma.
The following lemma gives a lower bound of µ c .
Lemma 4.2 Let µ l be a positive number defined by Then µ c > µ l for any c ∈ (0, c * ).
Remark.It seems that the optimal lower bound of µ c is µ * , defined by (1.9).To show this, it suffices to prove the monotonicity of µ c with respect to c.

Rigorous proof by using a variational characterization
If λ 2 < 0, then we expect that λ is not locally minimal.In order to prove this rigorously, we use the following well-known variational characterization of the positive principal eigenvalue [1,3,13,19]: Lemma 5.1 The positive principal eigenvalue λ of (1.1) is given by where Moreover, λ is simple, and the infimum is attained only by associated eigenfunctions that do not change sign in Ω.
Proof of Theorem 1.1.
Let g(y) be an eigenfunction of (1.7) associated with a positive eigenvalue µ > 0, and λ ε be an eigenvalue of (3.2) for such g.
Using Lemma 5.1, we compare λ and λ ε .To this purpose, we define a functional and will show that J ε [U ] > 0 for some U ∈ S(m ε ) and small ε.In the argument in Section 3, ϕ 2 does not play any role in determining λ 2 .Hence we take ϕ 2 = 0 and choose U = ϕ + εϕ 1 as a test function, where ϕ 1 is the solution of (3.9), (3.10) and (3.11) constructed in Subsection 4.1.Since ϕ ∈ S(m), we have U ∈ S(m ε ) by continuity for sufficiently small ε.

Formal analysis in the case of µ > µ c
In this section, we show formally that if µ > µ c , then λ is locally minimal in the class of Let {V j } be an orthonormal basis which consists of eigenfunctions of (1.7): In particular, we set µ 0 = 0 and V 0 = 1/|D|.Since g is orthogonal to V 0 by (3.1), we expand g as Then ϕ 1 is computed as We compute where we used Here P j (c) + ϕ (c) < 0 if µ j > µ c .Since µ j ≥ µ > µ c for all j, we obtain I[g] < 0 and hence λ 2 > 0. Thus it is shown that λ ε > λ holds for any g = 0 and sufficiently small ε = 0.This strongly suggests that λ is locally minimal.We believe that the above formal analysis can be verified rigorously.However, even if this is the case, it does not immediately mean that λ is minimal in the class M as we assumed that the boundary of the favorable region Ω + is a graph of some function.Some further mathematical analysis will be needed to show the minimality of λ in the larger class M.

Numerical Simulations
In this section, we show the numerical approach to find the optimal configuration m(x) which minimizes the principle eigenvalue λ of (1.1) and satisfies (A1) Previously both theoretical and numerical approaches have been investigated to minimize the first eigenvalue value for positive m(x) [10,11,9,18].Here we focus on the indefinite weight m(x).Our approach is based on the projection gradient method [18].The idea is to start from an initial guess for m(x), evolve it along the gradient direction until it reaches the optimal configuration.However, the gradient direction may result in the violation of the constraint.A projection approach is used to project m(x) back to the feasible set.Furthermore, we propose a new binary update for m(x) to preserve the bang-bang structure.
At each iteration, we need to compute the principal eigenvalue and its corresponding eigenfunction.This is done by expanding ϕ in the FEM (finite element method) [14] basis, multiplying with a basis element, and integrating on the domain Ω.It yields a generalized eigenvalue equation which can be solved by Arnoldi algorithm [16].This can be implemented easily by using Matlab Partial Differential Equation Toolbox.

Gradient Descent Approach
Consider a variation in m(x) by an amount δm, which causes variations in ϕ and λ by δϕ and δλ respectively.The formula for the gradient of λ with respect to m(x) is The descent direction can be chosen as δm = aϕ 2 with a > 0 because This implies that if m(x) increases linearly w.r.t ϕ 2 , the principle eigenvalue decreases.However, this descent direction increases m(x) everywhere and results in the violation of the mass constraint.The way to resolve this is by using the Lagrange Multiplier Method.The descent direction is then modified to be where the constant b can be obtained by enforcing the constraint However, an interesting thing happens when M is increased.In Fig. (3), we choose the initial distribution having M = −1.The optimal domain of Ω + tends to stay as a strip with eigenvalue λ = 0.83 without going to the corner.It seems that there exist a threshold M * such that the optimal domain of Ω + is a strip if M > M * while the optimal domain is at the corner if M < M * .
In Table (1 The conclusion we have here is that the optimal configuration of Ω + on a rectangular domain is either a circular type domain at one of the four corners or a strip at the one end with shorter edge.

Discussions
We are interested in the minimization of the positive principal eigenvalue under the constraint that the weight is bounded by a positive and a negative constant and the total weight is a fixed negative constant.Biologically, this minimization problem is motivated by the question of determining the optimal spatial arrangement of favorable and unfavorable regions for a species to survive.This issue is important for public policy decisions on conservation of species with limited resources.
It was shown [17] that the global minimizers of λ(m) in M must be of "bang-bang" type.When N = 1 and Ω is an interval, it is further shown in [17] that there are exactly two global minimizers of λ(m) (up to change of a set of measure zero).The major open problem is the characterization of the optimal set E in M for higherdimensional domains.We show that for cylindrical domains, strip type domain is not locally minimal if µ < µ c , and numerical simulation suggests that the optimal favorable region is a circular type domain located at one of the four corners of the rectangles.Quite interestingly and in strong contrast, when µ > µ c , both our analysis and numerical simulation strongly indicate that the strip located at the one end with shorter edge is locally minimal, at least in a sufficiently wide class of perturbation.In particular, we conjecture that such strip is global minimal for a thin cylinder.
For general domains, more study need to be done.At the end of this paper, we show the simulation on the ellipse domain Ω = {x 2 + y 2 /4 = 1} with the initial guess for Ω + = {(x, y)|(x − 0.3) 2 + (y − 0.1) 2 < (0.2) 2 } is a disk.In Fig. (6a,b,c, and d), the domain of Ω + moves toward the upper boundary and the principle eigenvalue decreases dramatically.After it reaches the boundary, it moves slowly to the right until it reaches the optimal configuration shown in Fig. (6f).The final eigenvalue becomes 17.09.So far we do not know what the optimal configuration is for general domain Ω, e.g., does the optimal Ω + prefer to stay at the higher curvature regions if the boundary is smooth?We will do further study and report it in a future paper.

Figure 1 :
Figure 1: The update process for m(x) in one dimension.
), we list the eigenvalues of a quarter of an exact circle at one of the four corners and a strip at one end of the square domain Ω = [−1, 1] × [−1, 1].Since the eigenvalues are computed numerically, we only list the first few digits.It can be seen that a quarter of a circle has a smaller principle eigenvalue when M < M * ≈ −1.15.So far we do not know the threshold value theoretically.We will study this in future work.Furthermore, we show the numerical results on the rectangular domain in Fig. (4) and Fig. (5).The domains of Ω + both occupy a quarter of the full domain Ω as the first example in the square domain.However, Fig. (4c) shows that the optimal configuration stays as a strip at the left while Fig. (5i) shows that the optimal configuration stays as a strip at the top.Notice that both of them are at the one end with shorter edge.In Fig. (5), Ω + evolves to the upper left corner first at the sixth iteration with eigenvalue 0.70 and keep changing until it reaches the optimal configuration at the top with eigenvalue 0.60.The corresponding principal eigenvalue versus the number of iterations is shown in Fig. (5j).
The principal eigenvalue v.s. the number of iterations