We consider the problem of minimizing the kth eigenvalue of rectangles with unit area and Dirichlet boundary conditions. This problem corresponds to finding the ellipse centred at the origin with axes on the horizontal and vertical axes with the smallest area containing k integer lattice points in the first quadrant. We show that, as k goes to infinity, the optimal rectangle approaches the square and, correspondingly, the optimal ellipse approaches the circle. We also provide a computational method for determining optimal rectangles for any k and relate the rate of convergence to the square with the conjectured error term for Gauss's circle problem.
The optimization of eigenvalues of the Laplace operator among domains of equal volume is a classical problem in spectral theory that dates at least as far back as the work of Rayleigh at the end of the nineteenth century . In spite of the fact that a lot of progress has been made on this problem, it was only recently, for instance, that existence among quasi-open sets was shown for the general kth eigenvalue [2,3].
Furthermore, most of the work done so far in terms of identifying the optimal shapes for different boundary conditions, namely Dirichlet, Neumann and Robin, is related to the first two frequencies and has spanned more than a century [1,4–12]. In this case, the optimizers are one ball and two balls of equal volume for the first and second eigenvalues, respectively. Recently, other related optimization problems with different restrictions or boundary conditions have also started appearing in the literature (e.g. [13–16]).
Within the last few years, approaches based on a combination of numerical and analytical methods have shown that it is not to be expected that optimizers of higher eigenvalues will, in general, have an analytical description in closed form in terms of elementary functions [17–19]. Furthermore, some expected properties such as the existence of symmetry for optimizers might also fail .
On the other hand, one may consider what happens at the other end of the spectrum, namely for the high frequencies, and see if it is possible to uncover any structure there. The purpose of this paper is to give a first step in this direction by studying the particular case of optimization over rectangles.
A first remark is in order here—namely, this should not be dismissed as a mere toy problem. In fact, and as we recall in §3c, eigenvalues of rectangles are related to the integer lattice point problem on ellipses, a classical problem in analytic number theory dating at least back to Gauss . More precisely, minimizing the kth eigenvalue over rectangles of unit area with Dirichlet boundary conditions is equivalent to finding the ellipse of minimal area which encloses k integer lattice points in the first quadrant, excluding the axes, but possibly including points on the boundary of the ellipse. As far as we are aware, this is the first time that a problem of this type is being addressed, in either the spectral or the number theory literature—see also the article by Bucur & Freitas , where this problem is studied for different classes of domains in two dimensions but with a perimeter restriction; in this case, and because this automatically ensures the necessary uniform condition on optimizers, it is possible to show that these converge to the domain in the given class having the largest area.
In order to get an idea for what is at stake here, let us begin by looking at the asymptotic behaviour of eigenvalues for a specific domain. The two-term Weyl asymptotics for the eigenvalues of the Dirichlet Laplacian on a two-dimensional domain Ω in is given by 1.1where |Ω| and |∂Ω| denote, respectively, the area and perimeter of Ω .
Naively, we might expect from this asymptotic behaviour that, in order to minimize λk among domains of equal area, we would want to minimize the second term in the asymptotics, which would correspond to minimizing the perimeter. Thus, in the case of rectangles, we would get the square as the minimizer, and for general domains, the disc.
It is clear, however, that an argument just based on asymptotics is not correct, as the dependence of the remaining term on quantities such as the perimeter is not known, and we do not have any control over this behaviour, which is uniform in all possible domains.
Furthermore, it has been shown recently that, in general, this expected result will not actually hold . Indeed, in the case of Robin boundary conditions, not only is the optimal asymptotic shape, if it exists, not the disc but, as was shown in that paper, the sequence of optimal domains will not satisfy the Weyl asymptotics corresponding to a single domain—more precisely, and while the leading term in the Weyl asymptotics in dimension d is of order k2/d, the optimal sequence will grow at most with k1/d in this case.
There is, however, a main difference between Robin and Dirichlet boundary conditions in this respect. More precisely, eigenvalues of the latter satisfy a lower bound which is a multiple of the leading term in (1.1)—in fact, the famous Pólya conjecture states that this lower bound should be exactly the leading term .
Although Pólya's conjecture has only been proved for tiling domains so far, lower bounds of a similar type are known to be true, although not with the optimal constant (e.g. [24,25]). In the two-dimensional case, this bound becomes and, although possibly not optimal, this is enough to ensure that the leading term in the asymptotics of any optimal sequence must also be of order k. For the special case of rectangles, and because these are tiling domains, we actually have
The main result of the paper, theorem 2.1 in §2, states that, for rectangles, as k approaches infinity the optimal shape approaches the square or that, equivalently, the optimal ellipse described earlier converges to the circle.
The crucial ingredient in the proof is to ensure that the sequence of optimal rectangles remains bounded as k goes to infinity, as once this is done we have reduced our problem to minimizing over a compact set and we may then use the known estimates for the integer lattice point problem for planar convex domains. To achieve this step, we had to derive an improved Pólya bound, as the available results for the integer lattice problem either just asserted the existence of a constant without explicitly indicating its dependence on the perimeter of the rectangle or the constants involved were too large for our purposes—this is a consequence of the fact that in integer lattice problems the emphasis is put on the order of the error terms, while here we do not need the best results in that direction.
Another aspect related to this problem that we consider in the second part of the paper is the issue of actually determining the kth optimal domain explicitly. As can be seen in §4c, clearly the branching of this function grows very rapidly, making it a very complex problem which we believe to be both interesting and challenging in its own right. To this end, we have devised a method to determine optimal domains computationally. This is based on the interpretation of the problem of calculating the kth Dirichlet eigenvalue of a given rectangle with unit area as that of determining the integer lattice points between two particular ellipses. This approach reduces significantly the number of eigenvalues that need to be calculated at each step. The optimization problem itself is solved by an algorithm involving the Nelder–Mead simplex method.
The structure of the paper is as follows. In §2, we introduce some notation and state the main result. The proof is then included in §3, which is divided into subsections as described in the main steps above. The computational part occupies most of the remaining part of the paper. Section 4 contains the description of the algorithm used with the following section presenting the results obtained. Finally, in the last section, we discuss the results obtained.
2. Notation and main result
Let denote a rectangle with sides a and 1/a (a≥1) and denote its Dirichlet eigenvalues by λk=λk(a), k=1,2,… . These are given by
We are interested in the problem of minimizing λk(a) as a function of a for all k, that is, 2.1
The above problem determines a sequence of values of a for which the minimizer is attained, which we shall denote by . While the first three or four optimal rectangles are not difficult to compute explicitly, the task becomes rather complex as k increases and the number of branchings that it is necessary to examine increases quite fast. In fact, the function to be minimized develops a large number of points where it is not differentiable, resulting from the intersections of eigencurves corresponding to different rectangles—see §4c, for instance. The first 15 optimal values of are found to be
We are interested in the behaviour of the sequence as k goes to infinity, and our main result in this direction is the following.
that is, the asymptotic optimal domain is the square.
Note that we do not exclude the possibility of having non-uniqueness of for each k.
The above result may be interpreted in broad terms as follows. Given a positive number f, the rectangle which has more modes with frequency below f gets closer to the square as f goes to infinity.
3. Proof of theorem 2.1
The proof of the convergence to the square is divided into several parts. We begin by obtaining a lower bound for the eigenvalues of general rectangles, which improves upon Pólya's inequality in this particular case. This allows us in turn to show that optimal rectangles remain uniformly bounded with k, which is the key ingredient in the proof of theorem 2.1. We may then use the relations between the spectral problem and the lattice problem for convex domains for which much is known in terms of asymptotic behaviour in order to prove the desired convergence result.
(a) A lower bound for eigenvalues of Ra
For all rectangles Ra and all k we have 3.1
We note that the term is in fact always positive, and thus the above is always stronger than Pólya's lower bound in the particular case of rectangles. To see that this is the case, note that this is equivalent to λk>64a2π2/729 and, on the other hand, we have
Since the inequality in the above theorem holds for all rectangles, it will also hold for the sequence of optimal domains with a and λk replaced by and , respectively.
In a way similar to Pólya's proof for the case of rectangles , §§6 and 7, we write where The function thus defined is continuous and strictly increasing for y greater than or equal to 1. We may thus define its inverse, L, which will also be continuous and increasing.
Furthermore, given a circle centred at the origin with radius y1/2, the function L−1(y) yields the sum of the areas of rectangles belonging to the first quadrant which are inscribed on this circle and with one side on the x-axis at the segments (0,1),(1,2),…,(⌊y1/2⌋−1,⌊y1/2⌋).
We thus have that the area of the portion of circle in the first quadrant is larger than the sums of the areas of these rectangles, that is, 3.2
For the more general case of the ellipse with radii and given by we have that the number of integer lattice points N inside or on this ellipse (but not on the axes) is smaller than the area of the inscribed rectangles and thus may be bounded in terms of L−1 as follows: 3.3for y=λk/(a2π2). Hence the order k of the eigenvalue, which will, in general, be smaller than N owing to possible multiplicities, satisfies 3.4This is equivalent to equation (7.5) in Pólya's paper and, together with inequality (3.2), yields the well-known Pólya bound λk≥4kπ.
In order to improve this, we shall improve upon inequality (3.2). To do this, we add to the right-hand side the areas of the right-angled triangles for which one side is the top side of each of the above rectangles and another is the segment connecting the right top vertices of two consecutive rectangles—the remaining side is that which closes the triangle.
Furthermore, we shall also consider the right-angled triangle situated at the right-hand end of the circle whose basis is the segment (⌊y1/2⌋,y1/2) and the third vertex is located at the point .
We thus have that the sum of the areas of all the rectangles plus the sum of the areas of the triangles described earlier is still smaller than the area of the quarter circle and so as the sum on the right-hand side telescopes.
Let us now write y1/2=⌊y1/2⌋+y0, where y0∈[0,1). Then and we may write the above as 3.5We shall now prove that
This is equivalent to and since for y0 on [0,1], it follows that the term inside the square brackets on the left-hand side is always positive in this range. Since the term on the right-hand side is non-positive, the inequality follows.
We thus have which, again letting y=λk/(a2π2) and using (3.4), yields the desired result.
It is possible to obtain a stronger lower bound by using the maximum of the function of y0 on the right-hand side in (3.5) instead of the above. However, this will yield a more complicated expression while the above result is sufficient for our purposes.
(b) Proof of the boundedness of
As we saw in §1, the main obstacle to proving convergence to the square is the a priori possibility that the sequence will not be bounded. In this section, we prove that this is not the case.
There exists a value of a, say aM, such that for all positive integers k.
It is enough to prove that remains bounded as k goes to infinity. From theorem 3.1, we have Since for all k, the function is an increasing function of and thus the same is true for the left-hand side in the above inequality. Hence, this will also hold with replaced by the kth eigenvalue of the square, which we shall denote by , We thus have Since , we have , and thus This implies for all k. Since , we obtain
(c) Relation to lattice problems
Let us consider again the ellipse E with radii and . The expression for the number of integer lattice points inside E that are in the first quadrant (excluding the axes) was given in equation (3.3). If we consider instead the number of integer lattice points N0 inside the whole of E, we have the relation where r1 and r2 are as above and .
The problem of bounding the difference in absolute terms between N0 and the area of E, sometimes referred to as the lattice point discrepancy of the domain, is a classical one, and may be traced back to Gauss's estimate for the case of the circle . Since then, much progress has been made along these lines, and the best known estimates for the particular case of ellipses are given by Krätzel and Nowak [26,27]. From these papers, we have, for instance, that which, when translated into eigenvalue bounds, yields
The problem here is that the extra term on the right-hand side is not necessarily of one sign, and it thus becomes difficult to bound a uniformly.
Since we have already proved the boundedness of optimal rectangles, it will now be more convenient to use those results whose emphasis lies on the exponents controlling the asymptotic behaviour. In this respect, the best known result is that obtained by Huxley  for a general convex domain, which states that for some positive constant C. For our present purposes, it is enough to know that there exist constants θ smaller than 1 and C positive such that although we will get back to this issue later in the paper.
(d) Proof of theorem 2.1
We are now ready to prove the main theorem. Using the above estimates for the lattice problem we obtain, in terms of eigenvalues and for general rectangles, and 3.6where the Ci(i=1,2,3) are suitable positive constants.
If we apply the above to optimal eigenvalues, we obtain 3.7where now C* denotes the supremum of the constants C3 in (3.6) for each value of ; because this latter sequence is uniformly bounded, the supremum of such constants exists.
We again use the fact that the function x1/2−4kπx−1/2 is an increasing function of x to replace by in the above expression. Taking limits as k goes to infinity yields implying that converges to 1 as desired.
4. Determining optimal eigenvalues computationally
Let us now consider the problem of finding the rectangle minimizing the Dirichlet eigenvalue λk. For eigenvalues of very low order, it is possible to find the solution by hand, although the number of computations involved will be seen to grow quite quickly. As a consequence, the problem soon becomes highly non-trivial from a computational point of view, and, for large order, its solution does require the usage of significant computational resources. Moreover, the time needed for solving the problem will depend strongly on the choice of algorithm.
With this in mind, we shall now propose an algorithm to solve the large-order problem and then present some computational results.
Let be the rectangle defined in §2 with eigenvalues given by (2.1). The order of an eigenvalue is related to the position of the corresponding expression 4.1in the sorted list of terms obtained for all pairs of integers i and j.
For fixed integers m and p, we want to solve the optimization problem 4.2
(a) Computational algorithm for solving the optimization problem (eqn4.2)
The solution of problem (4.2) can be divided into two tasks. First, we use an efficient algorithm for the calculation of the vector of eigenvalues vm,p=(λm,λm+1,…,λm+p) of a given rectangle . We then determine the values of a that minimize each of the components of vm,p. In the next two sections, we describe the algorithms for solving these two tasks in detail.
(b) Calculation of the vector vm,p for a given rectangle
From the considerations made above, for a fixed rectangle , the determination of the corresponding mth eigenvalue is equivalent to finding the minimum value of cm such that the region contains (at least) m points in the set It is clear that for each m there exists (at least) one point on the boundary of the region ; otherwise, it would be possible to decrease the value of cm until some point hits the boundary. Moreover, the multiplicity of the eigenvalue λm coincides with the number of points in lying on the boundary of the region . It is possible to handle cm+p in a similar way and thus, to solve the optimization problem (4.2), it is sufficient to compute and then sort the list of eigenvalues corresponding to the integer lattice points contained in the region . The nth component of the vector vm,p=(λm,λm+1,…,λm+p) is the nth component of the sorted list. For each value of c, we define the regions and , which are (respectively) the largest and smallest polygonal regions with corners in such that and . We have and where ⌊x⌋ and ⌈x⌉ denote the greatest integer less than or equal to x and the smallest integer greater than or equal to x, respectively. In figure 1a, we show part of the boundaries of the regions and with c1<c2 and the corresponding polygonal regions and .
Let and Λ=λi,j be the νth eigenvalue.
(i) If , then ν<⌊πc2/4⌋.
(ii) If , then ν>⌈πc2/4⌉.
Observe first that from the equation i2/a2+j2a2=c we see that the order of the eigenvalue is non-decreasing with c. Since , the strict inclusion implies that . Moreover, is an integer and, because is the finite union of disjoint squares with unit area, it gives precisely the number of points . It then follows that ν<⌊πc2/4⌋, proving (i). The proof of (ii) is similar using the inclusion .
Lemma 4.1 suggests an algorithm for the solution of problem (4.2). Let c1 and c2 be two constants with c1<c2 such that but and but . The eigenvalues corresponding to points in do not need to be computed because their order is smaller than m, by lemma 4.1(i). Thus, we only need to know how many such points there are, which may be obtained by computing the area of . Now by lemma 4.1(ii) the components of the vector solution of problem (4.2) are eigenvalues corresponding to points in . Thus, the computational effort needed to solve problem (4.2) is reduced to the calculation of terms of the form (4.1) corresponding to points in this last region. We then sort those values to build vm,p.
Figure 1b shows the points considered by the algorithm. Those shown in black are points in that need not be computed, while those in grey are related to the points in .
Note that and thus and , which are good initial estimates for the determination of c1 and c2.
(c) Solving the optimization problem
Each component of the vector vm,p determined above for a particular rectangle may be seen as a function of a and we are interested in the corresponding minimization problem.
To give an idea of the complexity of the problem under consideration, we show λ100 000 as a function of a in figure 2a. From this, it is clear that the solution of the optimization problem is non-trivial because not only does the objective function have a highly oscillatory behaviour, but there are also many points where the function is not differentiable. Zooming in on a smaller interval (figure 2b) shows that we continue to have non-smooth variations, which make the optimization procedure quite difficult.
Our strategy to handle this problem can be divided into two steps, which will be described in the next sections.
(i) Initialization of the optimization procedure
We shall begin by computing the vectors vm,p for values of ai in the interval [1,Amax], obtained with increments of a small value δ. For each component of the vector, we define a set by picking M local minima where the function takes its lowest values. Thus, we run several minimization processes with the Nelder–Mead method (see ), starting with a sample of points randomly generated in the neighbourhood of each element in . Finally, we define two vectors and with length p whose ith component is equal to the minimum among all ith components of the vectors vm,p and to the corresponding value of a, respectively.
(ii) Refining the solutions
In this section, we describe how to refine the vectors , with the approximated optimal values of a and the corresponding minimal eigenvalues. For each eigenvalue, the optimal value a satisfies one of two things: either a is the minimizer of a function as in (4.1) for integers i and j, or there are (at least) two of these functions associated with two pairs of integers (i,j) and (p,q) that intersect at a. In the first case, we have and, if j≥i and a≥1, it follows that In this case, we can only have a minimum at the square (a=1). On the other hand, if j<i, Now, for each component in the vector , we have an approximated value of a, say , and in we have the corresponding approximated eigenvalue . Thus, there exists a pair of integers (i,j) such that It is straightforward to check whether for each admissible integer j there is such that is sufficiently close to . If that is the case, we compute the eigenvalue π2(i2/a2+j2a2), compare with the corresponding value in and, if necessary, update and .
The remaining optimal values of a are located at the intersection of two functions of the form (4.1), 4.3In a similar fashion, we can build a list with pairs of integers (i,j) such that the corresponding eigenvalue is sufficiently close to . Then, for each with (i,j)≠(p,q), we use (4.3) to calculate a candidate for the minimizer. Again, we determine the eigenvalue π2(i2/a2+j2a2), compare with the corresponding value in and, if smaller, update and .
5. Computational results
We shall now present the main computational results obtained with the algorithm described earlier. Figures 3 and 4 show the first 300 000 terms of the sequence and , respectively. Although the two figures together do indeed suggest that the sequence converges to 1, that is, optimal rectangles do approach the square as k goes to infinity, figure 3 by itself is not as conclusive and it is compatible with the optimal values of a oscillating between 1 and some value between 1.1 and 1.2.
This is a consequence of the fact that the decay as k goes to infinity is expected to be, in fact, quite slow. Indeed, if we expand the right-hand side in (3.7) as k goes to infinity—after having replaced the two occurrences of inside the square brackets by —we see that the rate of decay of to 2 is bounded from above by the term coming from the bound for the remaining term in the integer lattice point problem. We thus have that if the bound given by Hardy  is optimal, then this should be of order k−1/4+ε for all positive ε, while for we would get k−1/8+ε. Although this is an upper bound, we do see that such a decay is compatible with the values in figure 4.
Another aspect that is of interest to consider are the multiplicities of the optimizers. Regarding this, the optimal values can occur if the corresponding minimizes one function as in (4.1) for some integer numbers i and j, or if we have a multiplicity larger than 1. Let us define Mk as the multiplicity of the optimal eigenvalue at the point , that is, Mk=q if and only if Figure 5 shows the multiplicities Mk, k=1,2,…,100 000.
It is worth noting that, in most cases, the multiplicity of the optimal eigenvalue equals 2. However, there are two optimal eigenvalues (apart from and ) with multiplicity 1 in the range considered, namely M55 211=M87 366=1, and quite a large number of optimal eigenvalues do have larger multiplicities. For instance, from the computational results obtained, we see that and has multiplicity 12. Moreover, we have . Figure 6a,b shows λj, j=4265,4266,…,4276 as a function of a∈[1,1.001] and λj, j=5032,5033,…,5043 for a in a neighbourhood of , respectively.
The determination of the multiplicity of an optimal eigenvalue associated with an optimal is equivalent to determining the number of pairs such that In figure 7, we show a density plot of the pairs (i,j) involved in the calculation of the first 300 000 optimal eigenvalues. We can observe that the distribution of these pairs generates certain patterns with some regions having a high concentration of points while others are almost empty.
The results presented suggest the extension of this problem to other families of domains with the purpose of determining whether or not an asymptotic optimal domain exists and, if so, to characterize it and study the properties of the sequence of optimal eigenvalues. Of these problems, perhaps the most obvious (and interesting) to consider is that when domains are allowed to vary in the class of general bounded domains in with fixed volume. From the results in Bucur  and Mazzoleni & Pratelli , we know that a sequence of optimizers, say , exists within the context of quasi-open sets, and that each optimizer has a bounded perimeter. However, most of the ingredients used in our proof for rectangles are still missing in this case, such as the uniform control of the perimeter. We remark that, in this more general situation, it is not even clear that such control needs to exist, as even if the sequence of minimizers remains uniformly bounded the surface area may increase to infinity with k. This, or a similar effect, may cause the remaining term in (1.1) to become very large and negative for subsequent domains, causing the sequence of optimizers to not even converge.
An added complexity of this problem is the fact that now there will also exist issues about what type of convergence should be considered.
In spite of this, and although there are no numerical results to corroborate a conjecture either way at this point in time, the most natural guess here would be that does indeed approach the ball as k goes to infinity. We note that this is the case for planar domains when the restriction is imposed directly on the perimeter .
P.R.S.A. was partially supported by FCT, Portugal, through scholarship SFRH/BPD/47595/2008 and project PTDC/MAT/105475/2008 and by Fundação Calouste Gulbenkian through programme Estímulo à Investigação 2009. Both authors were partially supported by FCT's project PTDC/MAT/101007/2008.
- Received August 21, 2012.
- Accepted October 23, 2012.
- © 2012 The Author(s) Published by the Royal Society. All rights reserved.