Can you hear the shape of a drum?

5. Can you hear the shape of a drum?#

https://www2.math.upenn.edu/~kazdan/425S11/Drum-Gordon-Webb.pdf

You can’t. This paper shows two different regions leading to the same spectrum. We numerically verify the result.

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

First, we compute eigenvalues highly accurate, with high order finite elements and geometric refinement towards corners of the geometry.

shape1 = MoveTo(0,0).Line(1,1).Line(0,1).Line(-2,0).Line(0,1).Line(-1,-1).Line(1,-1).Line(1,0).Close().Face()
shape2 = MoveTo(0,0).Line(1,0).Line(0,1).Line(-2,2).Line(0,-1).Line(-1,0).Line(1,-1).Line(1,0).Close().Face()

Draw (shape1)
Draw (shape2);

def Calc(mesh):
    fes = H1(mesh, order=8, dirichlet=".*")
    u,v = fes.TnT()
    A = BilinearForm(grad(u)*grad(v)*dx).Assemble().mat
    M = BilinearForm(u*v*dx).Assemble().mat
    pre = A.Inverse(freedofs=fes.FreeDofs(), inverse="sparsecholesky")
    num = 10
    evals, evecs = solvers.LOBPCG(A, M, pre=pre, num=num, maxit=20, printrates=False)  
    gfu = GridFunction(fes, multidim=num)
    for i in range(num):
        gfu.vecs[i] = evecs[i]
    return evals, gfu

shape1.vertices.hpref=1
mesh1 = Mesh(OCCGeometry(shape1, dim=2).GenerateMesh(maxh=0.2))
mesh1.RefineHP(4)
Draw (mesh1)
evals1, evecs1 = Calc(mesh1)

shape2.vertices.hpref=1
mesh2 = Mesh(OCCGeometry(shape2, dim=2).GenerateMesh(maxh=0.2))
mesh2.RefineHP(4)
Draw (mesh2)
evals2, evecs2 = Calc(mesh2)

print ("eig1 = ", list(evals1))
print ("eig2 = ", list(evals2))

eig1 =  [10.151776047619883, 14.62203890731103, 20.702237531002538, 26.15022979608839, 28.992311539006707, 36.83718003643364, 42.38794282218638, 46.16558160042807, 49.34802200544566, 52.214616369441465]
eig2 =  [10.151776046760101, 14.622038908791549, 20.70223752336885, 26.150229797296547, 28.99231154171411, 36.83718003351098, 42.38794282853997, 46.165581599904115, 49.34802200544565, 52.21461635691456]

Draw (evecs1, animate=True, min=-1, max=1);
Draw (evecs2, animate=True, min=-1, max=1);

5.1. Computations with structured meshes#

We generate meshes manually, such that symmetries used in the proof of the iso-spectral property can be carried over to the discrete setting.

import netgen.meshing as ngm
mesh1s = ngm.Mesh(dim=2)

pnts = [ (0,0), (1,1), (1,2), (0,2), (-1,2), (-1,3), (-2,2), (-1,1), (0,1) ]
els = [(0,1,8), (8,1,3), (1,2,3), (7,8,3), (7,3,4), (6,7,4), (6,4,5)]
pis = [ mesh1s.Add (ngm.MeshPoint(ngm.Pnt(*p,0))) for p in pnts ]
idx_dom = mesh1s.AddRegion("vol", dim=2)
for el in els:
    mesh1s.Add(ngm.Element2D(idx_dom, [pis[i] for i in el]))
idx_bnd = mesh1s.AddRegion("bnd", dim=1)
for i in range(9):
    mesh1s.Add(ngm.Element1D([pis[i], pis[(i+1)%9]], index=idx_bnd))
for l in range(3):
    mesh1s.Refine()
mesh1s = Mesh(mesh1s)

print (mesh1s.GetMaterials(), mesh1s.GetBoundaries())
Draw (mesh1s);

('vol',) ('bnd',)

mesh2s = ngm.Mesh(dim=2)

pnts = [ (0,0), (1,0), (1,1), (0,2), (-1,3), (-1,2), (-2,2), (-1,1), (0,1) ]
els = [(0,1,2), (0,2,8), (8,2,3), (7,8,3), (7,3,5), (7,5,6), (5,3,4) ]

pis = [ mesh2s.Add (ngm.MeshPoint(ngm.Pnt(*p,0))) for p in pnts ]
idx_dom = mesh2s.AddRegion("vol", dim=2)
for el in els:
    mesh2s.Add(ngm.Element2D(idx_dom, [pis[i] for i in el]))
idx_bnd = mesh2s.AddRegion("bnd", dim=1)
for i in range(9):
    mesh2s.Add(ngm.Element1D([pis[i], pis[(i+1)%9]], index=idx_bnd))
for l in range(3):
     mesh2s.Refine()
mesh2s = Mesh(mesh2s)

Draw(mesh2s);

evals1s, evecs1s = Calc(mesh1s)
print (list(evals1s))

evals2s, evecs2s = Calc(mesh2s)
print (list(evals2s))

[10.152622114299056, 14.62287510462176, 20.703893282458427, 26.15056481024013, 28.99387721709089, 36.837905300538594, 42.388902631551446, 46.16584469874167, 49.348022005445515, 52.21710076899968]

[10.152622114299048, 14.622875104621796, 20.703893282458395, 26.150564810240333, 28.993877217090866, 36.83790530053864, 42.38890263155148, 46.165844698741694, 49.348022005445614, 52.21710076899978]

Although the eigenvalues are not highly accurate, the eigenvalues of both discrete problems match at floating-point accuracy.

Contents