---
jupyter:
  description: 'Tuning a marimba or vibraphone bar with planex2d and NGSolve: the free bar''s bending modes against beam
    theory (overtone at 2.76 times the fundamental, nodes at 0.224 L), and the depth of the arch under its middle that
    puts the first overtone two octaves up, at four times the fundamental.'
---

# Tuning a marimba bar

Why do marimba and vibraphone bars have an arch cut into their underside? A plain bar sounds out of tune with itself,
because the first overtone of a free bar rings at 2.756 times its fundamental, not at a musical interval. Thinning the
middle of the bar lowers the fundamental much more than the overtone, and an arch of the right depth puts the overtone
two octaves up, at four times the fundamental. A vibraphone bar is tuned to these two. Marimba makers also place the
next overtone, near ten times the fundamental, with a more shaped cut. The bar hangs on cords through the nodes of its
fundamental, where it does not move. On a plain bar, they lie 0.224 of the length from each end.

The example is an aluminum bar like a vibraphone's, 350 mm long and 13 mm thick, in plane stress. NGSolve and SciPy find
its bending modes, which are checked against beam theory for the plain bar, and SciPy's `brentq` finds the depth of the
arch that tunes the overtone to four times the fundamental.

This page is part of a collection of [physics examples](physics.md) using planex2d.

![A tuned bar swinging in its fundamental and, four times as fast, in its first overtone](../_static/gallery/marimba-bar-tuning.gif){.result}

```python
import math

import matplotlib.pyplot as plt
import ngsolve as ngs
import numpy as np
import planex as px
import scipy.sparse as sp
from scipy.optimize import brentq
from scipy.sparse.linalg import eigsh
```

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

## The bar

The bar is 350 long and 13 thick, in millimeters. It is aluminum, with $E$ = 70000 N/mm², Poisson's ratio 0.33 and the
density 2.7e-9 t/mm³, so that the frequencies come out in Hz. Its width does not enter, since a bar bending in the plane
of its length and thickness has the same frequencies however wide it is. The arch is a disc taken out of the underside.
Its circle passes through the bottom edge 0.3 of the length on either side of the middle and reaches `depth` up in the
middle:

```python
L, t = 350, 13
E, nu, rho = 70e3, 0.33, 2.7e-9


def bar(depth, chord=0.6 * L):
    """the bar with an arch of the given depth cut under its middle"""
    face = px.QuickSketch().Rectangle(L, t).Face()
    if depth > 0:
        R = (chord**2 / 4 + depth**2) / (2 * depth)
        face = face - px.Circle((L / 2, depth - R), R).Face()
    return face
```

## Its modes

The bar hangs free, without supports. Its modes are the eigenvectors of its stiffness and its mass, which SciPy's
`eigsh` finds nearest to a shift below zero. The first three are the bar moving as a rigid body, at zero frequency, and
are skipped. The shift is of the size of the eigenvalues, the squares of the circular frequencies. With a shift too
close to zero, the rigid-body modes would leave the shifted matrix nearly singular and the frequencies inaccurate:

```python
mu, lam = E / (2 * (1 + nu)), E * nu / (1 - nu**2)


def modes(face, k=8):
    """the bending frequencies in Hz and the modes of the free bar, lowest first"""
    mesh = ngs.Mesh(px.OCCGeometry(face).GenerateMesh(maxh=t / 3))
    V = ngs.VectorH1(mesh, order=3)
    u, v = V.TnT()

    def stress(w):
        strain = ngs.Sym(ngs.Grad(w))
        return 2 * mu * strain + lam * ngs.Trace(strain) * ngs.Id(2)

    a = ngs.BilinearForm(ngs.InnerProduct(stress(u), ngs.Sym(ngs.Grad(v))) * ngs.dx).Assemble()
    m = ngs.BilinearForm(rho * u * v * ngs.dx).Assemble()
    A = sp.csr_matrix(a.mat.CSR(), shape=(V.ndof, V.ndof))
    M = sp.csr_matrix(m.mat.CSR(), shape=(V.ndof, V.ndof))
    w2, vectors = eigsh(A, k=k, M=M, sigma=-1e5, which="LM")
    order = np.argsort(w2)[3:]
    return np.sqrt(w2[order]) / (2 * math.pi), vectors[:, order], mesh, V
```

## The plain bar

Beam theory gives a free bar the fundamental $f_1 = \dfrac{4.730^2}{2 \pi L^2} \sqrt{\dfrac{E I}{\rho A}}$ and its first
overtones at 2.7565 and 5.4039 times it. The model lies close to it. Its overtones lie a little lower, because the bar
also shears, which beam theory leaves out and which matters more for the shorter waves of the higher modes:

```python
f, vectors, mesh, V = modes(bar(0))
f_beam = 4.7300407**2 / (2 * math.pi * L**2) * math.sqrt(E * t**2 / 12 / rho)
print("f1", f[0], "| beam", f_beam)
print("f2 / f1", f[1] / f[0], "| beam 2.7565 | f3 / f1", f[2] / f[0], "| beam 5.4039")
```

The nodes of the fundamental are where its vertical motion along the middle of the bar changes sign:

```python
mode = ngs.GridFunction(V)
mode.vec.FV().NumPy()[:] = vectors[:, 0]
xs = np.linspace(0, L, 3501)
up = np.array([mode(mesh(x, t / 2))[1] for x in xs])
print("nodes", xs[:-1][np.sign(up[:-1]) != np.sign(up[1:])], "| 0.2242 L", 0.2242 * L, 0.7758 * L)
```

## Tuning the overtone

A deeper arch lowers both frequencies, the fundamental more, so the ratio of the two grows with the depth.
`brentq` finds the depth at which it is four (Figure 1):

```python
def ratio(depth):
    f = modes(bar(depth))[0]
    return f[1] / f[0]


depths = np.arange(0, 9.5, 0.5)
ratios = [ratio(d) for d in depths]
depth = brentq(lambda d: ratio(d) - 4, 2, 9, xtol=1e-4)
f, vectors, mesh, V = modes(bar(depth))
print("depth", depth, "| f", f[:3], "| f2 / f1", f[1] / f[0], "| f3 / f1", f[2] / f[0])
# plt.plot(depths, ratios); plt.show()
# px.Plot(bar(depth))
```

```python nbsphinx="hide-input"
Chart("Depth of the arch", "First overtone / fundamental", dots={"NGSolve": (depths, ratios)},
      marks={"Tuned, two octaves": (depth, 4)}, lines={"Two octaves": ([0, 9], [4, 4])},
      title="The overtone ratio against the depth of the arch")
```

The arch takes 7.7 of the 13 in the middle and lowers the fundamental to about two fifths of the plain bar's. The
next overtone lands at about 9.3 times the fundamental, short of the ten that marimba makers aim for with a more
shaped cut. Figure 2 shows the tuned bar and its first two modes.

```python nbsphinx="hide-input"
def dense(points, step=1.0):
    """a closed polyline with a point at least every step along each edge"""
    ring = np.vstack([points, points[:1]])
    out = []
    for p, q in zip(ring[:-1], ring[1:]):
        k = max(1, int(np.ceil(np.hypot(*(q - p)) / step)))
        out += [p + (q - p) * i / k for i in range(k)]
    return np.array(out)


tuned = bar(depth)
outline = np.array(tuned.wires[0].Polylines(0.05)[0])
if np.allclose(outline[0], outline[-1]):
    outline = outline[:-1]
outline = dense(outline)
along = np.linspace(0, L, 701)
shapes = []
for k in (0, 1):
    mode = ngs.GridFunction(V)
    mode.vec.FV().NumPy()[:] = vectors[:, k]
    lift = np.array([mode(mesh(x, t - 1))[1] for x in along])
    lift *= 0.08 * L / np.abs(lift).max()
    moved = outline + np.column_stack([np.zeros(len(outline)), np.interp(outline[:, 0], along, lift)])
    shapes.append(px.FromPolyline(moved, closed=True).Face())
Row((tuned, "The tuned bar"), (shapes[0], f"Fundamental, {f[0]:.0f} Hz"), (shapes[1], f"First overtone, {f[1]:.0f} Hz"),
    caption="The tuned bar and its first two modes, the motion magnified")
```

See [netgen and NGSolve](../integrations/netgen/index.md) for meshing a planex2d face with netgen.

```python nbsphinx="hidden"
import sys
from matplotlib.animation import FuncAnimation
sys.path.insert(0, "..")
from docs_figures import FACE, INK, PLANEX
from gallery import save_animation

# the tuned bar in its fundamental (top) and its first overtone (bottom), which swings four times as fast
lifts = []
for k in (0, 1):
    mode = ngs.GridFunction(V)
    mode.vec.FV().NumPy()[:] = vectors[:, k]
    lift = np.array([mode(mesh(x, t - 1))[1] for x in along])
    lifts.append(np.interp(outline[:, 0], along, lift * 0.08 * L / np.abs(lift).max()))
fig, ax = plt.subplots(figsize=(6, 4))


def frame(i):
    ax.clear()
    phase = 2 * math.pi * i / 48
    for lift, speed, y in ((lifts[0], 1, 45), (lifts[1], 4, -45)):
        tuned.Translated((0, y)).Plot(ax=ax, face_col=None, line_col=FACE, legend=False)
        moved = outline + np.column_stack([np.zeros(len(outline)), math.cos(speed * phase) * lift + y])
        px.FromPolyline(moved, closed=True).Face().Plot(ax=ax, face_col=PLANEX, line_col=INK, legend=False)
    ax.set_aspect("equal", adjustable="box")
    ax.set_xlim(-25, L + 25)
    ax.set_ylim(-133, 133)
    ax.axis("off")


animation = FuncAnimation(fig, frame, frames=48, interval=60)
save_animation(animation, "marimba-bar-tuning", fps=16, bare=True)
plt.close(fig)
```