---
jupyter:
  description: 'The two drums of Gordon, Webb and Wolpert as planex2d faces: their eigenvalues
    from NGSolve agree, so the shape cannot be heard. Area and perimeter can.'
---

# Isospectral drums: can one hear the shape of a drum?

Mark Kac asked in 1966 whether the tones of a drum, the eigenfrequencies of its membrane, determine its shape. In 1992
Gordon, Webb and Wolpert answered no: two different drums, each made of seven half-squares, have exactly the same
eigenvalues. The example is their pair of drums as planex2d faces: netgen meshes them, NGSolve computes their
spectra, and the spectra agree. What can be heard is the area and the perimeter, which the two drums share.

![The first modes of two different drums that sound the same](../_static/gallery/isospectral-drums.png){.result}

```python
import math
import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation
import ngsolve as ngs
from ngsolve import solvers
import planex as px
```

```python nbsphinx="hidden"
from docs_figures import Row
```

## The drums

The two drums are eight-sided polygons, and each covers seven right isosceles triangles with legs of 2 (Figure 1).

```python
drums = {
    "first": [(-3, -3), (-3, -1), (1, 3), (1, 1), (3, 1), (1, -1), (-1, -1), (-1, -3)],
    "second": [(-3, 1), (1, 1), (1, 3), (3, 1), (1, -1), (-1, -1), (-1, -3), (-3, -1)],
}
faces = {name: px.FromPolyline(points, closed=True).Face() for name, points in drums.items()}
```

```python nbsphinx="hide-input"
Row((faces["first"], "The first drum"), (faces["second"], "The second drum"),
    caption="The drums of Gordon, Webb and Wolpert")
```

**Check:** both have the area of seven halves of a 2 × 2 square, 14, and the same perimeter:

```python
for name, face in faces.items():
    print(name, "area", face.mass, "perimeter", face.outer.mass)
```

## The tones

The membrane's eigenvalues satisfy −Δu = λu inside and u = 0 on the rim. NGSolve's preconditioned inverse iteration
(PINVIT) finds the lowest twelve, and with the matrix's own inverse as the preconditioner, it is plain inverse
iteration on a block of vectors. The computation runs on a coarse and a finer mesh. On the finer one, the largest
relative difference between the two drums drops by almost an order of magnitude, while each eigenvalue itself
settles. The remaining difference is discretization error, because the two drums are meshed differently, so a single
eigenvalue can also move the other way.

```python
def spectrum(face, maxh, order, count=12):
    mesh = ngs.Mesh(px.OCCGeometry(face).GenerateMesh(maxh=maxh))
    fes = ngs.H1(mesh, order=order, dirichlet=".*")
    u, v = fes.TnT()
    a = ngs.BilinearForm(ngs.grad(u) * ngs.grad(v) * ngs.dx).Assemble()
    m = ngs.BilinearForm(u * v * ngs.dx).Assemble()
    pre = a.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky")
    values, vectors = solvers.PINVIT(a.mat, m.mat, pre, num=count, maxit=30, printrates=False)
    modes = []
    for vector in vectors:
        mode = ngs.GridFunction(fes)
        mode.vec.data = vector
        modes.append(mode)
    return list(values), mesh, modes
```

The coarse mesh has elements of 0.5 and order 3, the finer one elements of 0.25 and order 4:

```python
results = {}
for maxh, order in ((0.5, 3), (0.25, 4)):
    results[maxh] = {name: spectrum(face, maxh, order) for name, face in faces.items()}
    first, second = results[maxh]["first"][0], results[maxh]["second"][0]
    print("maxh", maxh, "order", order)
    for n, (x, y) in enumerate(zip(first, second)):
        print(" ", n + 1, round(x, 8), round(y, 8), "relative difference", f"{abs(x - y) / x:.1e}")
```

## Same tone, different shape

The figure shows the first four modes of each drum. The eigenvalues match, but the shapes of the vibration do
not.

```python nbsphinx="hide-input"
def picture(ax, mesh, mode):
    xy = [v.point for v in mesh.vertices]
    triangles = [[v.nr for v in el.vertices] for el in mesh.Elements(ngs.VOL)]
    values = mode.vec.FV().NumPy()[:mesh.nv]
    scale = max(abs(values.min()), abs(values.max()))
    grid = Triangulation([p[0] for p in xy], [p[1] for p in xy], triangles)
    ax.tripcolor(grid, values / scale, shading="gouraud", cmap="RdBu_r", vmin=-1, vmax=1)
    ax.set_aspect("equal")
    ax.axis("off")


fig, axes = plt.subplots(2, 4, figsize=(14, 7))
for row, name in zip(axes, faces):
    values, mesh, modes = results[0.25][name]
    for n, ax in enumerate(row):
        picture(ax, mesh, modes[n])
        ax.set_title(f"{name}, λ{n + 1} = {values[n]:.6f}", fontsize=9)
plt.show()
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from gallery import save_scene, save_figure, save_animation, HIDDEN
card, pair = plt.subplots(1, 2, figsize=(10, 5))
for ax, name in zip(pair, faces):
    values, mesh, modes = results[0.25][name]
    picture(ax, mesh, modes[0])
    faces[name].outer.Plot(ax=ax, legend=False)
save_figure(card, "isospectral-drums", bare=True)
```