A nega i e esul o hea ing he
shape o a iangle
A compu e -assis ed p oo
Ge a d O iols Gim´enez
Join Bachelo ’s hesis o he deg ees o
Ma hema ics and Enginee ing Physics
Supe iso (P ince on Uni e si y): Ja ie G´omez-Se ano
Supe iso (UPC): Xa ie Cab ´e
Uni e si a Poli `ecnica de Ca alunya
May 2019
Abs ac
We p o e ha he e exis wo dis inc iangles o which he i s , second and ou h eigen-
alues o he Laplace ope a o wi h null Di ichle bounda y condi ions coincide. This sol es a
conjec u e aised by An unes and F ei as and sugges ed by hei in o mal nume ical e idence.
We use a no el echnique o a compu e -assis ed p oo abou he spec um o an ope a o ,
which combines a Fini e Elemen Me hod, o loca e oughly he i s eigen alues keeping ack
o hei posi ion in he spec um, and he Me hod o Fundamen al Solu ions, o ge a much
mo e p ecise bound o hese eigen alues. Due o he ime cons ain s, some o he compu a ions
s ill emain o be inished.
Keywo ds— compu e -assis ed p oo , Laplace eigen alues, spec al geome y, Fini e Elemen Me hod,
Me hod o Fundamen al Solu ions.
1
Con en s
1 In oduc ion 3
2 S uc u e o he p oo o Theo em 1 4
3 Sepa a ion o he i s ou eigen alues 6
4 Rigo ous eigen alue bounds o indi idual iangles 8
5 Ex ension o he bounds o a egion o iangles 10
6 Implemen a ion 13
6.1 Gi ens o a ions ........................................... 13
6.2 Uppe bound o he MFS bounda y no m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
6.3 Lowe bound o he MFS in e io no m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Acknowledgmen s 15
A Code o he sepa a ion o he smalles eigen alues 17
B Code o alida ing an in e al o iangles 23
C Code o op imizing he eigen alue 37
2
1 In oduc ion
In his hesis, we p esen a compu e -assis ed p oo o he ollowing o iginal heo em.
Theo em 1. The i s , second and ou h eigen alues o he Laplace ope a o on an Euclidean iangle wi h
null Di ichle bounda y condi ions a e no enough o de e mine i up o isome y.
This is a conjec u e p oposed by An unes and F ei as in [AF11], sugges ed by nume ical e idence, bu
a igo ous p oo was equi ed. The Di ichle eigen alues o he Laplace ope a o o a iangle Ω a e eal
numbe s λsuch ha he e is a nonze o smoo h unc ion ude ined on Ω and con inuous on Ω such ha
(−∆u=λu in Ω
u= 0 on ∂Ω.
I is a well known ac ha he se o such λ o ms an inc easing sequence 0 < λ1< λ2≤λ3≤ ··· whose
only limi poin is ∞, an ha he co esponding eigen unc ions uj o m an o hono mal basis o L2(Ω).
The eigen alues o a domain a e closely ela ed wi h i s geome ic p ope ies, cons i u ing an ac i e a ea o
esea ch called spec al geome y. A classical example o his ela ionship is Weyl’s law, which ela es he
asymp o ics o he eigen alues o he olume o he domain, and a la e esul by McKean and Singe s a es
ha he pe ime e is also de e mined by he eigen alues [MS67]. Mo e esul s o his kind can be ound in
[BG90] and [ dBS88]. O he esul s abou how he geome y o a domain de e mines i s spec um can be
ound in Hen o ’s book [Hen06].
The ques ion o he de e mina ion o a domain gi en he se o i s Laplace eigen alues was posed by
Ma k Kac in his amous pape “Can one hea he shape o a d um?” [Kac66]. Since hen, he answe has
been ound o be nega i e in gene al; in pa icula , o euclidean polygons, he i s example o a pai o
non-isome ic polygons wi h he same spec um is due o Go don, Webb and Wolpe [GWW92]. Howe e ,
he e a e posi i e esul s when we es ic he de e mina ion o a class o domains, he mos success ul
o which was ound by Zeldi ch [Zel09], who p o ed spec al de e mina ion o analy ic domains wi h wo
classes o symme ies.
Less is known abou domains wi h less egula bounda ies, he simples o which a e polygons. In he
case o iangles, i has been p o en ha he whole spec um o he Laplace ope a o de e mines he shape
o a iangle ([Du 88], wi h a ecen simple p oo by [GM13]), and la e Chang and De Tu ck p o ed ha
only a ini e amoun o eigen alues, which depends on λ1and λ2, is enough [CD89]. I is na u al o y o
imp o e he esul o only a ini e and ixed amoun o eigen alues, answe ing he ques ion “Can a human
hea he shape o a iangula d um?”.
Since he space o iangles up o isome ies is has dimension 3, we would expec ha 3 eigen alues
should be enough, bu a p io i i is no clea which ones. An unes and F ei as [AF11] conjec u ed ha
indeed he h ee i s eigen alues λ1, λ2and λ3do de e mine he shape o a iangle. Ins ead, nume ical
e idence by hemsel es seems o indica e ha his is no he case o λ1, λ2and λ4, and in his pape we
will p o e his ac (Theo em 1). This will gi e an example o an obs uc ion o de e mining he shape o a
iangle om a ini e po ion o i s spec um.
3
2 S uc u e o he p oo o Theo em 1
By he scaling o he p oblem, we educe ou sea ch o he se o iangles wi h a ixed base leng h ( oge he
wi h addi ional condi ions ha ensu e ha we only conside one iangle o each simila i y class); ins ead
o looking o all h ee eigen alues λ1, λ2, λ4 o be equal, we jus equi e he quo ien s ξ21 =λ2/λ1and
ξ41 =λ4/λ1 o ake he same alue. Since he eigen alues scale by λ −2when he leng hs o a iangle a e
scaled by , he wo condi ions a e equi alen .
Fixing he i s wo e ices o he iangle o be (0,0) and (1,0), we use he coo dina es (cx, cy) o he
hi d e ex o pa ame ize he sea ch space. Ou app oach consis s in using a opological a gumen o
show ha in each o wo disjoin egions in his pa ame e space he e is a iangle in which ξ21 and ξ41
ake he same p esc ibed alue. Mo e p ecisely, we claim ha he e a e wo dis inc iangles o which
ξ21 =¯
ξ21 := 1.67675 and ξ41 =¯
ξ41 := 2.99372.
Since igo ous calcula ions wi h he compu e a e done using in e al a i hme ic, we need a opological
echnique o ans o m he closed condi ion in o an open condi ion which can be compa ible wi h he in e als.
Fo his we will use he Poinca ´e–Mi anda heo em (see [Mi 40]):
Theo em 2. Gi en wo con inuous unc ions , g : [−1,1]2→Rsuch ha (x, y)has he same sign as x
when x=±1and g(x, y)has he same sign as ywhen y=±1, he e exis s a poin (x, y)∈[−1,1]2such ha
(x, y) = g(x, y)=0.
The egions ha we will conside a e wo pa allelog ams a ound he poin s A= (0.63500,0.27500) and
B= (0.84906,0.31995), designed such ha ξ21 and ξ41 ha e app oxima ely a cons an alue each in a pai
o opposi e edges. Using he compu e we will e i y ha he unc ions ξ21 −¯
ξ21 and ξ41 −¯
ξ41 each ha e a
cons an and opposi e sign in opposi e edges o he pa allelog am, and hence by he heo em, oge he wi h
he well known con inui y o eigen alues, we will conclude ha such wo dis inc iangles exis .
The ec o s de ining he pa allelog am a e o cou se ob ained om he in e se o an app oxima ion o
he di e en ial o he R2- alued uncion (ξ21, ξ41) a he poin s Aand B, eescaled so ha he inc emen o
ξ21 and ξ41 along each ec o is app oxima ely 0.01. Explici ly, hey a e gi en by:
• 21,A = (0.04269082311,0.01148489707)
• 41,A = (−0.02082984105,−0.00255790585)
• 21,B = (0.01717048015,0.01266844996)
• 41,B = (0.03590363291,0.01065148385)
Mo eo e , o he poin Ain he nega i e 21,A di ec ion, he ec o is sh inked by a ac o 0.6, because
a e ha he beha io de ia es om linea and he bound becomes wo se a e his displacemen . Hence
he e ices o he pa allelog ams will be
•A+ 21,A + 41,A,A+ 21,A − 41,A,A−0.6 21,A − 41,A,A−0.6 21,A + 41,A.
•B+ 21,B + 41,B,B+ 21,B − 41,B,B− 21,B − 41,B,B− 21,B + 41,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
Figu e 1: Nume ical app oxima e plo o he quo ien s ξ21 (discon inuous lines) and ξ41 (con inuous lines)
a ound he egion o in e es . The poin s Aand Band he alida ed pa allelog ams a e shown in ed.
5
This se up is displayed in Figu e 1 oge he wi h a plo o non igo ous con ou lines o he eigen alue
quo ien s.
The poin wise e i ica ion o he alues ξ21 and ξ41, which depends on an accu a e calcula ion o λi o
i= 1,2,4, consis s o wo s eps. The i s one, ea ed in Sec ion 3, is abou showing ha he compu ed
eigen alues ac ually co espond o he o de ed ones λ1, λ2and λ4; in o de o do ha , we will p o e a lowe
bound o λ5combining echniques om he Fini e Elemen Me hod wi h igo ous bounds linking he ini e
dimensional p oblem o he in ini e dimensional one. The second s ep consis s o inding accu a e alues o
ou eigen alues ha lie below he h eshold ob ained in he i s pa , which implies ha hey will indeed
ha e o be he ou lowes ones. This is done using he Me hod o Fundamen al Solu ions and ecen igo ous
bounds based on he L2no m o he bounda y e o , explained in Sec ion 4.
We emphasize he di icul y o inding he o de o an eigen alue, which is a global p oblem, compa ed
o he local easie ask o e ining i s alue. To he bes knowledge o he au ho , his is he i s compu e -
assis ed p oo in which hese wo dis inc , local and global me hods a e used o e i y eigen alues o an
ope a o .
In o de o ex end he compu e -assis ed poin wise e i ica ions o he eigen alues λ1, λ2and λ4 o he
con inuous egion in which hey need o be alida ed, we will use explici con inui y esul s o he spec um
o an ellip ic ope a o wi h he domain o ex end he bounds o a neighbo hood o he e i ied poin s. The
de ails o his pa a e explained in Sec ion 5. Sec ion 6 con ains he de ails o he implemen a ion and
execu ion o he p oo . Finally, he Appendix con ains he codes used o e i y an uppe bound o he 4 h
eigen alue a a poin and o alida e a segmen o iangles.
We end his sec ion by in oducing he eade o how a compu e -assis ed p oo wo ks. In he ecen
yea s, he applica ion o calcula ions done by compu e s o ma hema ical p oo s ha e become mo e popula
due o he inc emen o compu a ional esou ces, bu in o de o make su e ha hei esul s a e igo ous,
we need o con ol he e o s ha loa ing poin a i hme ic can accumula e. This is usually done by means o
in e al a i hme ic, in which he da a ha a compu e s o es o a eal numbe is an in e al ( wo endpoin s,
o a midpoin and a adius) o eal numbe s, s o ed by wo loa ing poin numbe s, ins ead o jus one.
Ope a ions be ween in e als a e implemen ed o e u n in e als which a e gua an eed o con ain e e y
possible esul when he ope ands belong o he inpu in e als. Fo example, i [x]=[x, x] and [y]=[y,y]
a e wo in e als, hei sum will can be gi en by he in e al [x] + [y]=[x+y, x +y] and hei p oduc
by [x]·[y] = [min{xy, xy, xy, xy},max{xy, xy, xy, xy}]. The same ule applies o unc ion implemen a ions:
a unc ion e alua ed on [x] should e u n an in e al con aining e e y (x) o x∈[x]. We e e o
he book [Tuc11] o an in oduc ion o alida ed nume ics, in which mos o he echniques used he e a e
explained, and o [GS18] o a mo e speci ic ea men o compu e -assis ed p oo s in PDE.
3 Sepa a ion o he i s ou eigen alues
In o de o ind a igo ous lowe bound o he i h eigen alue o a iangle we will use a ecen bound
ound by Liu [Liu15], which is simila o he one in [CG14] bu simpli ies he hypo heses and imp o es he
cons an . Bo h use he non-con o ming Fini e Elemen Me hod o C ouzeix–Ra ia ; o he igo ous bounds
wi h con o ming ini e elemen s we e explo ed, like [LO13], bu he bound is wo se and he me hod is ha de
6
o implemen wi h alida ed nume ics because i s mass ma ix is no diagonal.
The C ouzeix–Ra ia ini e-elemen me hod uses a iangula ion o he domain Ω, which in ou case we
will ake o be he i ial iangula ion gi en by N2 iangles wi h sides equal o 1/N o he o iginal one and
simila o i . The basis unc ions a e indexed by in e io edges o he iangula ion: i Eis a common edge
o iangles T1,T2, he basis unc ion ψEis he unique unc ion suppo ed on T1∪T2such ha es ic ed o
each iangle is a ine, akes he alue 1 in he midpoin o Eand he alue 0 in he midpoin s o he o he
edges o T1and T2.
We de ine he coe icien s o he s i ness and mass ma ices A= (aEF ),B= (bEF ) by he bilinea o ms
aEF =ZΩ∇ψE·∇ψFbEF =ZΩ
ψEψF
Fo ou choice o iangula ion, Bis simply a mul iple o he iden i y 2|Ω|I/3N2, whe eas Ais a
spa se ma ix. This will allow us o wo k wi h a ma ix eigen alue p oblem ins ead o a gene alized one.
The main esul ha we will use is he ollowing [Liu15, Theo em 2.1 and Rema k 2.2]:
Theo em 3. Conside a polygonal domain Ωwi h a iangula ion so ha each iangle has diame e a
mos h. Le λkbe he k- h eigen alue o Ωand λk,h he k- h eigen alue o he C ouzeix–Ra ia disc e ized
p oblem o Ω. Then λh,k
1 + C2
hλh,k ≤λk,(1)
whe e Ch≤0.1893his a cons an .
In o de o be able o deal wi h app oxima e eigen alues we will need in addi ion he ollowing lemma
om [Pa 80, Theo em 15.9.1].
Lemma 4. Le (˜
λh,˜
uh)be an app oxima ed algeb aic eigenpai such ha ˜
λhis close o some λh han
o any o he disc e e eigen alue. Suppose ha he coe icien ec o ˜
uhis no malised wi h espec o B,
kB˜
uhkB−1=k˜
uhkB= 1. Then he algeb aic esidual := A˜
uh−˜
λhB˜
uhsa is ies
|λh−˜
λh| ≤ k kB−1.
Rema k 5. We can combine Theo em 3 wi h Lemma 4 using he mono onici y o (1), using λh,k −k kB−1
as a lowe bound o λh,k ins ead.
I is easy o ob ain es ima ions ˜
λhwi h a e y small esidual. The ha des pa he e be o e applying
he heo em is o check ha hey ha e indeed he co ec index, i.e., ha hey a e close o he app op ia e
λh han o any o he disc e e eigen alue. This is why we need o con ol he whole spec um o he disc e e
p oblem, and we will do ha by applying Ge shgo in’s disks heo em a e pe o ming some Gi ens o a ions
o he ma ix o ge use ul bounds.
Mo e p ecisely, in o de o ge a lowe bound o λ5we need o sepa a e he i s 5 eigen alues om
he es so ha we can con ol hem, and in o de o do ha we will pe o m Gi ens o a ions un il he
in e als p o ided by Ge shgo in’s heo em can be sepa a ed in wo disjoin componen s, one con aining he
5 smalles eigen alues and he o he con aining he es . I his holds, hen he s ong e sion o Ge shgo in’s
7
heo em will gua an ee an uppe bound o λh,5. The e o e, p o ided ha he esiduals a e all e y small
and ha all app oxima e eigen alues a e di e en (which happens in ou se ing), Lemma 4 will gua an ee
ha he e a e 5 dis inc disc e e eigen alues below he uppe bound and he e o e hey will be o ced o
ha e he co ec indices.
This allows us o e i y λh,5wi h an e o only depending on i s esidual, and using Rema k 5, ge a
lowe bound o λ5. The e o e, he i s ou eigen alues can be sepa a ed jus by checking ha hey a e
dis inc and smalle han his lowe bound.
4 Rigo ous eigen alue bounds o indi idual iangles
Ou app oach o ind igh bounds o he eigen alues o iangles uses he Me hod o Fundamen al Solu ions
(MFS), in oduced by Fox, Hen ici and Mole in [FHM67] and mo e ecen ly e i ed by Be cke and T e e hen
[BT05]. In his me hod, a unc ion uis w i en as a linea combina ion o unc ions φi(1 ≤i≤N) ha
sa is y poin wise he equa ion (∆ + λ)φi= 0 o a ixed λ. The coe icien s a e chosen o op imize he
p oximi y o he unc ion o he eigenspace o he ac ual eigen alue λj, in a sense made p ecise in [BT05],
and his is measu ed by he leas singula alue o a ce ain ma ix ha in ol es he alues o ua disc e e
poin s o he bounda y ∂Ω. This pa ame e is minimized wi h espec o λby using a golden a io sea ch.
This p o ides a candida e λ∈Rand coe icien s ci o which u(x) = PN
i=1 ciφi(x) can be compu ed wi h
a bi a y p ecision.
The unc ions ha we will use o he MFS consis o wo ypes: he i s ones a e o he o m φ(x) =
Y0(√λ|x−x0|), o x0a poin ou side Ω. The second ype o unc ions a e pa ame ized by a e ex o he
iangle and a posi i e in ege j, and ake he o m ψj( , θ) = Jjα(√λ ) sin(jαθ), whe e ( , θ) a e he pola
coo dina es o he poin wi h espec o a e ex in he iangle whose o al angle is π/α, and θis measu ed
om an adjacen side. The i s kind o unc ions allow us o app oxima e he unc ion in he in e io o he
iangle and nea he sides, while he second kind gi es he co ec asymp o ic beha io o he solu ion nea
he e ices o he iangle.
The main ool ha we will use o ind igo ous bounds o eigen alues is he L2bound gi en by Ba ne
and Hassell [BH11]. Howe e , hei me hod is op imized o high eigen alues, so we will ha e o adap some
o he s eps o ou case o small eigen alues. We summa ize he main esul s ha we will use. Le Ω be a
iangle and u∈C2(Ω) be nonze o such ha (∆ + λ)u= 0. Conside he ension
[u] = kukL2(∂Ω)
kukL2(Ω)
.
Le λj, ujbe he sequence o eigen alues and eigen unc ions o Ω, sa is ying (∆ + λj)uj= 0 wi h Di ichle
null bounda y condi ions. Le jbe he no mal de i a i e o uj, de ined on ∂Ω. We de ine he ope a o
A(λ) = X
λj
jh j,·i
(λ−λj)2,
and i s decomposi ion as a sum o h ee:
Anea (λ) = X
|λ−λj|≤√λ
jh j,·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
Figu e 2: Ge shgo in in e als be o e applying he Gi ens o a ions o iangle A.
−20,000 −10,000 0 10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000
Figu e 3: Ge shgo in in e als a e applying he Gi ens o a ions o iangle A. The dashed e ical line
a posi ion 790 sepa a es he i s 5 eigen alues om he es .
Acknowledgmen s
Fi s o all, I would like o hank my supe iso Ja ie G´omez-Se ano o hos ing me a P ince on Uni e si y,
sugges ing me he p oblem, p o iding me e e ences, gi ing me a lo o guidance and sha ing wi h me many
15

−200 0 200 400 600 800 1,000 1,200 1,400 1,600 1,800
Figu e 4: Zoom o he sepa a ion o he i s 5 eigen alues a e he Gi ens o a ions.
(a)
(b)
Figu e 5: G id used o alida e a lowe bound o kukL2(Ω) o iangle B, wi h (a) he i s eigen unc ion
and (b) he second eigen unc ion plo ed on op.
use ul discussions. I am e y g a e ul o he Ma hema ics Depa men o P ince on Uni e si y, in pa icula
he G adua e P og am, o unding my uni e si y ee and allowing my s ay as a isi o esea che . I would
also like o hank CFIS and he adminis a o s o he Mobili y p og am o gi ing me he opo uni y o
do his esea ch in ano he uni e si y and unding me. I also hank he MOBINT schola ship o p o iding
pa ial inancial suppo o my s ay.
16
A Code o he sepa a ion o he smalles eigen alues
This p og am eads he pa ame e Nand he coo dina es o he hi d e ex o he iangle and e u ns a
lowe bound o λ5.
#include "a b.h"
#include "a b_ma .h"
#include "linalg.h"
#include "a bcc.h"
#include <casse >
#include <cma h>
#include <ios eam>
#include <map>
#include < ec o >
// Compile wi h -DPREC=1024
using namespace s d;
using namespace alglib;
in di[6] = {0, 0, 1, -1, 1, -1};
in dj[6] = {1, -1, 0, 0, -1, 1};
in side(in i, in j) {
i (i % 2 == 1) e u n 2;
i (j % 2 == 0) e u n 0;
e u n 1;
}
double sq(double ) { e u n * ;}
/* s uc A Real de ini ion omi ed o b e i y */
in n, m;
oid gi ens( ec o < ec o <A bReal> >& A, double lowba , double upba ,
in nsmall) {
ec o <pai <A bReal, in > > s (m);
o (in i=0;i<m;++i){
s [i] = make_pai (A[i][i], i);
17
}
so (s .begin(), s .end());
ec o <bool> small(m, alse);
o (in i = 0; i < nsmall; ++i) {
small[s [i].second] = ue;
}
bool done = alse;
in i = 0;
while (!done and i < 10) {
++i ;
done = ue;
o (in p=0;p<m;++p){
A bReal ad = 0.;
o (in q=0;q<m;++q){
i (q != p) {
ad = ad + A[q][p].abs();
}
}
A bReal ma gin;
i (small[p]) {
ma gin = A bReal(upba ) - A[p][p];
}else {
ma gin = A[p][p] - A bReal(lowba );
}
i ( ad < ma gin) {
con inue;
}
done = alse;
p in ("i %d: changing column %d wi h ma gin ", i , p);
ma gin.p in ();
o (in q=0;q<m;++q){
i (q != p && ma gin / (m - 1) < A[p][q].abs()) {
A Real = (A Real(A[q][q]) - A Real(A[p][p])) /
(A Real(A[p][q]) * A Real(2));
A Real = ( ).sign() / ( .abs() + ( .sq() + 1).sq ());
A bReal _a b( .x);
A bReal c = 1_a / ( _a b.sq() + 1_a).sq ();
A bReal s = c * _a b;
o (in i=0;i<m;++i){
A bReal u = A[p][i];
18
A bReal = A[q][i];
A[p][i] = c * u - s * ;
A[q][i] = c * + s * u;
}
o (in i=0;i<m;++i){
A bReal u = A[i][p];
A bReal = A[i][q];
A[i][p] = c * u - s * ;
A[i][q] = c * + s * u;
}
}
}
}
}
}
ec o <pai <A bReal, A bReal> > ge shgo in( ec o < ec o <A bReal> >& A) {
ec o <pai <A bReal, A bReal> > e (m);
o (in i=0;i<m;++i){
e [i]. i s = A[i][i];
e [i].second = 0.;
o (in j=0;j<m;++j){
i (j != i) {
e [i].second = e [i].second + A[j][i].abs();
}
}
}
so ( e .begin(), e .end());
e u n e ;
}
in main() {
ios::sync_wi h_s dio( alse);
cou .se (ios:: ixed);
cou .p ecision(10);
cin >> n;
double x_s, y_s;
cin >> x_s >> y_s;
A bReal x(x_s), y(y_s);
19
A bReal a ea = y / 2_a;
A bReal l_a = (y * y + (x - 1_a) * (x - 1_a)).sq ();
A bReal l_b = (y * y + x * x).sq ();
A bReal l_c = 1;
A bReal len[3];
len[0] = l_a;
len[1] = l_b;
len[2] = l_c;
A bReal h_a = y / l_a;
A bReal h_b = y / l_b;
A bReal h_c = y;
A bReal heigh [3];
heigh [0] = h_a;
heigh [1] = h_b;
heigh [2] = h_c;
A bReal cos_a = x / l_b;
A bReal cos_b = (1_a - x) / l_a;
A bReal cos_c = (x * (x - 1_a) + y * y) / (l_a * l_b);
A bReal cosine[3];
cosine[0] = cos_a;
cosine[1] = cos_b;
cosine[2] = cos_c;
map<pai <in ,in >, in > encode;
o (in i=0;i<2*n-2;++i){
o (in j=0;i+j<2*n-2;++j){
i (i%2==0o j%2==0){
encode[make_pai (i, j)] = m++;
}
}
}
ec o < ec o <A bReal> > ma (m, ec o <A bReal>(m));
A bReal mass_diag = (a ea * 2) / (3 * n * n);
o (in i=0;i<2*n-2;++i){
o (in j=0;i+j<2*n-2;++j){
20

i (i%2==1and j%2==1)con inue;
in u = encode[make_pai (i, j)];
in side_u = side(i, j);
ma [u][u] = (a ea * 8) / heigh [side_u].sq();
o (in d=0;d<6;++d){
in i2 = i + di[d];
in j2 = j + dj[d];
i (i2 >= 0 and j2 >= 0 and i2+j2<2*n-2and
(i2 % 2 == 0 o j2 % 2 == 0)) {
in = encode[make_pai (i2, j2)];
in side_ = side(i2, j2);
A bReal mp = -(cosine[3 - side_u - side_ ] * a ea * 4) /
(heigh [side_u] * heigh [side_ ]);
ma [u][ ] = mp;
ma [ ][u] = mp;
}
}
}
}
o (in i=0;i<m;++i){
o (in j=0;j<m;++j){
ma [i][j] = ma [i][j] / mass_diag;
}
}
eal_2d_a ay ma _d;
ma _d.se leng h(m, m);
o (in i=0;i<m;++i){
o (in j=0;j<m;++j){
ma _d(i, j) = ma [i][j].ge _app ox_double();
}
}
ae_in _ num_eigs_ ound;
eal_1d_a ay eigs_d;
eal_2d_a ay eigen ec o s;
21
asse (sma ixe d (ma _d, m, 1, 1, 0.0, 1000.0, num_eigs_ ound, eigs_d,
eigen ec o s));
ce << "Found " << num_eigs_ ound << " below 1000 n";
ec o <A bReal> disc e e_eigs_lb(num_eigs_ ound);
o (in i = 0; i < num_eigs_ ound; ++i) {
cou << eigs_d(i) << endl;
A bReal lambda(eigs_d(i));
ec o <A bReal> ec(m);
o (in j=0;j<m;++j){
ec[j] = A bReal(eigen ec o s(j, i));
}
A bReal num = 0., den = 0.;
o (in j=0;j<m;++j){
A bReal esidue = -lambda * ec[j];
o (in k=0;k<m;++k){
esidue = esidue + A bReal(ma [j][k]) * ec[k];
}
num = num + esidue.sq();
den = den + ec[j].sq();
}
A bReal e o = (num / den).sq ();
p in ("E o : ");
e o .p in ();
disc e e_eigs_lb[i] = lambda - e o ;
}
asse (num_eigs_ ound >= 6);
double ba ie = (eigs_d[4] + eigs_d[5]) / 2;
// Successi ely apply Gi ens o a ions o educe he ma gin.
// I does no wo k i one ies o sepa a e hem wi h he
// ba ie om he beginning.
gi ens(ma , -100000, 1000000, 10);
gi ens(ma , -80000, 1000000, 10);
gi ens(ma , -60000, 1000000, 10);
gi ens(ma , -50000, 1000000, 10);
gi ens(ma , -20000, 1000000, 10);
gi ens(ma , -10000, 1000000, 10);
gi ens(ma , -1000, 100000, 10);
22
gi ens(ma , 0, 100000, 10);
gi ens(ma , 1000, 10000, 10);
gi ens(ma , 1000, 2000, 10);
gi ens(ma , ba ie , 2000, 5);
gi ens(ma , ba ie , ba ie , 5);
p in ("Disc e e lambda_5 success ully sepa a ed n");
ec o <pai <A bReal, A bReal> > ge = ge shgo in(ma );
A bReal liu_cons an = 0.1893;
A bReal lambda5lb =
disc e e_eigs_lb[4] /
(1_a + disc e e_eigs_lb[4] * liu_cons an .sq() / A bReal(n * n));
p in ("Ve i ied lowe bound o lambda_5: ");
lambda5lb.p in ();
}
B Code o alida ing an in e al o iangles
This p og am eads he name o he hi d e ex (Ao B), he indexs o he eigen alue quo ien o alida e
(“21” o “41”), a le e s anding o he sign (’p’ o plus, ’m’ o minus) o ξ−¯
ξ ha we ha e o alida e, he
cu en coo dina e c, and he o al numbe o coo dina es ncin which his side is di ided (c∈ {1, . . . ,nc}).
This will alida e a segmen a ound he posi ion gi en by 2(c−0.5)/nc−1 in he co esponding side o he
pa allelog am (whe e he posi ion is no malized so ha e ices a e a posi ions −1 and 1), and he adius
`o his segmen will be he alue e u ned by he p og am.
Fo example, i he inpu is “B 21 p 312 1000” and he p og am e u ns a alue o 0.0013, his means
ha o all he poin s in he segmen be ween B+ 21,B −0.1898 41,B and B+ 21,B −0.1872 41,B he
co esponding ξ21 has been e i ied o be g ea e han ¯
ξ21.
The p og am uses he A bTaylo class, which is no included he e o b e i y, bu compu es simul a-
neously he Taylo polynomial and i s esidue, and uses hem o ob ain an enclosu e o he unc ion in an
in e al, as desc ibed in [Tuc11].
#include "a b.h"
#include "a b_hypgeom.h"
#include "a bcc.h"
#include "a bse ies.h"
#include "a b aylo .h"
23
#include <algo i hm>
#include <casse >
#include < s eam>
#include <ios eam>
#include <map>
#include < ec o >
// Compile wi h -DPREC=128
using namespace s d;
// Decla ed he e and implemen ed in a sepa a e ile.
oid ill_in_ ec o s( ec o <double>& _coe s, double& _lambda,
ec o < ec o <double>>& _sou ces, double cx0, double cy0,
double lmin, double lmax);
namespace {
// Coo dina es o he hi d e ex. Se as global because hey a e used
// e e ywhe e.
A bReal cx, cy;
/* GENERAL FUNCTIONS */
// Fundamen al solu ion a ound an ex e nal cha ge poin ‘(xs, ys)‘.
A bTaylo und_sol_cha ge(cons A bTaylo & x, cons A bTaylo & y,
cons A bReal& xs, cons A bReal& ys,
cons A bReal& lambda, in o de ) {
A bTaylo a g = (x - A bTaylo ::cons an (xs, o de )).sq() +
(y - A bTaylo ::cons an (ys, o de )).sq();
a g = a g.sq ();
a g *= lambda.sq ();
e u n a g.bessely0();
}
// Angle be ween wo ec o s.
A bReal ccw_angle(cons A bReal& x1, cons A bReal& y1, cons A bReal& x2,
cons A bReal& y2) {
A bReal = (x1.sq() + y1.sq()).sq ();
A bReal c = x1 / ;
24
ans = A bReal(aux);
a _clea (aux);
}else {
A bReal al1 =
pa h_d s(coe , lambda, sou ces, 2 * num, 2 * den, e _goal, xx, yy);
A bReal al2 =
pa h_d s(coe , lambda, sou ces, 2 * num + 1, 2 * den, e _goal, xx, yy);
a b_min(ans, al1, al2, PREC);
}
mag_clea (e _ al);
e u n ans;
}
// Minimum o u^2 in he pa h gi en by (xx[i], yy[i]).
A bReal pa h_minsq(cons ec o <A bReal>& coe , cons A bReal& lambda,
cons ec o < ec o <A bReal>>& sou ces, double goal_d,
cons ec o <A bReal>& xx, cons ec o <A bReal>& yy) {
mag_ e _goal;
mag_ini (e _goal);
mag_se _d(e _goal, goal_d);
in den = 18;
A bReal lb;
a b_pos_in (lb);
o (in num = 0; num < den; ++num) {
A bReal al =
pa h_d s(coe , lambda, sou ces, num, den, e _goal, xx, yy).sq();
a b_min(lb, lb, al, PREC);
}
mag_clea (e _goal);
e u n lb;
}
A bReal compu e_l2meanlb(cons ec o <A bReal>& coe , A bReal lambda,
cons ec o < ec o <A bReal>>& sou ces) {
double goal = 5e-2;
in n = 8; // Bounding L2 no m using 8*8 iangles,
A bReal l2sum = 0_a;
A bReal = 0.8; // Using 0.64 o he o al a ea o a oid he bo de s.
31

o (in coun e = 0; coun e < n * n ; ++coun e ) {
au o pa h = ge _ iangle(coun e , n , );
A bReal pmsq =
pa h_minsq(coe , lambda, sou ces, goal, pa h. i s , pa h.second);
l2sum += pmsq * .sq() / A bReal(n ).sq();
}
A bReal l2meanlb = l2sum.sq ();
p in ("Lowe bound o he L2 mean: ");
l2meanlb.p in ();
e u n l2meanlb;
}
/* COMPUTE THE EIGENVALUE WITH AN ABSOLUTE ERROR */
// Re u ns he pai {lambda, abs_e }, wi h lambda in [lbound, ubound]. ‘index‘
// speci ies he posi ion o he eigen alue in he spec um (1,2,...).
pai <A bReal, A bReal> compu e_all(in index, double lbound, double ubound) {
ec o <double> coe _d;
ec o < ec o <double>> sou ces_d;
double lambda_d;
ill_in_ ec o s(coe _d, lambda_d, sou ces_d, cx.ge _app ox_double(),
cy.ge _app ox_double(), lbound, ubound);
in nc = coe _d.size();
in ns = sou ces_d.size();
ec o <A bReal> coe (nc);
ec o < ec o <A bReal>> sou ces(ns, ec o <A bReal>(2));
A bReal lambda = lambda_d;
o (in i = 0; i < nc; ++i) {
coe [i] = A bReal(coe _d[i]);
}
o (in i = 0; i < ns; ++i) {
sou ces[i][0] = A bReal(sou ces_d[i][0]);
sou ces[i][1] = A bReal(sou ces_d[i][1]);
}
p in ("Candida e lambda: ");
lambda.p in ();
32
A bReal meanl2lb = compu e_l2meanlb(coe , lambda, sou ces);
double goal_d = 1e-5;
a _ goal;
a _ini (goal);
a _se _d(goal, goal_d);
A bReal l2bd y = side_l2_ ec(coe , lambda, sou ces, goal);
p in ("Valida ed a goal %e n", a _ge _d(goal, 10));
p in ("L2 no m a he bounda y: ");
l2bd y.p in ();
A bReal a ea = cy / 2;
A bReal o all2lb = meanl2lb * a ea.sq ();
A bReal ension = l2bd y / o all2lb;
A bReal side_a = (cy.sq() + (1_a - cx).sq()).sq ();
A bReal side_b = (cy.sq() + cx.sq()).sq ();
A bReal semip = (side_a + side_b + 1_a) / 2_a;
A bReal in ad =
((semip - side_a) * (semip - side_b) * (semip - 1_a) / semip).sq ();
A bReal c_omega = 4_a * (1_a + in ad) / in ad;
A bReal ub_all_lambda = 800_a; // Rough uppe bound o all lambda (1 o 4).
A bReal a_ a = 7_a * c_omega;
A bReal a_ ail = 7_a * c_omega / lambda.sq ();
A bReal aux = 1_a / ension.sq();
aux -= a_ a + a_ ail;
A bReal abs_e = 2_a * A bReal(ubound) / (aux * in ad);
asse (0_a < abs_e );
abs_e = abs_e .sq ();
p in (
"Rigo ous alue o lambda ob ained. Follows lambda, i s ension and i s "
"absolu e e o . n");
lambda.p in ();
ension.p in ();
abs_e .p in ();
a _clea (goal);
e u n {lambda, abs_e };
33
}
}// namespace
in main() {
cha p name, signc;
s ing indexs;
in coo d, ncoo d;
cin >> p name >> indexs >> signc >> coo d >> ncoo d;
A bReal cx0, cy0, 21x, 21y, 41x, 41y, posi ion;
// Heu is ics o ind he app op i a e eigen alue.
// This will be alida ed la e by a sepa a e p og am.
ec o <double> lmin, lmax;
i (p name == ’A’) {
// Coo dina es:
cx0 = 0.635;
cy0 = 0.275;
// Vec o s o he pa allelog am:
21x = 4.269082311683548 / 100;
21y = 1.148489707350921 / 100;
41x = -2.082984105473002 / 100;
41y = -0.255790585210996 / 100;
i (indexs == "21") {
posi ion = (A bReal(coo d) - 0.5_a) / A bReal(ncoo d) * 2_a - 1_a;
lmin = {180, 320};
lmax = {280, 450};
}else {
// In his case, he pa allelog am is sh inked om below:
// he posi ion akes alues in [-0.6, 1].
posi ion = (A bReal(coo d) - 0.5_a) / A bReal(ncoo d) * 1.6_a - 0.6_a;
lmin = {180, 620};
lmax = {280, 727. - (2 * coo d - ncoo d) / ncoo d * 27};
}
}else {
cx0 = 0.849057346949971;
cy0 = 0.319950941965592;
21x = 1.717048015465781 / 100;
21y = 1.266844996461298 / 100;
34
41x = 3.590363291549854 / 100;
41y = 1.065148385561495 / 100;
posi ion = (A bReal(coo d) - 0.5_a) / A bReal(ncoo d) * 2_a - 1_a;
i (indexs == "21") {
i (signc == ’m’) {
lmin = {160, 280};
lmax = {255, 380};
}else {
lmin = {150, 280};
lmax = {190, 380};
}
}else {
i (signc == ’m’) {
double cen e _l1 = 214 - coo d * 23.0 / ncoo d;
double cen e _l4 = 637 - coo d * 65.0 / ncoo d;
lmin = {cen e _l1 - 20, cen e _l4 - 20};
lmax = {cen e _l1 + 20, cen e _l4 + 20};
}else {
double cen e _l1 = 197 - coo d * 20.0 / ncoo d;
double cen e _l4 = 593 - coo d * 60.0 / ncoo d;
lmin = {cen e _l1 - 20, cen e _l4 - 20};
lmax = {cen e _l1 + 20, cen e _l4 + 20};
}
}
}
i (indexs == "21") {
i (signc == ’m’) {
cx = cx0 - 21x + posi ion * 41x;
cy = cy0 - 21y + posi ion * 41y;
}else {
cx = cx0 + 21x + posi ion * 41x;
cy = cy0 + 21y + posi ion * 41y;
}
}else {
i (signc == ’m’) {
cx = cx0 + posi ion * 21x - 41x;
cy = cy0 + posi ion * 21y - 41y;
}else {
cx = cx0 + posi ion * 21x + 41x;
35
cy = cy0 + posi ion * 21y + 41y;
}
}
s ing ilename = "ou pu _";
ilename += p name;
ilename += indexs;
ilename += signc;
ilename += o_s ing(coo d);
ilename += "o ";
ilename += o_s ing(ncoo d);
ilename += ". x ";
eopen( ilename.c_s (), "w", s dou );
au o p1 = compu e_all(1, lmin[0], lmax[0]);
au o pk = compu e_all((indexs == "21" ? 2 : 4), lmin[1], lmax[1]);
A bReal xi = pk. i s / p1. i s ;
A bReal e xi = (pk.second + p1.second * xi) / (p1. i s - p1.second);
p in ("Xi and e o : n");
xi.p in ();
e xi.p in ();
A bReal goalxi = (indexs == "21" ? 1.67675 : 2.99372);
A bReal ma gin = (signc == ’p’ ? xi - goalxi - e xi : goalxi - xi - e xi);
A bReal eps = ma gin / (pk. i s * (1_a + xi + ma gin));
A bReal del a = ((2_a * A bReal::pi()).sq() * eps) /
(1_a + (2_a * A bReal::pi()).sq() * eps);
A bReal cc, kk; // Auxillia y a iables o calcula e ell.
i (indexs == "21") {
cc = cy * ( 41x.abs() + 2_a * 41y.abs());
kk = 41y.abs();
}else {
cc = cy * ( 21x.abs() + 2_a * 21y.abs());
kk = 21y.abs();
}
// Leng h alida ed a his poin h ough p opaga ion o he e o .
// Gi en as he coe icien o he ec o 21 o 41 ha one can mo e.
A bReal ell = (2_a * del a * cy * kk + cc -
(4_a * del a * cy * kk * cc + cc.sq()).sq ()) /
36

(2_a * del a * kk.sq());
i (p name == ’A’ && indexs == "41") {
// Reno malize ell so i ep esen s a leng h in [-1, 1]:
ell *= 1.25;
}
p in ("Value o ell: ");
ell.p in ();
}
C Code o op imizing he eigen alue
This code implemen s he unc ion ill_in_ ec o s decla ed in he ile abo e, which inds an app oxima e
eigen alue and eigen ec o in a non- igo ous ashion. I does so by op imizing he smalles nonze o singula
alue o a ma ix ha imposes he bounda y condi ions, as explained abo e.
#include <s dio.h>
#include <sys/ ime.h>
#include <algo i hm>
#include <u ili y>
#include "linalg.h"
#include "s da x.h"
#include "op imiza ion.h"
#include "special unc ions.h"
#include <boos /ma h/special_ unc ions/bessel.hpp>
#include <boos /ma h/ ools/minima.hpp>
using namespace alglib;
using namespace s d;
namespace {
double cx, cy;
inline double sq(double x) { e u n x*x;}
double no m(cons eal_1d_a ay& a) {
double s = 0;
o (in i = 0; i < a.leng h(); ++i) {
s += sq(a[i]);
}
37
e u n sq (s);
}
eal_1d_a ay ope a o *(cons eal_1d_a ay& a, double x) {
eal_1d_a ay c;
in n = a.leng h();
c.se leng h(n);
o (in i=0;i<n;++i){
c[i] = a[i] * x;
}
e u n c;
}
eal_1d_a ay ope a o /(cons eal_1d_a ay& a, double x) {
eal_1d_a ay c;
in n = a.leng h();
c.se leng h(n);
o (in i=0;i<n;++i){
c[i] = a[i] / x;
}
e u n c;
}
eal_1d_a ay ope a o +(cons eal_1d_a ay& a, cons eal_1d_a ay& b) {
eal_1d_a ay c;
in n = a.leng h();
c.se leng h(n);
o (in i=0;i<n;++i){
c[i] = a[i] + b[i];
}
e u n c;
}
eal_1d_a ay ope a o -(cons eal_1d_a ay& a, cons eal_1d_a ay& b) {
eal_1d_a ay c;
in n = a.leng h();
c.se leng h(n);
o (in i=0;i<n;++i){
c[i] = a[i] - b[i];
}
38
e u n c;
}
double und_sol_cha ge_d(double x, double y, double x0, double y0,
double lambda) {
double = sq (sq(x - x0) + sq(y - y0));
e u n bessely0( * sq (lambda));
}
double und_sol_ e ex_d(double x, double y, in , in j, double lambda) {
double x e [3] = {0, 1, cx};
double y e [3] = {0, 0, cy};
in nx = ( + 1) % 3;
in p = ( + 2) % 3;
double x = x e [ ];
double y = y e [ ];
double nx _[2] = {x e [nx ] - x , y e [nx ] - y };
double p _[2] = {x e [p ] - x , y e [p ] - y };
double cu _[2] = {x - x , y - y };
eal_1d_a ay nx , p , cu ;
nx .se con en (2, nx _);
p .se con en (2, p _);
cu .se con en (2, cu _);
double hmx = a an2( p [1], p [0]) - a an2( nx [1], nx [0]);
i ( hmx < 0) {
hmx += 2 * M_PI;
}
i ( hmx > 2 * M_PI) {
hmx -= 2 * M_PI;
}
double h = a an2( cu [1], cu [0]) - a an2( nx [1], nx [0]);
i ( h < 0) {
h += 2 * M_PI;
}
i ( h > 2 * M_PI) {
h -= 2 * M_PI;
}
double = no m( cu );
double alpha = M_PI * (j + 1) / hmx;
double sl = sq (lambda);
39
e u n boos ::ma h::cyl_bessel_j(alpha, sl * ) * sin(alpha * h);
}
// Gene a es sou ce poin s: ‘ns‘ pe side and an adi ional ‘nex a‘
// si ua ed a ound he op e ex. ‘nex a‘ should be odd.
eal_2d_a ay gen_sou ces(in ns, in nex a) {
double be a = 0.01;
eal_2d_a ay sou ces;
sou ces.se leng h(3 * ns + nex a, 2);
double x e [3] = {0, 1, cx};
double y e [3] = {0, 0, cy};
o (in =0; <3;++ ){
eal_1d_a ay cu , nx ;
double cu _[2] = {x e [ ], y e [ ]};
cu .se con en (2, cu _);
in nx = ( + 1) % 3;
double nx _[2] = {x e [nx ], y e [nx ]};
nx .se con en (2, nx _);
o (in i = 0; i < ns; ++i) {
double a = 0.5 * (1.0 - cos(double(i + 1) * M_PI / double(ns)));
eal_1d_a ay p = cu * (1.0 - a ) + nx * a ;
eal_1d_a ay ec = nx - cu ;
ec = ec / no m( ec);
double pe p_[2] = { ec[1], - ec[0]};
eal_1d_a ay pe p;
pe p.se con en (2, pe p_);
o (in c=0;c<2;++c){
sou ces( * ns + i, c) = p [c] + be a * pe p[c];
}
}
}
o (in i = 0; i < nex a; ++i) {
sou ces(3 * ns + i, 0) = cx + (double(i) - (nex a - 1) / 2) * 0.02;
sou ces(3 * ns + i, 1) = cy + 0.1;
}
e u n sou ces;
}
// Gene a es ‘ni‘ in e io poin s a andom.
40
[Pa 80] B. N. Pa le . The Symme ic Eigen alue P oblem. Classics in Applied Ma hema ics. Socie y o
Indus ial and Applied Ma hema ics, 1980.
[Rel40] F. Rellich. Da s ellung de Eigenwe e on ∆u+λu = 0 du ch ein Randin eg al. Ma h. Z.,
46:635–636, 1940.
[Tuc11] W. Tucke . Valida ed nume ics: a sho in oduc ion o igo ous compu a ions. P ince on Uni-
e si y P ess, 2011.
[ dBS88] M. an den Be g and S. S isa kuna ajah. Hea equa ion o a egion in R2wi h a polygonal
bounda y. J. London Ma h. Soc. (2), 37(1):119–127, 1988.
[Zel09] S. Zeldi ch. In e se spec al p oblem o analy ic domains. II. Z2-symme ic domains. Ann. o
Ma h. (2), 170(1):205–269, 2009.
47