Dirichlet Eigenvalues (100-digits) of the Bunimovich mushroom
Robert S Jones
rsjones7 at yahoo dot com
August 22, 2016 (completed)
Posted below are the lowest few Dirichlet eigenvalues of the Bunimovich mushroom with dimensions shown in Fig. 1, precise to at least 100 correct digits.
Fig 1. My Bunimovich mushroom. The shape is a unit square foot centered under an R=3/2 radius semicircle hat. This calculation was inspired entirely by [1] "Quantum mushroom billiards" by Barnett and Betcke (2007) arXiv:nlin/0611059v2. Very specifically, I extend their numerical calculations of the lowest eigenvalues as displayed in their Table I to at least 100 digits. The blue digits in my Table 1 line up with and match their ten-digit results.
The set of eigenvalues are the lowest values of $\lambda$ such that $$\Delta\Psi({\mathbf r})+\lambda\Psi({\mathbf r})=0$$ where $\Psi=0$ on the indicated mushroom boundary. See below for details.
The solutions [eigenmodes] are classified as even or odd according to the line of symmetry that divides the mushroom in half, and so I work in the de-symmetrized half mushroom.
By carefully tuning the parameters (i.e., distribution of matching points), I am able to complete the calculations (all 20 eigenvalues) within one day elapsed time on one laptop computer (i7 quad core 8 thread w/ 16GB RAM). Without other interfering processes, each eigenvalue takes about one hour. The reported digits are all correct digits, not rounded.
Table 1. Mushroom eigenvalues. Even and odd mushroom modes are designated "ME" and "MO", respectively. A link to the calculation output is provided. n is the overall Dirichlet eigenvalue order, as far as I can determine from Ref [1].
n Mode $\lambda$ 1 ME-01 3 ME-02 4 ME-03 5 ME-04 8 ME-05 9 ME-06 10 ME-07 13 ME-08 15 ME-09 17 ME-10 2 MO-01 6 MO-02 7 MO-03 11 MO-04 12 MO-05 14 MO-06 16 MO-07 ? MO-08 ? MO-09 ? MO-10 General method of solution
I believe the following is the only practical method (using a typical laptop) to obtain significantly more than a dozen or so digits in the eigenvalue for this problem.
Refer to the Fig. 2, which shows the de-symmetrized mushroom, i.e., the fundamental region in which I work. In this diagram, I show the coordinates and a typical arrangement of matching points in the chosen proportion, as described below.
Fig 2. Geometry of fundamental region. Note that there are matching points close to the two vertices (at either end of the line of symmetry), but not on them. On $\partial\Omega_2$, $\partial\Omega_3$, and $\partial\Omega_4$, there are 5, 18, and 17 matching points, respectively. To begin the solution, expand the eigenfunction about the re-entrant corner with the usual Fourier-Bessel series, $$ \Psi(r,\theta)=\sum_{i=1}^{N}c_i J_{m_i}(kr)\sin(m_i\theta) $$ where $m_i=2i/3$, $k^2=\lambda$, and $\theta$ is adjusted to vary counter-clockwise from 0 to $3\pi/2$. $N$ is both the number of terms in this equation, and the number of matching points on the part of the boundary of the half-mushroom including ($\partial\Omega_2$) half the foot bottom, ($\partial\Omega_3$) the line of symmetry, and ($\partial\Omega_4$) the quarter-circle curve, with total length $$S=\frac12+\frac52+\frac{3\pi}{4}=3(1+\pi/4)\approx 5.356$$
The distribution of boundary matching points is very important. Evenly-spaced matching points give very poor convergence rates, about 1/50 (one extra eigenvalue digit to every 50 added matching points). It is well known that crowding near the distant vertices yields much better convergence rates. But crowding near the vertices at the ends of the adjacent edges (intersecting the re-entrant angle) is not necessary. With large N-values, this can be quite helpful in speeding up the computations.
I choose Chebyshev-like distributions and obtain convergence rates on the order of better than 1/10. At one hundred digits, the global convergence rate is about one digit per 8.5 matching points, and slowly worsening.
Along the foot bottom ($y=-1$), I use specifically: $$ x_j = \frac12 \sin \left( \frac{(j-\frac12)\pi}{2N_f} \right) $$ where $j=1$ to $N_f$ and where $N_f$ is the number of matching on that segment. These points are crowded near (1/2,-1) but not at (0,-1).
Along the symmetry line ($x=1/2$), they are spread out according to: $$ s_{j} = \frac12 \left[ 1 - \cos\left(\frac{(j-\frac12)\pi}{N_s}\right)\right]$$ where $j=1$ to $N_s$ and where $N_s$ is the number of matching on the symmetry line. These points are crowded at both ends.
Along the curve (quarter circle), they are spread out similar to the foot with the crowding at the top of the hat but not at its base.
By manually adjusting the proportion of matching points, I found that the ratio 5:18:17 (foot:symmetry-line:curve) yields a good convergence rate. I'm pretty sure there is a better optimal proportion, but this is close enough.
For comparison, as shown in the Table 2, two other proportions might be (2) to make them proportional to the length of the boundary portion (e.g., as in equally-spaced points) or (3) proportional to the subtended angle of the boundary portion (e.g., as in equal-angle points). Table 2 compares them to (1) my chosen proportion by rescaling to about 100 matching points.
Table 2. Comparison of proportions for different configurations of 100 matching points. The proportion order is "half of foot base" to "line of symmetry" to "quarter circle of half hat". proportion method 1 12.5 : 45.0 : 42.5 my chosen proportion (5:18:17) 2 9.3 : 46.7 : 44.0 "equal-spaced" 3 9.8 : 50.0 : 40.2 "equal-angular" The paper (Ref [1]) mentioned above does not specify the number of or the distribution of boundary matching points. (However, there are indications within the MPSpack and vergini source-code that equal-spaced points were used.) I think the convergence rate of all boundary methods (this method, GSVD, scaling) is improved when crowding points near the vertices, but the maximum gap between them should be a little less than half the free wavelength $\Lambda=2\pi/\sqrt{\lambda}$. This leads to a tradeoff for very high eigenvalues, as they had calculated.
This is indeed a critical aspect of this project: By carefully choosing the matching points (by trial-and-error), I was able to turn a difficult computation into a one-hour computation. Indeed, equally-spaced points, which still yields exponential convergence, may require a 5000x5000 matrix for a 100-digit computation. I can estimate that such a calculation would take about one month to compute (with my computer). Indeed, each time I increase N by 20%, the time doubles. As it is, my biggest matrix was 924x924, and this one matrix calculation takes about 30 minutes.
With a given N-value, and a set of N boundary matching points $\{ {\mathbf r}_i\}$, I construct the NxN point-matching matrix $M$. The matrix is considered a function of a free [eigen-]parameter $k=\sqrt{\lambda}$. For the odd modes, it is $$ M_{i,j}(k) = J_{m_i}(kr_j)\sin(m_i\theta_j) $$ where both $i$ and $j$ go from 1 to N. The determinant of this matrix must be zero for there to exist a non-trivial solution (in terms of the expansion coefficients).
For the even modes, we can use the same formula above for the matrix elements on the foot bottom and the curved portion. On the symmetry line, we must set $\partial\Psi/\partial x=0$. This leads to matrix elements on that symmetry line in the form $$ M_{i,j}(k) = m_i\left(\frac12 \sin(m_i\theta_j)-y_j\cos(m_i\theta_j)\right) J_{m_i}(kr_j) - \frac12 kr_j \sin(m_i\theta_j) J_{m_i+1}(kr_j) $$ where the $x_j=1/2$ comes from the width of half the foot. Remember in these matrix elements, $0 < \theta_j < 3\pi/2$, and the origin is at the re-entrant corner.
Of course, the next step is to find the root of the determinant of that point-matching matrix by varying the eigen-parameter, which I call the "approximate eigenvalue" for that N-value. Getting the estimates to alternate with increasing N is not reliable or requires special increments unique to each eigenvalue, so instead, I use the simpler fact that as N increases, the approximate eigenvalue gets closer and closer to the true eigenvalue.
I start with $N=40$ and increment by (about) 20%, each time picking up about 20% more eigenvalue digits. With each estimate, I guess a bracket (i.e., very conservative upper and lower bounds - much larger than they have to be) for starting the next iteration. Specifically, that conservative bracket is $$ L_{\left({\mbox{up}}\atop{\mbox{down}}\right)}=L_{N_i}\pm 2 \left| L_{N_i}-L_{N_{i-1}}\right| $$ I stop when I get at least 100 matching digits in $L_{N_i}$ and $L_{N_{i-1}}$.
To address the ill-conditioning, I find that working with real precision of about $0.75N$ is sufficient.
ME-01: Lowest eigenvalue to over 200 digits
Going a little further, letting the program run for a whole day on the lowest eigenvalue yields 221 correct digits:
References
[1] "Quantum mushroom billiards" by Barnett and Betcke (2007). arXiv:nlin/0611059v2.
