A negative result for hearing the shape of a triangle

A negative result for hearing the shape of a triangle A computer-assisted proof Gerard Orriols Gim´enez Joint Bachelor’s thesis for the degrees of Mathematics and Engineering Physics Supervisor (Princeton University): Javier G´omez-Serrano Supervisor (UPC): Xavier Cabr´e Universitat Polit`ecnica de Catalunya May 2019 Abstract We prove that there exist two distinct triangles for which the first, second and fourth eigenvalues of the Laplace operator with null Dirichlet boundary conditions coincide. This solves a conjecture raised by Antunes and Freitas and suggested by their informal numerical evidence. We use a novel technique for a computer-assisted proof about the spectrum of an operator, which combines a Finite Element Method, to locate roughly the first eigenvalues keeping track of their position in the spectrum, and the Method of Fundamental Solutions, to get a much more precise bound of these eigenvalues. Due to the time constraints, some of the computations still remain to be finished. Keywords— computer-assisted proof, Laplace eigenvalues, spectral geometry, Finite Element Method, Method of Fundamental Solutions. 1 Contents 1 Introduction 3 2 Structure of the proof of Theorem 1 4 3 Separation of the first four eigenvalues 6 4 Rigorous eigenvalue bounds for individual triangles 8 5 Extension of the bounds to a region of triangles 10 6 Implementation 13 6.1 Givensrotations ........................................... 13 6.2 Upper bound of the MFS boundary norm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 6.3 Lower bound of the MFS interior norm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Acknowledgments 15 A Code for the separation of the smallest eigenvalues 17 B Code for validating an interval of triangles 23 C Code for optimizing the eigenvalue 37 2 1 Introduction In this thesis, we present a computer-assisted proof of the following original theorem. Theorem 1. The first, second and fourth eigenvalues of the Laplace operator on an Euclidean triangle with null Dirichlet boundary conditions are not enough to determine it up to isometry. This is a conjecture proposed by Antunes and Freitas in [AF11], suggested by numerical evidence, but a rigorous proof was required. The Dirichlet eigenvalues of the Laplace operator for a triangle Ω are real numbers λsuch that there is a nonzero smooth function udefined on Ω and continuous on Ω such that (−∆u=λu in Ω u= 0 on ∂Ω. It is a well known fact that the set of such λforms an increasing sequence 0 < λ1< λ2≤λ3≤ ··· whose only limit point is ∞, an that the corresponding eigenfunctions ujform an orthonormal basis of L2(Ω). The eigenvalues of a domain are closely related with its geometric properties, constituting an active area of research called spectral geometry. A classical example of this relationship is Weyl’s law, which relates the asymptotics of the eigenvalues to the volume of the domain, and a later result by McKean and Singer states that the perimeter is also determined by the eigenvalues [MS67]. More results of this kind can be found in [BG90] and [vdBS88]. Other results about how the geometry of a domain determines its spectrum can be found in Henrot’s book [Hen06]. The question of the determination of a domain given the set of its Laplace eigenvalues was posed by Mark Kac in his famous paper “Can one hear the shape of a drum?” [Kac66]. Since then, the answer has been found to be negative in general; in particular, for euclidean polygons, the first example of a pair of non-isometric polygons with the same spectrum is due to Gordon, Webb and Wolpert [GWW92]. However, there are positive results when we restrict the determination to a class of domains, the most successful of which was found by Zelditch [Zel09], who proved spectral determination for analytic domains with two classes of symmetries. Less is known about domains with less regular boundaries, the simplest of which are polygons. In the case of triangles, it has been proven that the whole spectrum of the Laplace operator determines the shape of a triangle ([Dur88], with a recent simple proof by [GM13]), and later Chang and De Turck proved that only a finite amount of eigenvalues, which depends on λ1and λ2, is enough [CD89]. It is natural to try to improve the result to only a finite and fixed amount of eigenvalues, answering the question “Can a human hear the shape of a triangular drum?”. Since the space of triangles up to isometries is has dimension 3, we would expect that 3 eigenvalues should be enough, but a priori it is not clear which ones. Antunes and Freitas [AF11] conjectured that indeed the three first eigenvalues λ1, λ2and λ3do determine the shape of a triangle. Instead, numerical evidence by themselves seems to indicate that this is not the case for λ1, λ2and λ4, and in this paper we will prove this fact (Theorem 1). This will give an example of an obstruction to determining the shape of a triangle from a finite portion of its spectrum. 3 2 Structure of the proof of Theorem 1 By the scaling of the problem, we reduce our search to the set of triangles with a fixed base length (together with additional conditions that ensure that we only consider one triangle for each similarity class); instead of looking for all three eigenvalues λ1, λ2, λ4to be equal, we just require the quotients ξ21 =λ2/λ1and ξ41 =λ4/λ1to take the same value. Since the eigenvalues scale by λr−2when the lengths of a triangle are scaled by r, the two conditions are equivalent. Fixing the first two vertices of the triangle to be (0,0) and (1,0), we use the coordinates (cx, cy) of the third vertex to parametrize the search space. Our approach consists in using a topological argument to show that in each of two disjoint regions in this parameter space there is a triangle in which ξ21 and ξ41 take the same prescribed value. More precisely, we claim that there are two distinct triangles for which ξ21 =¯ ξ21 := 1.67675 and ξ41 =¯ ξ41 := 2.99372. Since rigorous calculations with the computer are done using interval arithmetic, we need a topological technique to transform the closed condition into an open condition which can be compatible with the intervals. For this we will use the Poincar´e–Miranda theorem (see [Mir40]): Theorem 2. Given two continuous functions f, g : [−1,1]2→Rsuch that f(x, y)has the same sign as x when x=±1and g(x, y)has the same sign as ywhen y=±1, there exists a point (x, y)∈[−1,1]2such that f(x, y) = g(x, y)=0. The regions that we will consider are two parallelograms around the points A= (0.63500,0.27500) and B= (0.84906,0.31995), designed such that ξ21 and ξ41 have approximately a constant value each in a pair of opposite edges. Using the computer we will verify that the functions ξ21 −¯ ξ21 and ξ41 −¯ ξ41 each have a constant and opposite sign in opposite edges of the parallelogram, and hence by the theorem, together with the well known continuity of eigenvalues, we will conclude that such two distinct triangles exist. The vectors defining the parallelogram are of course obtained from the inverse of an approximation of the differential of the R2-valued funcion (ξ21, ξ41) at the points Aand B, reescaled so that the increment of ξ21 and ξ41 along each vector is approximately 0.01. Explicitly, they are given by: •v21,A = (0.04269082311,0.01148489707) •v41,A = (−0.02082984105,−0.00255790585) •v21,B = (0.01717048015,0.01266844996) •v41,B = (0.03590363291,0.01065148385) Moreover, for the point Ain the negative v21,A direction, the vector is shrinked by a factor 0.6, because after that the behavior deviates from linear and the bound becomes worse after this displacement. Hence the vertices of the parallelograms will be •A+v21,A +v41,A,A+v21,A −v41,A,A−0.6v21,A −v41,A,A−0.6v21,A +v41,A. •B+v21,B +v41,B,B+v21,B −v41,B,B−v21,B −v41,B,B−v21,B +v41,B. 4 0.56 0.58 0.6 0.62 0.64 0.66 0.68 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.25 0.26 0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 A B ξ21 ξ41 Figure 1: Numerical approximate plot of the quotients ξ21 (discontinuous lines) and ξ41 (continuous lines) around the region of interest. The points Aand Band the validated parallelograms are shown in red. 5 This setup is displayed in Figure 1 together with a plot of non rigorous contour lines of the eigenvalue quotients. The pointwise verification of the values ξ21 and ξ41, which depends on an accurate calculation of λifor i= 1,2,4, consists of two steps. The first one, treated in Section 3, is about showing that the computed eigenvalues actually correspond to the ordered ones λ1, λ2and λ4; in order to do that, we will prove a lower bound for λ5combining techniques from the Finite Element Method with rigorous bounds linking the finite dimensional problem to the infinite dimensional one. The second step consists of finding accurate values of four eigenvalues that lie below the threshold obtained in the first part, which implies that they will indeed have to be the four lowest ones. This is done using the Method of Fundamental Solutions and recent rigorous bounds based on the L2norm of the boundary error, explained in Section 4. We emphasize the difficulty of finding the order of an eigenvalue, which is a global problem, compared to the local easier task of refining its value. To the best knowledge of the author, this is the first computerassisted proof in which these two distinct, local and global methods are used to verify eigenvalues of an operator. In order to extend the computer-assisted pointwise verifications of the eigenvalues λ1, λ2and λ4to the continuous region in which they need to be validated, we will use explicit continuity results of the spectrum of an elliptic operator with the domain to extend the bounds to a neighborhood of the verified points. The details of this part are explained in Section 5. Section 6 contains the details of the implementation and execution of the proof. Finally, the Appendix contains the codes used to verify an upper bound for the 4th eigenvalue at a point and to validate a segment of triangles. We end this section by introducing the reader to how a computer-assisted proof works. In the recent years, the application of calculations done by computers to mathematical proofs have become more popular due to the increment of computational resources, but in order to make sure that their results are rigorous, we need to control the errors that floating point arithmetic can accumulate. This is usually done by means of interval arithmetic, in which the data that a computer stores for a real number is an interval (two endpoints, or a midpoint and a radius) of real numbers, stored by two floating point numbers, instead of just one. Operations between intervals are implemented to return intervals which are guaranteed to contain every possible result when the operands belong to the input intervals. For example, if [x]=[x, x] and [y]=[y,y] are two intervals, their sum will can be given by the interval [x] + [y]=[x+y, x +y] and their product by [x]·[y] = [min{xy, xy, xy, xy},max{xy, xy, xy, xy}]. The same rule applies to function implementations: a function fevaluated on [x] should return an interval containing every f(x) for x∈[x]. We refer to the book [Tuc11] for an introduction to validated numerics, in which most of the techniques used here are explained, and to [GS18] for a more specific treatment of computer-assisted proofs in PDE. 3 Separation of the first four eigenvalues In order to find a rigorous lower bound for the fifth eigenvalue of a triangle we will use a recent bound found by Liu [Liu15], which is similar to the one in [CG14] but simplifies the hypotheses and improves the constant. Both use the non-conforming Finite Element Method of Crouzeix–Raviart; other rigorous bounds with conforming finite elements were explored, like [LO13], but the bound is worse and the method is harder 6 to implement with validated numerics because its mass matrix is not diagonal. The Crouzeix–Raviart finite-element method uses a triangulation of the domain Ω, which in our case we will take to be the trivial triangulation given by N2triangles with sides equal to 1/N of the original one and similar to it. The basis functions are indexed by interior edges of the triangulation: if Eis a common edge of triangles T1,T2, the basis function ψEis the unique function supported on T1∪T2such that restricted to each triangle is affine, takes the value 1 in the midpoint of Eand the value 0 in the midpoints of the other edges of T1and T2. We define the coefficients of the stiffness and mass matrices A= (aEF ),B= (bEF ) by the bilinear forms aEF =ZΩ∇ψE·∇ψFbEF =ZΩ ψEψF For our choice of triangulation, Bis simply a multiple of the identity 2|Ω|I/3N2, whereas Ais a sparse matrix. This will allow us to work with a matrix eigenvalue problem instead of a generalized one. The main result that we will use is the following [Liu15, Theorem 2.1 and Remark 2.2]: Theorem 3. Consider a polygonal domain Ωwith a triangulation so that each triangle has diameter at most h. Let λkbe the k-th eigenvalue of Ωand λk,h the k-th eigenvalue of the Crouzeix–Raviart discretized problem for Ω. Then λh,k 1 + C2 hλh,k ≤λk,(1) where Ch≤0.1893his a constant. In order to be able to deal with approximate eigenvalues we will need in addition the following lemma from [Par80, Theorem 15.9.1]. Lemma 4. Let (˜ λh,˜ uh)be an approximated algebraic eigenpair such that ˜ λhis closer to some λhthan to any other discrete eigenvalue. Suppose that the coefficient vector ˜ uhis normalised with respect to B, kB˜ uhkB−1=k˜ uhkB= 1. Then the algebraic residual r:= A˜ uh−˜ λhB˜ uhsatisfies |λh−˜ λh| ≤ krkB−1. Remark 5. We can combine Theorem 3 with Lemma 4 using the monotonicity of (1), using λh,k −krkB−1 as a lower bound of λh,k instead. It is easy to obtain estimations ˜ λhwith a very small residual. The hardest part here before applying the theorem is to check that they have indeed the correct index, i.e., that they are closer to the appropriate λhthan to any other discrete eigenvalue. This is why we need to control the whole spectrum of the discrete problem, and we will do that by applying Gershgorin’s disks theorem after performing some Givens rotations to the matrix to get useful bounds. More precisely, in order to get a lower bound for λ5we need to separate the first 5 eigenvalues from the rest so that we can control them, and in order to do that we will perform Givens rotations until the intervals provided by Gershgorin’s theorem can be separated in two disjoint components, one containing the 5 smallest eigenvalues and the other containing the rest. If this holds, then the strong version of Gershgorin’s 7 theorem will guarantee an upper bound for λh,5. Therefore, provided that the residuals are all very small and that all approximate eigenvalues are different (which happens in our setting), Lemma 4 will guarantee that there are 5 distinct discrete eigenvalues below the upper bound and therefore they will be forced to have the correct indices. This allows us to verify λh,5with an error only depending on its residual, and using Remark 5, get a lower bound of λ5. Therefore, the first four eigenvalues can be separated just by checking that they are distinct and smaller than this lower bound. 4 Rigorous eigenvalue bounds for individual triangles Our approach to find tight bounds for the eigenvalues of triangles uses the Method of Fundamental Solutions (MFS), introduced by Fox, Henrici and Moler in [FHM67] and more recently revived by Betcke and Trefethen [BT05]. In this method, a function uis written as a linear combination of functions φi(1 ≤i≤N) that satisfy pointwise the equation (∆ + λ)φi= 0 for a fixed λ. The coefficients are chosen to optimize the proximity of the function to the eigenspace for the actual eigenvalue λj, in a sense made precise in [BT05], and this is measured by the least singular value of a certain matrix that involves the values of uat discrete points of the boundary ∂Ω. This parameter is minimized with respect to λby using a golden ratio search. This provides a candidate λ∈Rand coefficients cifor which u(x) = PN i=1 ciφi(x) can be computed with arbitrary precision. The functions that we will use for the MFS consist of two types: the first ones are of the form φ(x) = Y0(√λ|x−x0|), for x0a point outside Ω. The second type of functions are parametrized by a vertex of the triangle and a positive integer j, and take the form ψj(r, θ) = Jjα(√λr) sin(jαθ), where (r, θ) are the polar coordinates of the point with respect to a vertex in the triangle whose total angle is π/α, and θis measured from an adjacent side. The first kind of functions allow us to approximate the function in the interior of the triangle and near the sides, while the second kind gives the correct asymptotic behavior of the solution near the vertices of the triangle. The main tool that we will use to find rigorous bounds for eigenvalues is the L2bound given by Barnett and Hassell [BH11]. However, their method is optimized for high eigenvalues, so we will have to adapt some of the steps to our case of small eigenvalues. We summarize the main results that we will use. Let Ω be a triangle and u∈C2(Ω) be nonzero such that (∆ + λ)u= 0. Consider the tension t[u] = kukL2(∂Ω) kukL2(Ω) . Let λj, ujbe the sequence of eigenvalues and eigenfunctions of Ω, satisfying (∆ + λj)uj= 0 with Dirichlet null boundary conditions. Let vjbe the normal derivative of uj, defined on ∂Ω. We define the operator A(λ) = X λj vjhvj,·i (λ−λj)2, and its decomposition as a sum of three: Anear(λ) = X |λ−λj|≤√λ vjhvj,·i (λ−λj)2, 8 −20,000 −10,000 0 10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000 Figure 2: Gershgorin intervals before applying the Givens rotations for triangle A. −20,000 −10,000 0 10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000 Figure 3: Gershgorin intervals after applying the Givens rotations for triangle A. The dashed vertical line at position 790 separates the first 5 eigenvalues from the rest. Acknowledgments First of all, I would like to thank my supervisor Javier G´omez-Serrano for hosting me at Princeton University, suggesting me the problem, providing me references, giving me a lot of guidance and sharing with me many 15 −200 0 200 400 600 800 1,000 1,200 1,400 1,600 1,800 Figure 4: Zoom of the separation of the first 5 eigenvalues after the Givens rotations. (a) (b) Figure 5: Grid used to validate a lower bound for kukL2(Ω) for triangle B, with (a) the first eigenfunction and (b) the second eigenfunction plotted on top. useful discussions. I am very grateful to the Mathematics Department of Princeton University, in particular the Graduate Program, for funding my university fee and allowing my stay as a visitor researcher. I would also like to thank CFIS and the administrators of the Mobility program for giving me the oportunity to do this research in another university and funding me. I also thank the MOBINT scholarship for providing partial financial support for my stay. 16 A Code for the separation of the smallest eigenvalues This program reads the parameter Nand the coordinates of the third vertex of the triangle and returns a lower bound for λ5. #include "arb.h" #include "arb_mat.h" #include "linalg.h" #include "arbcc.h" #include <cassert> #include <cmath> #include <iostream> #include <map> #include <vector> // Compile with -DPREC=1024 using namespace std; using namespace alglib; int di[6] = {0, 0, 1, -1, 1, -1}; int dj[6] = {1, -1, 0, 0, -1, 1}; int side(int i, int j) { if (i % 2 == 1) return 2; if (j % 2 == 0) return 0; return 1; } double sq(double t) { return t*t;} /* struct ArfReal definition omitted for brevity */ int n, m; void givens(vector<vector<ArbReal> >& A, double lowbar, double upbar, int nsmall) { vector<pair<ArbReal, int> > srt(m); for (int i=0;i<m;++i){ srt[i] = make_pair(A[i][i], i); 17 } sort(srt.begin(), srt.end()); vector<bool> small(m, false); for (int i = 0; i < nsmall; ++i) { small[srt[i].second] = true; } bool done = false; int it = 0; while (!done and it < 10) { ++it; done = true; for (int p=0;p<m;++p){ ArbReal rad = 0.; for (int q=0;q<m;++q){ if (q != p) { rad = rad + A[q][p].abs(); } } ArbReal margin; if (small[p]) { margin = ArbReal(upbar) - A[p][p]; }else { margin = A[p][p] - ArbReal(lowbar); } if (rad < margin) { continue; } done = false; printf("it %d: changing column %d with margin ", it, p); margin.print(); for (int q=0;q<m;++q){ if (q != p && margin / (m - 1) < A[p][q].abs()) { ArfReal tt = (ArfReal(A[q][q]) - ArfReal(A[p][p])) / (ArfReal(A[p][q]) * ArfReal(2)); ArfReal t = (tt).sign() / (tt.abs() + (tt.sq() + 1).sqrt()); ArbReal t_arb(t.x); ArbReal c = 1_a / (t_arb.sq() + 1_a).sqrt(); ArbReal s = c * t_arb; for (int i=0;i<m;++i){ ArbReal u = A[p][i]; 18 ArbReal v = A[q][i]; A[p][i] = c * u - s * v; A[q][i] = c * v + s * u; } for (int i=0;i<m;++i){ ArbReal u = A[i][p]; ArbReal v = A[i][q]; A[i][p] = c * u - s * v; A[i][q] = c * v + s * u; } } } } } } vector<pair<ArbReal, ArbReal> > gershgorin(vector<vector<ArbReal> >& A) { vector<pair<ArbReal, ArbReal> > ret(m); for (int i=0;i<m;++i){ ret[i].first = A[i][i]; ret[i].second = 0.; for (int j=0;j<m;++j){ if (j != i) { ret[i].second = ret[i].second + A[j][i].abs(); } } } sort(ret.begin(), ret.end()); return ret; } int main() { ios::sync_with_stdio(false); cout.setf(ios::fixed); cout.precision(10); cin >> n; double x_s, y_s; cin >> x_s >> y_s; ArbReal x(x_s), y(y_s); 19 ArbReal area = y / 2_a; ArbReal l_a = (y * y + (x - 1_a) * (x - 1_a)).sqrt(); ArbReal l_b = (y * y + x * x).sqrt(); ArbReal l_c = 1; ArbReal len[3]; len[0] = l_a; len[1] = l_b; len[2] = l_c; ArbReal h_a = y / l_a; ArbReal h_b = y / l_b; ArbReal h_c = y; ArbReal height[3]; height[0] = h_a; height[1] = h_b; height[2] = h_c; ArbReal cos_a = x / l_b; ArbReal cos_b = (1_a - x) / l_a; ArbReal cos_c = (x * (x - 1_a) + y * y) / (l_a * l_b); ArbReal cosine[3]; cosine[0] = cos_a; cosine[1] = cos_b; cosine[2] = cos_c; map<pair<int,int>, int> encode; for (int i=0;i<2*n-2;++i){ for (int j=0;i+j<2*n-2;++j){ if (i%2==0or j%2==0){ encode[make_pair(i, j)] = m++; } } } vector<vector<ArbReal> > mat(m, vector<ArbReal>(m)); ArbReal mass_diag = (area * 2) / (3 * n * n); for (int i=0;i<2*n-2;++i){ for (int j=0;i+j<2*n-2;++j){ 20 if (i%2==1and j%2==1)continue; int u = encode[make_pair(i, j)]; int side_u = side(i, j); mat[u][u] = (area * 8) / height[side_u].sq(); for (int d=0;d<6;++d){ int i2 = i + di[d]; int j2 = j + dj[d]; if (i2 >= 0 and j2 >= 0 and i2+j2<2*n-2and (i2 % 2 == 0 or j2 % 2 == 0)) { int v = encode[make_pair(i2, j2)]; int side_v = side(i2, j2); ArbReal tmp = -(cosine[3 - side_u - side_v] * area * 4) / (height[side_u] * height[side_v]); mat[u][v] = tmp; mat[v][u] = tmp; } } } } for (int i=0;i<m;++i){ for (int j=0;j<m;++j){ mat[i][j] = mat[i][j] / mass_diag; } } real_2d_array mat_d; mat_d.setlength(m, m); for (int i=0;i<m;++i){ for (int j=0;j<m;++j){ mat_d(i, j) = mat[i][j].get_approx_double(); } } ae_int_t num_eigs_found; real_1d_array eigs_d; real_2d_array eigenvectors; 21 assert(smatrixevdr(mat_d, m, 1, 1, 0.0, 1000.0, num_eigs_found, eigs_d, eigenvectors)); cerr << "Found " << num_eigs_found << " below 1000\n"; vector<ArbReal> discrete_eigs_lb(num_eigs_found); for (int i = 0; i < num_eigs_found; ++i) { cout << eigs_d(i) << endl; ArbReal lambda(eigs_d(i)); vector<ArbReal> vec(m); for (int j=0;j<m;++j){ vec[j] = ArbReal(eigenvectors(j, i)); } ArbReal num = 0., den = 0.; for (int j=0;j<m;++j){ ArbReal residue = -lambda * vec[j]; for (int k=0;k<m;++k){ residue = residue + ArbReal(mat[j][k]) * vec[k]; } num = num + residue.sq(); den = den + vec[j].sq(); } ArbReal error = (num / den).sqrt(); printf("Error: "); error.print(); discrete_eigs_lb[i] = lambda - error; } assert(num_eigs_found >= 6); double barrier = (eigs_d[4] + eigs_d[5]) / 2; // Successively apply Givens rotations to reduce the margin. // It does not work if one tries to separate them with the // barrier from the beginning. givens(mat, -100000, 1000000, 10); givens(mat, -80000, 1000000, 10); givens(mat, -60000, 1000000, 10); givens(mat, -50000, 1000000, 10); givens(mat, -20000, 1000000, 10); givens(mat, -10000, 1000000, 10); givens(mat, -1000, 100000, 10); 22 givens(mat, 0, 100000, 10); givens(mat, 1000, 10000, 10); givens(mat, 1000, 2000, 10); givens(mat, barrier, 2000, 5); givens(mat, barrier, barrier, 5); printf("Discrete lambda_5 successfully separated\n"); vector<pair<ArbReal, ArbReal> > ger = gershgorin(mat); ArbReal liu_constant = 0.1893; ArbReal lambda5lb = discrete_eigs_lb[4] / (1_a + discrete_eigs_lb[4] * liu_constant.sq() / ArbReal(n * n)); printf("Verified lower bound for lambda_5: "); lambda5lb.print(); } B Code for validating an interval of triangles This program reads the name of the third vertex (Aor B), the indexs of the eigenvalue quotient to validate (“21” or “41”), a letter standing for the sign (’p’ for plus, ’m’ for minus) of ξ−¯ ξthat we have to validate, the current coordinate c, and the total number of coordinates ncin which this side is divided (c∈ {1, . . . ,nc}). This will validate a segment around the position given by 2(c−0.5)/nc−1 in the corresponding side of the parallelogram (where the position is normalized so that vertices are at positions −1 and 1), and the radius `of this segment will be the value returned by the program. For example, if the input is “B 21 p 312 1000” and the program returns a value of 0.0013, this means that for all the points in the segment between B+v21,B −0.1898v41,B and B+v21,B −0.1872v41,B the corresponding ξ21 has been verified to be greater than ¯ ξ21. The program uses the ArbTaylor class, which is not included here for brevity, but computes simultaneously the Taylor polynomial and its residue, and uses them to obtain an enclosure of the function in an interval, as described in [Tuc11]. #include "arb.h" #include "arb_hypgeom.h" #include "arbcc.h" #include "arbseries.h" #include "arbtaylor.h" 23 #include <algorithm> #include <cassert> #include <fstream> #include <iostream> #include <map> #include <vector> // Compile with -DPREC=128 using namespace std; // Declared here and implemented in a separate file. void fill_in_vectors(vector<double>& r_coefs, double& r_lambda, vector<vector<double>>& r_sources, double cx0, double cy0, double lmin, double lmax); namespace { // Coordinates of the third vertex. Set as global because they are used // everywhere. ArbReal cx, cy; /* GENERAL FUNCTIONS */ // Fundamental solution around an external charge point ‘(xs, ys)‘. ArbTaylor fund_sol_charge(const ArbTaylor& x, const ArbTaylor& y, const ArbReal& xs, const ArbReal& ys, const ArbReal& lambda, int order) { ArbTaylor arg = (x - ArbTaylor::constant(xs, order)).sq() + (y - ArbTaylor::constant(ys, order)).sq(); arg = arg.sqrt(); arg *= lambda.sqrt(); return arg.bessely0(); } // Angle between two vectors. ArbReal ccw_angle(const ArbReal& x1, const ArbReal& y1, const ArbReal& x2, const ArbReal& y2) { ArbReal r = (x1.sq() + y1.sq()).sqrt(); ArbReal c = x1 / r; 24 ans = ArbReal(aux); arf_clear(aux); }else { ArbReal val1 = path_dfs(coef, lambda, sources, 2 * num, 2 * den, err_goal, xx, yy); ArbReal val2 = path_dfs(coef, lambda, sources, 2 * num + 1, 2 * den, err_goal, xx, yy); arb_min(ans, val1, val2, PREC); } mag_clear(err_val); return ans; } // Minimum of u^2 in the path given by (xx[i], yy[i]). ArbReal path_minsq(const vector<ArbReal>& coef, const ArbReal& lambda, const vector<vector<ArbReal>>& sources, double goal_d, const vector<ArbReal>& xx, const vector<ArbReal>& yy) { mag_t err_goal; mag_init(err_goal); mag_set_d(err_goal, goal_d); int den = 18; ArbReal lb; arb_pos_inf(lb); for (int num = 0; num < den; ++num) { ArbReal val = path_dfs(coef, lambda, sources, num, den, err_goal, xx, yy).sq(); arb_min(lb, lb, val, PREC); } mag_clear(err_goal); return lb; } ArbReal compute_l2meanlb(const vector<ArbReal>& coef, ArbReal lambda, const vector<vector<ArbReal>>& sources) { double goal = 5e-2; int nt = 8; // Bounding L2 norm using 8*8 triangles, ArbReal l2sum = 0_a; ArbReal fr = 0.8; // Using 0.64 of the total area to avoid the borders. 31 for (int counter = 0; counter < nt * nt; ++counter) { auto path = get_triangle(counter, nt, fr); ArbReal pmsq = path_minsq(coef, lambda, sources, goal, path.first, path.second); l2sum += pmsq * fr.sq() / ArbReal(nt).sq(); } ArbReal l2meanlb = l2sum.sqrt(); printf("Lower bound of the L2 mean: "); l2meanlb.print(); return l2meanlb; } /* COMPUTE THE EIGENVALUE WITH AN ABSOLUTE ERROR */ // Returns the pair {lambda, abs_err}, with lambda in [lbound, ubound]. ‘index‘ // specifies the position of the eigenvalue in the spectrum (1,2,...). pair<ArbReal, ArbReal> compute_all(int index, double lbound, double ubound) { vector<double> coef_d; vector<vector<double>> sources_d; double lambda_d; fill_in_vectors(coef_d, lambda_d, sources_d, cx.get_approx_double(), cy.get_approx_double(), lbound, ubound); int nc = coef_d.size(); int ns = sources_d.size(); vector<ArbReal> coef(nc); vector<vector<ArbReal>> sources(ns, vector<ArbReal>(2)); ArbReal lambda = lambda_d; for (int i = 0; i < nc; ++i) { coef[i] = ArbReal(coef_d[i]); } for (int i = 0; i < ns; ++i) { sources[i][0] = ArbReal(sources_d[i][0]); sources[i][1] = ArbReal(sources_d[i][1]); } printf("Candidate lambda: "); lambda.print(); 32 ArbReal meanl2lb = compute_l2meanlb(coef, lambda, sources); double goal_d = 1e-5; arf_t goal; arf_init(goal); arf_set_d(goal, goal_d); ArbReal l2bdry = side_l2_rec(coef, lambda, sources, goal); printf("Validated at goal %e\n", arf_get_d(goal, 10)); printf("L2 norm at the boundary: "); l2bdry.print(); ArbReal area = cy / 2; ArbReal totall2lb = meanl2lb * area.sqrt(); ArbReal tension = l2bdry / totall2lb; ArbReal side_a = (cy.sq() + (1_a - cx).sq()).sqrt(); ArbReal side_b = (cy.sq() + cx.sq()).sqrt(); ArbReal semip = (side_a + side_b + 1_a) / 2_a; ArbReal inrad = ((semip - side_a) * (semip - side_b) * (semip - 1_a) / semip).sqrt(); ArbReal c_omega = 4_a * (1_a + inrad) / inrad; ArbReal ub_all_lambda = 800_a; // Rough upper bound for all lambda (1 to 4). ArbReal a_far = 7_a * c_omega; ArbReal a_tail = 7_a * c_omega / lambda.sqrt(); ArbReal aux = 1_a / tension.sq(); aux -= a_far + a_tail; ArbReal abs_err = 2_a * ArbReal(ubound) / (aux * inrad); assert(0_a < abs_err); abs_err = abs_err.sqrt(); printf( "Rigorous value of lambda obtained. Follows lambda, its tension and its " "absolute error.\n"); lambda.print(); tension.print(); abs_err.print(); arf_clear(goal); return {lambda, abs_err}; 33 } }// namespace int main() { char ptname, signc; string indexs; int coord, ncoord; cin >> ptname >> indexs >> signc >> coord >> ncoord; ArbReal cx0, cy0, v21x, v21y, v41x, v41y, position; // Heuristics to find the appropritate eigenvalue. // This will be validated later by a separate program. vector<double> lmin, lmax; if (ptname == ’A’) { // Coordinates: cx0 = 0.635; cy0 = 0.275; // Vectors of the parallelogram: v21x = 4.269082311683548 / 100; v21y = 1.148489707350921 / 100; v41x = -2.082984105473002 / 100; v41y = -0.255790585210996 / 100; if (indexs == "21") { position = (ArbReal(coord) - 0.5_a) / ArbReal(ncoord) * 2_a - 1_a; lmin = {180, 320}; lmax = {280, 450}; }else { // In this case, the parallelogram is shrinked from below: // the position takes values in [-0.6, 1]. position = (ArbReal(coord) - 0.5_a) / ArbReal(ncoord) * 1.6_a - 0.6_a; lmin = {180, 620}; lmax = {280, 727. - (2 * coord - ncoord) / ncoord * 27}; } }else { cx0 = 0.849057346949971; cy0 = 0.319950941965592; v21x = 1.717048015465781 / 100; v21y = 1.266844996461298 / 100; 34 v41x = 3.590363291549854 / 100; v41y = 1.065148385561495 / 100; position = (ArbReal(coord) - 0.5_a) / ArbReal(ncoord) * 2_a - 1_a; if (indexs == "21") { if (signc == ’m’) { lmin = {160, 280}; lmax = {255, 380}; }else { lmin = {150, 280}; lmax = {190, 380}; } }else { if (signc == ’m’) { double center_l1 = 214 - coord * 23.0 / ncoord; double center_l4 = 637 - coord * 65.0 / ncoord; lmin = {center_l1 - 20, center_l4 - 20}; lmax = {center_l1 + 20, center_l4 + 20}; }else { double center_l1 = 197 - coord * 20.0 / ncoord; double center_l4 = 593 - coord * 60.0 / ncoord; lmin = {center_l1 - 20, center_l4 - 20}; lmax = {center_l1 + 20, center_l4 + 20}; } } } if (indexs == "21") { if (signc == ’m’) { cx = cx0 - v21x + position * v41x; cy = cy0 - v21y + position * v41y; }else { cx = cx0 + v21x + position * v41x; cy = cy0 + v21y + position * v41y; } }else { if (signc == ’m’) { cx = cx0 + position * v21x - v41x; cy = cy0 + position * v21y - v41y; }else { cx = cx0 + position * v21x + v41x; 35 cy = cy0 + position * v21y + v41y; } } string filename = "output_"; filename += ptname; filename += indexs; filename += signc; filename += to_string(coord); filename += "of"; filename += to_string(ncoord); filename += ".txt"; freopen(filename.c_str(), "w", stdout); auto p1 = compute_all(1, lmin[0], lmax[0]); auto pk = compute_all((indexs == "21" ? 2 : 4), lmin[1], lmax[1]); ArbReal xi = pk.first / p1.first; ArbReal errxi = (pk.second + p1.second * xi) / (p1.first - p1.second); printf("Xi and error:\n"); xi.print(); errxi.print(); ArbReal goalxi = (indexs == "21" ? 1.67675 : 2.99372); ArbReal margin = (signc == ’p’ ? xi - goalxi - errxi : goalxi - xi - errxi); ArbReal eps = margin / (pk.first * (1_a + xi + margin)); ArbReal delta = ((2_a * ArbReal::pi()).sq() * eps) / (1_a + (2_a * ArbReal::pi()).sq() * eps); ArbReal cc, kk; // Auxilliary variables to calculate ell. if (indexs == "21") { cc = cy * (v41x.abs() + 2_a * v41y.abs()); kk = v41y.abs(); }else { cc = cy * (v21x.abs() + 2_a * v21y.abs()); kk = v21y.abs(); } // Length validated at this point through propagation of the error. // Given as the coefficient of the vector v21 or v41 that one can move. ArbReal ell = (2_a * delta * cy * kk + cc - (4_a * delta * cy * kk * cc + cc.sq()).sqrt()) / 36 (2_a * delta * kk.sq()); if (ptname == ’A’ && indexs == "41") { // Renormalize ell so it represents a length in [-1, 1]: ell *= 1.25; } printf("Value of ell: "); ell.print(); } C Code for optimizing the eigenvalue This code implements the function fill_in_vectors declared in the file above, which finds an approximate eigenvalue and eigenvector in a non-rigorous fashion. It does so by optimizing the smallest nonzero singular value of a matrix that imposes the boundary conditions, as explained above. #include <stdio.h> #include <sys/time.h> #include <algorithm> #include <utility> #include "linalg.h" #include "stdafx.h" #include "optimization.h" #include "specialfunctions.h" #include <boost/math/special_functions/bessel.hpp> #include <boost/math/tools/minima.hpp> using namespace alglib; using namespace std; namespace { double cx, cy; inline double sq(double x) { return x*x;} double norm(const real_1d_array& a) { double s = 0; for (int i = 0; i < a.length(); ++i) { s += sq(a[i]); } 37 return sqrt(s); } real_1d_array operator*(const real_1d_array& a, double x) { real_1d_array c; int n = a.length(); c.setlength(n); for (int i=0;i<n;++i){ c[i] = a[i] * x; } return c; } real_1d_array operator/(const real_1d_array& a, double x) { real_1d_array c; int n = a.length(); c.setlength(n); for (int i=0;i<n;++i){ c[i] = a[i] / x; } return c; } real_1d_array operator+(const real_1d_array& a, const real_1d_array& b) { real_1d_array c; int n = a.length(); c.setlength(n); for (int i=0;i<n;++i){ c[i] = a[i] + b[i]; } return c; } real_1d_array operator-(const real_1d_array& a, const real_1d_array& b) { real_1d_array c; int n = a.length(); c.setlength(n); for (int i=0;i<n;++i){ c[i] = a[i] - b[i]; } 38 return c; } double fund_sol_charge_d(double x, double y, double x0, double y0, double lambda) { double r = sqrt(sq(x - x0) + sq(y - y0)); return bessely0(r * sqrt(lambda)); } double fund_sol_vertex_d(double x, double y, int v, int j, double lambda) { double xvert[3] = {0, 1, cx}; double yvert[3] = {0, 0, cy}; int nxt = (v + 1) % 3; int prv = (v + 2) % 3; double xv = xvert[v]; double yv = yvert[v]; double vnxt_[2] = {xvert[nxt] - xv, yvert[nxt] - yv}; double vprv_[2] = {xvert[prv] - xv, yvert[prv] - yv}; double vcur_[2] = {x - xv, y - yv}; real_1d_array vnxt, vprv, vcur; vnxt.setcontent(2, vnxt_); vprv.setcontent(2, vprv_); vcur.setcontent(2, vcur_); double thmx = atan2(vprv[1], vprv[0]) - atan2(vnxt[1], vnxt[0]); if (thmx < 0) { thmx += 2 * M_PI; } if (thmx > 2 * M_PI) { thmx -= 2 * M_PI; } double th = atan2(vcur[1], vcur[0]) - atan2(vnxt[1], vnxt[0]); if (th < 0) { th += 2 * M_PI; } if (th > 2 * M_PI) { th -= 2 * M_PI; } double r = norm(vcur); double alpha = M_PI * (j + 1) / thmx; double sl = sqrt(lambda); 39 return boost::math::cyl_bessel_j(alpha, sl * r) * sin(alpha * th); } // Generates source points: ‘ns‘ per side and an aditional ‘nextra‘ // situated around the top vertex. ‘nextra‘ should be odd. real_2d_array gen_sources(int ns, int nextra) { double beta = 0.01; real_2d_array sources; sources.setlength(3 * ns + nextra, 2); double xvert[3] = {0, 1, cx}; double yvert[3] = {0, 0, cy}; for (int v=0;v<3;++v){ real_1d_array vcur, vnxt; double vcur_[2] = {xvert[v], yvert[v]}; vcur.setcontent(2, vcur_); int nxt = (v + 1) % 3; double vnxt_[2] = {xvert[nxt], yvert[nxt]}; vnxt.setcontent(2, vnxt_); for (int i = 0; i < ns; ++i) { double rat = 0.5 * (1.0 - cos(double(i + 1) * M_PI / double(ns))); real_1d_array pt = vcur * (1.0 - rat) + vnxt * rat; real_1d_array vec = vnxt - vcur; vec = vec / norm(vec); double perp_[2] = {vec[1], -vec[0]}; real_1d_array perp; perp.setcontent(2, perp_); for (int c=0;c<2;++c){ sources(v * ns + i, c) = pt[c] + beta * perp[c]; } } } for (int i = 0; i < nextra; ++i) { sources(3 * ns + i, 0) = cx + (double(i) - (nextra - 1) / 2) * 0.02; sources(3 * ns + i, 1) = cy + 0.1; } return sources; } // Generates ‘ni‘ interior points at random. 40 [Par80] B. N. Parlett. The Symmetric Eigenvalue Problem. Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, 1980. [Rel40] F. Rellich. Darstellung der Eigenwerte von ∆u+λu = 0 durch ein Randintegral. Math. Z., 46:635–636, 1940. [Tuc11] W. Tucker. Validated numerics: a short introduction to rigorous computations. Princeton University Press, 2011. [vdBS88] M. van den Berg and S. Srisatkunarajah. Heat equation for a region in R2with a polygonal boundary. J. London Math. Soc. (2), 37(1):119–127, 1988. [Zel09] S. Zelditch. Inverse spectral problem for analytic domains. II. Z2-symmetric domains. Ann. of Math. (2), 170(1):205–269, 2009. 47