On the behaviour of clamped plates under large compression

TL;DR

This paper analyzes how the eigenvalues of clamped plates behave under large compression, relating the problem to Laplacian eigenvalues with Robin boundary conditions, and studies extremal domain shapes numerically.

Contribution

It establishes the asymptotic behavior of eigenvalues under compression and explores the shape and nodal domain changes of extremal domains numerically.

Findings

Eigenvalues relate to Laplacian with Robin boundary conditions under large compression.

Extremal domain shapes depend on compression level and develop boundary structures.

Number of nodal domains increases with compression.

Abstract

We determine the asymptotic behaviour of eigenvalues of clamped plates under large compression, by relating this problem to eigenvalues of the Laplacian with Robin boundary conditions. Using the method of fundamental solutions, we then carry out a numerical study of the extremal domains for the first eigenvalue, from which we see that these depend on the value of the compression, and start developing a boundary structure as this parameter is increased. The corresponding number of nodal domains of the first eigenfunction of the extremal domain also increases with the compression.

Tables

Table 1. Table 1. Numerical approximations obtained for the first eigenvalue of the minimizer when α = 110 𝛼 110 \alpha=110 , for different values of m 𝑚 m .

$m$	${\tilde{\lambda }}_1$
1000	-1786.3537774
1500	-1786.3537779
1800	-1786.3537762
2000	-1786.3537753

Equations

$$\left\{\begin{array}[]{ll}\Delta^{2}u+\alpha\Delta u=\lambda u,&{\rm\ in\ }\Omega,\vspace*{2mm}\\ u=\frac{\displaystyle\partial u}{\displaystyle\partial\nu}=0,&{\rm\ on\ }\partial\Omega,\end{array}\right.$$

$${\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_k\le \dots ,$$

$$\left\{\begin{array}[]{ll}\Delta^{2}v+av+\gamma\Delta v=0,&{\rm\ in\ }\Omega,\vspace*{2mm}\\ v=\frac{\displaystyle\partial v}{\displaystyle\partial\nu}=0,&{\rm\ on\ }\partial\Omega,\end{array}\right.$$

$$\underset{a\to +\infty }{\lim }\frac{{\gamma }_1(a)}{\sqrt{a}}=2,$$

$$\underset{\alpha \to +\infty }{\lim }\frac{{\lambda }_1(\alpha )}{{\alpha }^2}=-\frac{1}{4}.$$

$${\lambda }_k(\Omega ,\alpha )=-\frac{{\alpha }^2}{4}+\mathrm{o}({\alpha }^2),$$

$${\lambda }_1(\Omega ,\alpha )=-\frac{{\alpha }^2}{4}+\mathrm{O}(\alpha ),$$

$$\lambda =-\frac{{\alpha }^2}{4}+\frac{c_1\alpha }{R^2}+\frac{c_2}{R^4}+\mathrm{o}(1),$$

$${\lambda }_1=-\frac{{\alpha }^2}{4}+\frac{{\pi }^2\alpha }{2R^2}+\frac{{\pi }^2(N^2-1-{\pi }^2)}{4R^4}+\mathrm{o}(1),$$

$${\int }_{\Omega }\Delta u\Delta \varphi -\alpha \nabla u\nabla \varphi =\lambda {\int }_{\Omega }u\varphi ,\forall \varphi \in H_0^2(\Omega ),$$

$${\lambda }_k^D(\Omega ,\alpha )=\underset{\begin{array}{c}V\subset H_0^2(\Omega )\\ \dim V=k\end{array}}{\min }\underset{0\ne u\in V}{\max }\frac{{\int }_{\Omega }(\Delta u{)}^2-\alpha |\nabla u{|}^2}{{\int }_{\Omega }{u}^2}.$$

$$(\Delta +{\alpha }_{+})(\Delta +{\alpha }_{-})u=0,$$

$${\alpha }_{+}=\frac{\alpha }{2}+\sqrt{\frac{{\alpha }^2}{4}+\lambda },{\alpha }_{-}=\frac{\alpha }{2}-\sqrt{\frac{{\alpha }^2}{4}+\lambda }.$$

$$u(r,\theta )=r^{1-\frac{N}{2}}\left[A{J}_{k+\frac{N}{2}-1}\left(r\sqrt{{\alpha }_{+}}\right)+B{I}_{k+\frac{N}{2}-1}\left(r\sqrt{-{\alpha }_{-}}\right)\right]S_k(\theta ),$$

$$\left\{\begin{array}[]{l}Af_{k}(R)+Bg_{k}(R)=0,\\ Af_{k}^{\prime}(R)+Bg_{k}^{\prime}(R)=0,\end{array}\right.$$

$$f_k(r)=r^{1-\frac{N}{2}}{J}_{k+\frac{N}{2}-1}\left(r\sqrt{{\alpha }_{+}}\right)\text{mbox}and{g}_k(r)=r^{1-\frac{N}{2}}{I}_{k+\frac{N}{2}-1}\left(r\sqrt{-{\alpha }_{-}}\right).$$

$$f_k(R)g_k^{\prime }(R)-g_k(R)f_k^{\prime }(R)=0,$$

$$\begin{array}[]{l}\frac{\displaystyle RJ_{k+\frac{N}{2}-1}\left(R\sqrt{\alpha_{+}}\right)}{\displaystyle kJ_{k+\frac{N}{2}-1}\left(R\sqrt{\alpha_{+}}\right)-R\sqrt{\alpha_{+}}J_{k+\frac{N}{2}}\left(R\sqrt{\alpha_{+}}\right)}\vspace*{2mm}\\ \hskip 85.35826pt=\frac{\displaystyle RI_{k+\frac{N}{2}-1}\left(R\sqrt{-\alpha_{-}}\right)}{\displaystyle kI_{k+\frac{N}{2}-1}\left(R\sqrt{-\alpha_{-}}\right)+R\sqrt{-\alpha_{-}}I_{k+\frac{N}{2}}\left(R\sqrt{-\alpha_{-}}\right)}.\end{array}$$

$$u(r,\theta )=r^{1-\frac{N}{2}}\left[A{J}_{k+\frac{N}{2}-1}\left(r\sqrt{{\alpha }_{+}}\right)+B{J}_{k+\frac{N}{2}-1}\left(r\sqrt{{\alpha }_{-}}\right)\right]S_k(\theta ),$$

$$\begin{array}[]{ll}\frac{\displaystyle J_{k+\frac{N}{2}-1}\left(R\sqrt{\alpha_{+}}\right)}{\displaystyle kRJ_{k+\frac{N}{2}-1}\left(R\sqrt{\alpha_{+}}\right)-\sqrt{\alpha_{+}}J_{k+\frac{N}{2}}\left(R\sqrt{\alpha_{+}}\right)}\vspace*{2mm}\\ \hskip 85.35826pt=\frac{\displaystyle J_{k+\frac{N}{2}-1}\left(R\sqrt{\alpha_{-}}\right)}{\displaystyle kRJ_{k+\frac{N}{2}-1}\left(R\sqrt{\alpha_{-}}\right)-\sqrt{\alpha_{-}}J_{k+\frac{N}{2}}\left(R\sqrt{\alpha_{-}}\right)}.\end{array}$$

$$\left\{\begin{array}[]{ll}\Delta^{2}u=-\Lambda\Delta u,&{\rm\ in\ }\Omega,\\ u=\frac{\partial u}{\partial\nu}=0,&{\rm\ on\ }\partial\Omega,\end{array}\right.$$

$$u(r,\theta )=\left[A{r}^{1-\frac{N}{2}}{J}_{k+\frac{N}{2}-1}\left(r\sqrt{\alpha }\right)+B{r}^k\right]S_k(\theta ).$$

$$J_{k+\frac{N}{2}}\left(R\sqrt{\alpha }\right)=0.$$

$$\left\{\begin{array}[]{ll}-\Delta u=\sigma u,&{\rm\ in\ }\Omega,\\ \frac{\partial u}{\partial\nu}+\beta u=0,&{\rm\ on\ }\partial\Omega.\end{array}\right.$$

$$\left\{\begin{array}[]{ll}-\Delta u=\gamma u,&{\rm\ in\ }\Omega,\\ u=0,&{\rm\ on\ }\partial\Omega,\end{array}\right.$$

$$0<{\gamma }_1\le {\gamma }_2\le \cdots \to +\infty .$$

$$r^{1-\frac{N}{2}}{J}_{k+\frac{N}{2}-1}(r\sqrt{\sigma })S_k(\theta ),$$

$$\left(\frac{k}{R}+\beta \right)J_{k+\frac{N}{2}-1}(R\sqrt{\sigma })=\sqrt{\sigma }{J}_{k+\frac{N}{2}}(R\sqrt{\sigma });$$

$$r^{1-\frac{N}{2}}{I}_{k+\frac{N}{2}-1}(r\sqrt{-\sigma })S_k(\theta ),$$

$$\left(\frac{k}{R}+\beta \right)I_{k+\frac{N}{2}-1}(R\sqrt{-\sigma })=-\sqrt{-\sigma }{I}_{k+\frac{N}{2}}(R\sqrt{-\sigma });$$

Full text

1. Introduction

Let $\Omega$ be a smooth bounded domain in $\mathbb{R}^{N}$, $N\geq 2$. We are interested in the following eigenvalue problem

[TABLE]

considered as a model for a clamped plate. Here $\alpha$ is a real parameter corresponding to the quotient between the tension and the flexural rigidity and, depending on its sign, represents whether the plate is under tension ($\alpha<0$) or compression ($\alpha>0$). For domains $\Omega$ as described above, the eigenvalues of (1.1) form an infinite sequence

[TABLE]

where $\lambda_{k}=\lambda_{k}(\Omega,\alpha)$ approaches $+\infty$ as $k$ goes to infinity.

The study of this and similar problems has been considered in the literature continuously over time since the works of Lord Rayleigh [18] and Love [21] on clamped plates. We refer the reader to the book [15] for an extensive historical and scientific overview on the mechanics of plates through the Kirchhoff-Love model, which leads to problem (1.1).

In this paper, we are concerned with two issues related to (1.1), namely, the asymptotic behaviour of the eigenvalues $\lambda_{k}$ as the parameter $\alpha$ approaches $+\infty$ (the case of $-\infty$ was considered in [14]), and the extremal domains of such eigenvalues as $\alpha$ varies. In the first instance, the above problem is closely related to

[TABLE]

where now the eigenvalue parameter is $\gamma=\gamma(a)$, and the (positive) parameter $a$ stands for the elasticity constant of the medium surrounding the plate. We know from a result in [19] that for (1.2)

[TABLE]

which when translated into the eigenvalue problem in (1.1) yields

[TABLE]

Our main result along these lines is to extend this to all eigenvalues $\lambda_{k}$. This is achieved by a different approach from that used in [19], involving now a connection which, to the best of our knowledge, is new, between the eigenvalues of the clamped plate problem (1.1) and those of a Robin eigenvalue problem for the Dirichlet Laplacian in the case where $\Omega$ is a ball – see Section 3 for the details. To be more precise, we prove the following

Theorem 1.1** (Asymptotic behaviour of the $k^{\rm th}$ eigenvalue).**

Let $\Omega$ be a bounded domain in $\mathbb{R}^{N}$. Then, for any positive integer $k$, the eigenvalues of (1.1) satisfy

[TABLE]

as $\alpha\to+\infty$. Moreover,

[TABLE]

as $\alpha\to+\infty$.

For positive values of $\alpha$, each of the eigenvalue curves $\lambda_{k}=\lambda_{k}(\alpha)$ is, in fact, made up of analytic eigenvalue branches which intersect each other - see Figure 1, where to illustrate this effect we plotted the quantity $\lambda_{k}(\Omega,\alpha)+\frac{\displaystyle\alpha^{2}}{\displaystyle 4}$ for the disk and for ellipses.

This branch-switching phenomenon makes it much more difficult to obtain further terms in the asymptotic expansion and it is the independence of the first term on the order of the eigenvalue which allows us to derive the expansion for all $k$. In the particular case where $\Omega$ is a ball of radius $R$, which is at the heart of the proof of Theorem 1.1, we are able to prove that the number of such eigenvalue branches which make up the $k^{\rm th}$ eigencurve is finite, and we determine further terms in the asymptotic expansion of these analytic branches. These results are summarised in the following

Theorem 1.2** (Asymptotic behaviour of analytic eigenvalue branches for balls).**

For any analytical branch of the eigenvalues of problem (1.1) when $\Omega$ is a ball of radius $R$, we have

[TABLE]

as $\alpha\to+\infty$, where $c_{1}$ and $c_{2}$ are constants depending on the eigenvalue branch, with $c_{1}$ being positive. In the case of the first eigenvalue we have

[TABLE]

The full description of the coefficients $c_{1}$ and $c_{2}$ may be found in Theorem 3.3 in Section 3.

It is possible to consider problem (1.1) with other boundary conditions, such as the Navier setting. This is not as interesting from a mathematical perspective, since the problem then reduces directly to the study of the second order elliptic operator $\Delta+\alpha/2$. However, and as we show in Section 4, there is a major difference between the Dirichlet and Navier cases in that for the Navier problem the number of crossings of analytic branches to make up an eigencurve corresponding to the $k^{\rm th}$ eigenvalue is actually infinite for each $k$. Complex crossing and avoided-crossing patterns seem to be a characteristic of such systems in the large compression regime, and they have also been identified in the one-dimensional fourth-order problem with different boundary conditions studied in [11].

Concerning our second topic of study, namely extremal domains for eigenvalues of problem (1.1), even in the case where the parameter $\alpha$ vanishes the problem is known to be extremely difficult with results available only in two and three dimensions – see [22, 5], respectively; for (1.2) there are no complete results in any dimension. Once $\alpha$ is taken to be nonzero in (1.1), the only existing result is an extension to sufficiently small positive values of $\alpha$ in two dimensions [4]. Our purpose in this part is thus mainly to provide a numerical exploration of the different types of extremal domains under an area restriction, showing in particular that the ball is no longer a minimizer for large compression.

We consider the numerical solution of the eigenvalue problem (1.1) using the Method of Fundamental Solutions (see e.g., [1, 3]). This is a meshless numerical method where the approximation is made by a discretization of an expansion in terms of the single and double layer potentials. In particular, by construction, the numerical approximation satisfies the fourth-order partial differential equation and we can focus just on the approximation of the boundary conditions of the problem. The computational implementation of this numerical method is described in Section 6.1 and some numerical results for the shape optimization problem are presented in Section 6.3. In particular, we will study minimizers of the first eigenvalue of problem (1.1) subject to an area constraint. The obtained numerical results suggest that the minimizer depends on the parameter $\alpha$, with the ball being the minimizer for all negative $\alpha$ and then extending to $\alpha\in[0,\alpha^{\star}]$, for some positive $\alpha^{\star}$. Note that this last result corresponds to that proved in [4] for sufficiently small $\alpha$, with our numerical simulations suggesting that, in fact, one may take $\alpha^{\star}$ to be at least as large as the first buckling eigenvalue. For large values of the parameter $\alpha$ we obtain some non-trivial minimizers - see Figure 3.

This numerical study has been performed mainly among general simply connected domains. However, we performed also the optimization of the first eigenvalue of problem (1.1) among annuli having unit area and compared the optimal values that were obtained with the corresponding values of the ball. These results suggest that the first eigenvalue of the disk is always smaller than the corresponding eigenvalue of the optimal annulus, independently of the parameter $\alpha$.

2. Statement of the problem

We start by observing that problem (1.1) has the following weak formulation

[TABLE]

and its eigenvalues may be described through their variational characterizations

[TABLE]

In what follows, whenever the meaning is clear from the context, we will drop either argument in $\lambda_{k}^{D}(\Omega,\alpha)$ for the sake of simplicity.

In order to determine the eigenfunctions of problem (1.1) when $\Omega=B_{R}(0)$, we rewrite equation (1.1) as

[TABLE]

where

[TABLE]

Both $\alpha_{+}$ and $\alpha_{-}$ are always real, as may be seen from inequality (4.6), and $\alpha_{+}$ is always positive while the sign of $\alpha_{-}$ depends on the sign of the eigenvalue $\lambda$.

For positive $\lambda$ it is known that the solution of (2.3) can be written as (cf. [5])

[TABLE]

where $J_{\nu}$ and $I_{\nu}$ are the Bessel and the modified Bessel functions, respectively, of the first kind of order $\nu$, and $S_{k}$ are the spherical harmonic functions of order $k$. The boundary conditions then yield the following system of equations

[TABLE]

where we have set

[TABLE]

Since we are interested in the existence of nontrivial solutions of system (2.6), we impose the corresponding determinant to be zero, namely

[TABLE]

from which we obtain the corresponding eigenvalues and, as a consequence, the general form of the eigenfunctions. Using standard Bessel function identities, equation (2.7) may be rewritten as

[TABLE]

When $\lambda$ is strictly negative, then $\alpha_{-}$ is strictly positive and, in place of (2.5), we now have

[TABLE]

where the coefficients are given by a system similar to (2.6), and the eigenvalues are now solutions of

[TABLE]

Finally, it remains to consider the case $\lambda=0$, which behaves in a slightly different way. We note that then $\alpha_{+}=\alpha$ while $\alpha_{-}=0$, and in particular this means that $\alpha$ has to be an eigenvalue of the following buckling problem

[TABLE]

for which the eigenfunctions are known to be of the form (they can be derived in a similar way as for the other cases)

[TABLE]

In particular, $\alpha$ has to be a solution of the following

[TABLE]

Moreover, if $\alpha$ is the $k$-th eigenvalue $\Lambda_{k}$ of the buckling problem (2.11), we immediately deduce that the vanishing eigenvalue of problem (1.1) is exactly the $k$-th one $\lambda_{k}$, and the multiplicity will be the same of $\Lambda_{k}$.

3. Connection to the Robin Laplacian

Even though there are no simple relations between the Laplacian and the Bilaplacian in general (apart for the Navier problem (4.1)), if we consider the generic situation of problem (1.1) in the ball in $\mathbb{R}^{N}$ with $\alpha\in\mathbb{R}$, we can draw a very precise connection to the Robin Laplacian.

To this end we recall that the Robin problem for the Laplace operator is as follows

[TABLE]

For any real value of $\beta$ the corresponding spectrum consists of a non-decreasing sequence of eigenvalues with finite multiplicities diverging to plus infinity. In particular, for positive values of $\beta$ the eigenvalues are all strictly positive, while for $\beta=0$ the Robin problem (3.1) becomes the Neumann problem. It is also known that, as $\beta\to+\infty$, problem (3.1) converges to the Dirichlet problem for the Laplace operator, namely

[TABLE]

whose eigenvalues we will denote by

[TABLE]

We recall that the eigenfunctions of (3.1) on a ball can be sorted out into three categories (cf. Section 2):

• i.

if the eigenvalue $\sigma$ is positive, the eigenfunction is of the form

[TABLE]

and eigenvalues are solutions of

[TABLE]

• ii.

if the eigenvalue $\sigma$ is negative, the eigenfunction is of the form

[TABLE]

and eigenvalues are solutions of

[TABLE]

• iii.

if the eigenvalue $\sigma$ is zero, the eigenfunction is of the form $r^{k}S_{k}(\theta)$, and in particular this occurs when $\beta=-\frac{k}{R}$.

At this point it is clear that any eigenfunction of the clamped plate problem (1.1) on the ball can be thought of as the sum of two different eigenfunctions of the Robin problem for the Laplacian (3.1). A first condition that these two Robin eigenfunctions have to satisfy is that their spherical parts coincide. This implies that they must come from two different eigenvalues and, in particular, we have that these two eigenvalues are $\alpha_{+}$ and $\alpha_{-}$, and the multiplicities must coincide. Furthermore, such eigenfunctions must have the same index $\nu=k+\frac{N}{2}-1$ in their Bessel function part. Let us call $v_{1},v_{2}$ two such eigenfunctions (with associated eigenvalues $\sigma_{1},\sigma_{2}$) and let

[TABLE]

i.e., we denote by $v_{j}^{R}$ the radial part. Since we want the boundary conditions in the clamped plate problem (1.1) to be satisfied, the only way to combine $v_{1}$ and $v_{2}$ is to set

[TABLE]

as can be easily checked from the boundary conditions in the Robin problem for the Laplace operator (3.1). In particular, $v_{1}$ and $v_{2}$ must be Robin eigenfunctions associated with the same parameter $\beta$. As for the equation, we observe that

[TABLE]

and the equality $\sigma_{1}^{2}-\alpha\sigma_{1}=\sigma_{2}^{2}-\alpha\sigma_{2}$ is naturally satisfied since

[TABLE]

On the other hand, letting $\beta$ go to infinity we get that, for specific values of $\alpha$ and $\lambda$, we should consider eigenvalues of the Dirichlet Laplacian (3.2) instead. From the literature (see e.g., [8] and the references therein, and also [25] for a study of the first two Robin eigenvalues) we know that all the analytical branches related to the same Bessel function $J_{k+\frac{N}{2}-1}$ ($I_{k+\frac{N}{2}-1}$ when the eigenvalue is negative, $r^{k}$ if zero) can be continued at $\beta=\infty$ generating a function which wraps around $\mathbb{R}$ infinitely many times. If we call $\sigma_{k,j}(\beta)$ the $j$-th eigenvalue associated with $J_{k+\frac{N}{2}-1}$, then the analytical branches of eigenvalues of problem (1.1) are given by

[TABLE]

for some $t\in\mathbb{N}$, where the parameter $j$ is of no relevance here since any time $\beta$ reaches infinity $j$ has to be replaced by $j+1$ as

[TABLE]

where $\gamma_{k,j}$ is the $j$-th eigenvalue of the Dirichlet problem for the Laplace operator (3.2) associated with $J_{k+\frac{N}{2}-1}$. In particular, different branches of eigenvalues of the clamped plate problem (1.1) associated with the Bessel index $k+\frac{N}{2}-1$ are indexed by the parameter $t$ in (3.5). We remark that all the branches are of this type, hence no other branches are present. We sum up all these arguments in the following

Theorem 3.1**.**

Let $\beta\in\mathbb{R}$ and $v_{1}$ and $v_{2}$ be any two eigenfunctions of problem (3.1) in a ball $B_{R}(0)$ associated with the eigenvalues $\sigma_{1}$ and $\sigma_{2}$, respectively, and having the same spherical part, namely

[TABLE]

Then the function defined in (3.3) is an eigenfunction of problem (1.1) in the ball $B_{R}(0)$ associated with the eigenvalue $\lambda=-\sigma_{1}\sigma_{2}$ and with the parameter $\alpha=\sigma_{1}+\sigma_{2}$.

This representation completely characterizes the analytic branch of the eigenvalue $\lambda=-\sigma_{1}\sigma_{2}$ (for $\alpha\in\mathbb{R}$) as the parameter $\beta$ varies. In the limits $\beta\to\pm\infty$ we have that the eigenvalue $\lambda$ can be written as a product of eigenvalues of the Dirichlet problem for the Laplace operator (3.2), and the corresponding eigenfunction can also be written as a combination of eigenfunctions of problem (3.2).

In addition, all analytic branches of eigenvalues of the clamped plate problem (1.1) can be represented in this fashion.

Theorem 3.1 allows us to study the behaviour of the eigenvalues as $\alpha=\pm\infty$. Actually, for the case $\alpha=-\infty$, the convergence is well known in the literature for any smooth domain (see e.g., [14, p. 392]).

Theorem 3.2**.**

Let $w_{\alpha}$ be the eigenfunction associated with $\lambda_{k}(\alpha)$, and suppose that there exists a point $\alpha_{0}\in\mathbb{R}$ such that $w_{\alpha}\in C^{5}(\Omega)$ for any $\alpha<\alpha_{0}$. Then

[TABLE]

as $\alpha\to-\infty$, where $u_{k}$ is an eigenfunction of the Dirichlet problem for the Laplace operator (3.2) associated with $\gamma_{k}$.

We recall that, thanks to classical regularity theory for elliptic operators (cf.[15]), if $\Omega\in C^{5,\delta}$ then $w_{\alpha}\in C^{5,\delta}(\Omega)$ for any $\alpha\in\mathbb{R}$ and, in particular, balls satisfy the hypotheses of Theorem 3.2. It is easily seen then that we can recover the first term of the asymptotics (3.6) using the known asymptotics for the Robin problem (see [26] and the references therein).

Regarding the asymptotics as $\alpha\to+\infty$, we compute it using the knowledge that for any given branch when we get to $\beta=\infty$ we obtain that both $\alpha$ and $\lambda$ can be expressed in terms of zeros of Bessel functions:

[TABLE]

where $j_{\nu,m}$ is the $m$-th zero of $J_{\nu}$, whose asymptotic behaviour is known to be (cf. [23, formula (10.21.19)])

[TABLE]

as $m\to\infty$.

Let us now denote by $\psi_{m}$ and $\psi_{m+t}$ two eigenfunctions of the Dirichlet problem for the Laplace operator (3.2) associated with $\gamma_{m}=R^{-2}j_{k+\frac{N}{2}-1,m}^{2}$ and $\gamma_{m+t}=R^{-2}j_{k+\frac{N}{2}-1,m+t}^{2}$, respectively, having the same spherical part, and normalized such that $\psi_{m}+\psi_{m+t}$ is an eigenfunction of the clamped plate problem (1.1) under condition (3.7). Then using the Rayleigh quotient representation of $\lambda$ we have

[TABLE]

We recall that in this particular case we have $\alpha=(j_{\nu,m}^{2}+j_{\nu,m+t}^{2})R^{-2}$, where we set $\nu=k+\frac{N}{2}-1$ for simplicity. We now compute the asymptotics for the remainder in (3.9) and get

[TABLE]

telling us that $\lambda\sim-\frac{\displaystyle\alpha^{2}}{\displaystyle 4}+\frac{\displaystyle\alpha t^{2}\pi^{2}}{\displaystyle 2R^{2}}+{\rm o}(\alpha)$. Going further we can get

[TABLE]

and hence we have

Theorem 3.3**.**

For any analytical branch of the eigenvalues of problem (1.1) on a ball $B_{R}$ of radius $R$, we have

[TABLE]

as $\alpha\to+\infty$, where $\nu=k+\frac{N}{2}-1$ is the index of the associated Bessel functions, and $t$ is the parameter introduced in (3.5).

We observe that, even if at a first glance the presence of the parameter $t$ may seem unnatural, it may be compared for example with the ordering number for zeros of Bessel functions $j_{\nu,k}$. From this perspective, it is natural that it appears in formula (3.10).

We are now ready to prove Theorem 1.1.

Proof of Theorem 1.1.

In the case of a general domain $\Omega$, we shall denote the radius of the largest inscribed ball and that of the smallest ball containing $\Omega$ by $R_{i},R_{c}$ respectively. By the inclusion properties for problem (1.1), we know that any eigenvalue of $\Omega$ is bounded from above and from below by the corresponding eigenvalues of the inscribed and circumscribed balls, respectively. This immediately proves (1.3). For higher eigenvalues it will, in general, be difficult to determine the precise order of each eigenvalue, but in the case of the first eigenvalue it is possible to identify the corresponding branch, namely that obtained by making $t=1$ and $k=0$, and in turn obtain the following (asymptotic) expression

[TABLE]

which implies (1.4). ∎

4. The Navier problem

We now turn our attention to the following eigenvalue problem

[TABLE]

for any $\alpha\in\mathbb{R}$. We immediately notice the resemblance of the Navier problem (4.1) with problem (1.1), as its weak formulation reads

[TABLE]

the only difference between this and (2.1) being the ambient space. In particular, comparing the variational characterization (2.2) of the eigenvalues of problem (1.1) with that of the eigenvalues of the Navier problem (4.1)

[TABLE]

yields

[TABLE]

Now we want to compute eigenfunctions and eigenvalues of the Navier problem (4.1). We can of course proceed as for the Dirichlet case in Section 2. However, we observe that we can modify the problem as follows

[TABLE]

which tells us immediately that, if the domain has the cone property, the Navier operator in (4.4) is the square of the translated Dirichlet Laplace operator $\Delta+\frac{\alpha}{2}$ (cf. [15]). In particular, if we denote by $\gamma_{k}$ the $k$-th eigenvalue of the Dirichlet Laplacian (3.2), we get that the spectrum of (4.1) is given by

[TABLE]

for any $\alpha\in\mathbb{R}$ and for any (smooth enough) domain $\Omega$. We remark that, for $\alpha<0$ (actually, for $\alpha<2\gamma_{1}$) we have

[TABLE]

for any $k$, while on the other hand we actually have intersections of the branches (the intersection points will depend on $\Omega$). However, we can still say that

[TABLE]

that is

[TABLE]

Theorem 4.1**.**

Let $\Omega$ be a bounded open set in $\mathbb{R}^{N}$ with the cone property. Then, for any $k\in\mathbb{N}$

[TABLE]

as $\alpha\to+\infty$. Moreover,

[TABLE]

as $\alpha\to+\infty$.

Proof.

Equality (4.7) easily follows from the inequality chain (4.6) coupled with the asymptotic expansion (1.3).

As for (4.8), we first observe that

[TABLE]

and for the choice $\alpha=2\gamma_{k}$ we have

[TABLE]

This alone is not enough to prove the asymptotic behaviour. However, we know that $\lambda_{1}^{N}(\alpha)$ is a polygonal line and that each and every segment is tangent to the asymptotic curve (thanks to (4.9)). It is thus enough to show that the vertices have the same asymptotic behaviour, i.e., the points $\alpha=\gamma_{k}+\gamma_{k+1}$ for which

[TABLE]

or equivalently

[TABLE]

therefore we have to show that

[TABLE]

To this end we recall the Weyl asymptotics for the Dirichlet eigenvalue problem for the Laplacian (3.2), namely,

[TABLE]

where $C_{1}$ and $C_{2}$ are (known) constants depending only on $\Omega$ and the dimension $N$. From the binomial Taylor expansion

[TABLE]

we have

[TABLE]

which clearly goes to zero for $N$ larger than one. ∎

Remark 4.2*.*

If the domain is not bounded, it is still possible to prove (4.8) without using the asymptotics (1.3) while following the same strategy we used in the previous proof. In particular, in order to get the term $-\frac{\alpha^{2}}{4}$, it suffices to show that

[TABLE]

which follows from the equality

[TABLE]

and the fact that the ratio of consecutive eigenvalues converges to $1$, thanks to Weyl’s asymptotics (4.11).

Also, it is clear from (4.12) that the term ${\rm o}(\alpha)$ in (4.8) is sharp, since a different exponent in the denominator in the limit (4.10) would not go to zero as $k\to\infty$.

Remark 4.3*.*

We observe that the behaviour of the clamped plate problem (1.1) and that of the Navier problem (4.1) are substantially different. On the one hand, from the asymptotics (3.10) we have that the branches of eigenvalues of the clamped plate problem (1.1) will stop intersecting for some sufficiently large value of $\alpha$, at least in the case of balls where the parameters $t$ and $k$ provide a clear ordering of the branches, so that it is in principle possible to see which branch will eventually be the $k$-th eigenvalue. On the other hand, we know a priori that the branches of eigenvalues of the Navier problem (4.1) will have an infinite number of intersections, making it quite complicated to decide which is the $k$-th eigenvalue. In particular, the knowledge of the behaviour of each individual branch does not provide sufficient information on the asymptotics of the eigenvalues. Similarly, even though the eigenspaces do not depend on $\alpha$, that associated with the $k^{\rm th}$ eigenvalue will keep on changing, creating a strange phenomenon of non-convergence.

5. Shape derivatives

We will now consider the problem of finding extremal domains for the $k$-th eigenvalue of problem (1.1), namely,

Problem 1**.**

Determine

[TABLE]

We observe that proving existence for Problem 1 within a specific class of domains can be quite difficult and, to the best of our knowledge, there are no results available in general. To gauge the difficulties involved, we refer the reader to [7] for a survey on existence results for the Laplacian case, for which it is still not known if existence holds within the class of open sets.

We will focus now on Problem 1 with $k=1$. We begin by deriving the formula for the Hadarmard shape derivative of an eigenvalue of (1.1). Note that the formula in the case $\alpha=0$ was already derived in a general setting and for multiple eigenvalues, see [10, 24]. We also refer to [9] and the references therein for a complete discussion on Hadamard formulas for the Biharmonic operator, also in the case $\alpha\neq 0$. Nevertheless, for the sake of simplicity we show here how to derive it in our specific case.

Consider an application $\Psi(t)$ such that $\Psi:t\in[0,T[\rightarrow W^{1,\infty}(\mathbb{R}^{N},\mathbb{R}^{N})$ is differentiable at 0 with $\Psi(0)=I,\ \Psi^{\prime}(0)=V$, where $W^{1,\infty}(\mathbb{R}^{N},\mathbb{R}^{N})$ is the set of bounded Lipschitz maps from $\mathbb{R}^{N}$ into itself, $I$ is the identity and $V$ is a given deformation field.

We will use the notation $\Omega_{t}=\Psi(t)(\Omega)$, for a given set $\Omega$, $\lambda_{n}(t):=\lambda_{n}(\Omega_{t},\alpha)$, $u_{t}$ is an associated eigenfunction with unitary $L^{2}$ norm, and $u^{\prime}$ will denote the derivative of $u_{t}$ at $t=0$. Moreover, we assume that $\lambda_{n}(0)$ is simple.

It is well known (see e.g., [13]) that if we define

[TABLE]

for some function $y$, then the Hadamard shape derivative is given by

[TABLE]

As a consequence we have

Theorem 5.1**.**

Let $\Omega$ be a bounded open set of class $C^{4}$. The Hadamard shape derivative for a simple eigenvalue $\lambda$ of problem (1.1) with corresponding eigenfunction $u$ is given by

[TABLE]

Proof.

We have

[TABLE]

and the eigenfunction is normalized,

[TABLE]

The function $u^{\prime}$ can be calculated by solving the following boundary value problem (c.f. [16, 17])

[TABLE]

Since the case $\alpha=0$ can be recovered from [24] (and can be done similarly to what follows), we assume $\alpha\neq 0$ and the eigenvalue equation can be written as

[TABLE]

so that we have

[TABLE]

Applying now formula (5.1) to the equation (5.3) and using (5.6) we obtain

[TABLE]

The proof is concluded once we observe that $u\in H^{4}(\Omega)$ (cf. [15]), and since $u=\frac{\partial u}{\partial\nu}=0$ on $\partial\Omega$, we have that

[TABLE]

∎

Remark 5.2*.*

Using formula (5.2), we may try to attack Problem 1 via the Lagrange Multiplier Theorem. Since the constraint here is $|\Omega|=1$, we obtain the following condition

[TABLE]

Note that condition (5.7) has then to be added to problem (1.1), yielding an overdetermined problem resembling the Serrin problem (see [27]). However, problem (1.1) coupled with condition (5.7) is a more difficult problem, and the only partial result available in the literature can be found in [12].

It is worth observing that solving the overdetermined problem (1.1), (5.7) is not equivalent to solving Problem 1: in fact, the former provides just a critical point, that may be only a local minimizer, or even a local maximizer. Interestingly enough, though, eigenfunctions on the ball always satisfy condition (5.7). For a more detailed analysis of this fact, we refer to [9, 10].

6. Numerical Methods

6.1. Numerical solution of the eigenvalue problem

In this section we will describe a numerical method for solving (1.1).

A fundamental solution $\Phi_{\lambda}$ of the partial differential equation of the eigenvalue problem (1.1) is given by (see e.g., [20])

[TABLE]

where $H_{0}^{(1)}$ is a Hankel function of the first kind.

We will consider particular solutions of the partial differential equation of the eigenvalue problem (1.1), by defining the boundary integral operators $(x\in\Omega)$

[TABLE]

where $\hat{\Gamma}$ is an artificial boundary that surrounds $\partial\Omega$ (see e.g., [1, 3]), and $\varphi$ and $\psi$ are densities. The numerical approximation of an arbitrary solution of the PDE of the eigenvalue problem (1.1) can be justified by density results e.g., [1, 2]. Moreover, we will assume that $\hat{\Gamma}$ does not intersect $\bar{\Omega}$. Thus, we can discretise the boundary integral operators by considering the linear combinations

[TABLE]

where $y_{m,j}$ are some points on $\hat{\Gamma}.$ Note that the functions $u_{m}$ are particular solutions of the partial differential equation involved in the eigenvalue problem (1.1) and the coefficients can be determined by fitting the boundary conditions of the problem.

We consider some collocation points $x_{1},...,x_{m}$, (almost) uniformly distributed on $\partial\Omega$ and impose the boundary conditions of the problem which leads to the $(2m)\times(2m)$ system

[TABLE]

We will consider the choice for source points $y_{m,j}$ described in [1], assume that $\nu_{y_{m,j}}=\nu_{x_{j}}$, and denote this vector simply by $\nu_{j}$. Using the notation $d_{i,j}=x_{i}-y_{m,j}$, the system (6.3) can be rewritten as

[TABLE]

The approximations of the eigenvalues can be calculated by adapting the Betcke-Trefethen method (see [6]) to this context. We consider $p$ points $z_{1},z_{2},...,z_{p}$, randomly chosen in $\Omega$ and define the following six blocks

\begin{array}[c]{lll}A(\lambda)=\left[\Phi_{\lambda}(d_{i,j})\right]_{m\times m},&&B(\lambda)=\left[\nu_{j}\cdot\nabla\Phi_{\lambda}(d_{i,j})\right]_{m\times m},\\ &&\\ C(\lambda)=\left[\nu_{i}\cdot\nabla\Phi_{\lambda}(d_{i,j})\right]_{m\times m},&&D(\lambda)=\left[\nu_{i}\cdot\nabla\left(\nu_{j}\cdot\nabla\Phi_{\lambda}(d_{i,j})\right)\right]_{m\times m},\\ &&\\ E(\lambda)=\left[\Phi_{\lambda}(\tilde{d}_{i,j})\right]_{p\times m},&&F(\lambda)=\left[\nu_{j}\cdot\nabla\Phi_{\lambda}(\tilde{d}_{i,j})\right]_{p\times m},\end{array}

where $\tilde{d}_{i,j}=z_{i}-y_{m,j}$. Then, we define the matrix

[TABLE]

compute the $\mathbf{QR}$ decomposition of $\mathbf{M}(\lambda)$, and calculate the minimal eigenvalue of the first $(2m)\times(2m)$ block of the matrix $\mathbf{M}(\lambda)$ that will be denoted by $\sigma_{1}(\lambda)$. The approximations for the eigenvalues of problem (1.1) are the values $\lambda$, for which $\sigma_{1}(\lambda)\approx 0$.

6.2. Numerical shape optimization

In this section we will consider the shape optimization Problem 1 among general simply connected planar domains, whose boundary can be parametrized by

[TABLE]

for some continuous and $(2\pi)$-periodic functions $\Gamma_{1}$ and $\Gamma_{2}$. We will consider the (truncated) Fourier expansions

[TABLE]

and

[TABLE]

for a sufficiently large $P\in\mathbb{N}$, and the optimization procedure consists in finding optimal coefficients $a_{j}^{(1)}$, $b_{j}^{(1)}$, $a_{j}^{(2)}$, $b_{j}^{(2)}$. The optimization is performed by a gradient-type method, using the Hadamard shape derivative given by Theorem 5.1 to calculate the derivative of the eigenvalue with respect to perturbations of the coefficients $a_{j}^{(1)}$, $b_{j}^{(1)}$, $a_{j}^{(2)}$, $b_{j}^{(2)}$.

6.3. Numerical results

In this section we present the main results that we gathered with our numerical procedure for solving Problem 1.

As was mentioned in the Introduction, each of the eigenvalue curves $\lambda_{k}(\alpha)$ is made up of analytic eigenvalue branches which intersect each other. We illustrate this fact in Figure 1. As was shown in Theorem 1.1, all the eigenvalues have the following asymptotic behaviour

[TABLE]

Thus, in order to produce more convenient pictures, instead of plotting the first eigenvalues as functions of $\alpha$, we will extract the first term of the expansion, which is the same for all eigenvalues, i.e., in Figure 1 we plot the quantities

[TABLE]

for a disk of unit area and similar results for an ellipse with unit area and eccentricity equal to $\sqrt{3}/2$.

Figure 2 shows the curve of the quantity $\lambda_{1}^{\ast}(\alpha)+\frac{\alpha^{2}}{4}$. We can observe several branches corresponding to different types of minimizers. Some of them, obtained for $\alpha=110,170,230,400$, are plotted in Figure 3. The optimal eigenvalue $\lambda_{1}^{\ast}(\alpha)$ is the minimum among the values obtained for all the branches. We calculated the critical value of $\alpha$, which is the maximal value of $\alpha$ for which the ball is the minimizer and obtained $\alpha^{\star}\approx 102.23$. In [4] it was proved that the ball is the minimizer for $\alpha\in[0,a],$ for some $a<\Lambda$, where $\Lambda=\pi j_{1,1}^{2}\approx 12.0377$ is the first buckling eigenvalue of the disk with unit area. Our numerical results suggest that actually the result may be true for a larger range of values of $\alpha$ and we conjecture that the ball is the minimizer for $\alpha\in[0,\alpha^{\star}]$. On the other hand, we have numerical evidence to support the conjecture that for $\alpha>\alpha^{\star}$, the ball is no longer the minimizer. For instance, for $\alpha=110$, the first eigenvalue of the ball of unit area can be directly calculated by solving (2.9) and is equal to -1622.16613… In Table 1 we show some numerical approximations for the first eigenvalue of the minimizer that we obtained with our algorithm when $\alpha=110$, which is plotted in Figure 3, for different values of $m$. These results suggest that the first eigenvalue of this domain is equal to -1786.35377…, which is significantly smaller than the first eigenvalue of the disk.

In Figure 5 we plot the eigenfunctions associated to the first three eigenvalues of the optimizers of $\lambda_{1}$, obtained for $\alpha=110,170,230$. In this work we considered just the optimization of the first eigenvalue. However, we observed that, besides the fact that the eigenfunction associated with the first eigenvalue changes sign, it also has different number of nodal domains, depending on the parameter $\alpha$. Moreover, ’similar’ eigenfunctions appear associated with eigenvalues of different orders. For instance, the eigenfunction associated with the first eigenvalue for $\alpha=110$ is antisymmetric with respect to the first axis. However, the eigenfunction associated with the first eigenvalue for $\alpha=170$ is symmetric with respect to the first axis and the first antisymmetric eigenfunction with respect to the first axis is associated not with the first eigenvalue, but with the second eigenvalue.

Figure 4 shows a zoom of the boundary of the optimizer obtained numerically for $\alpha=110$, in a neighbourhood of the re-entrant part of the boundary. Note that the boundary of the domains considered in the optimization procedure was parameterized by a (truncated) Fourier expansion. In particular the domains considered are always smooth and it is not clear how to obtain information on the regularity of the boundary of the optimizer from this. In particular, it is not possible to deduce whether this corresponds to a smooth boundary, a corner, or even a cusp.

Acknowledgements

This work was partially supported by the Fundação para a Ciência e a Tecnologia (Portugal) through the program “Investigador FCT” with reference IF/00177/2013 and the project Extremal spectral quantities and related problems (PTDC/MAT-CAL/4334/2014). Most of the research in this paper was carried out while the second author held a post-doctoral position at the University of Lisbon within the scope of this project. The second author is a member of the Gruppo Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM).