The Quest to Finding Chladni Patterns, Part 1: Theory and Algorithms

Background and Mathematical Formulation

Take a solid metal plate (or any other rigid material) and sprinkle some particles on it. It can either be sand, salt, or even baking soda. Then, generate some vibration on the plate by using either a violin bow to excite the plate, or place it on a wave generator and set it to specific frequencies. What you will observe are the beautiful patterns generated as shown in the figure above. This experiment sounds relatively simple, but as simple as it sounds, this experiment reveals some extremely complex and intensive underlying mathematical principles. In fact, the problem of understanding these patterns challenged some of the brightest mathematical minds of the time, and the research to understand these patterns spans for hundreds of years, and still continues till the modern day.

How does this work? Well, it turns out these are examples of standing waves in the two-dimensional space. In one dimension, a “node” in a standing wave is a point along the wave where the amplitude is at minimum. Similarly, in two dimensions, there are “nodal lines,” which are static points on the plate where the particles settle.

In the spring semester of 2024, our physics project required selecting an experiment to demonstrate a scientific law. Chladni's law was a compelling choice, and one good idea was to prove the constant wave speed across a plate using the two-dimensional wave equation model

\begin{equation} u_{tt} = c^2\nabla^2u + Q(\mathbf{x}, t),\quad\quad \mathbf{x}\in\Omega, \end{equation}

where $c$ here denotes the wave speed of the plate, $u_{tt}$ is the plate's displacement, and $Q$ is the forcing function (the wave generator).

By assuming that the plate is supported by tension, one can derive the two-dimensional wave equation in the Cartersian coordinate to be

In the next grand section, we will discuss the simulation process of this model on the square and circular domain, and to discuss the flaws of this model (lots of flaws, but no spoiler yet). Basically, the emphasized, initial goal of this project is to simulate Chladni patterns, prove that the wave speed is constant, and validate Chladni's law (see the section Chladni's Law below). For a much more detailed explanation, please refer to my project report.

Read more: Historical Notes

The study of wave and vibration theory dates its way back to the time of Pythagoras in the fifth century BC, where the Ionian School of Natural Philosophy first introduced the scientific methods to study the real world and natural phenomena. But just as with many other scientific discoveries before Newton, we lacked the tools to rigorously understand the underlying theory behind them. Fast forward to the Renaissance period. Galileo also noticed that pieces of bristle on the sounding board only move in certain areas. Leonardo da Vinci also observed the strange patterns generated by particles in response to vibration.

The first rigorous investigation of these patterns so called “Chladni patterns” is of course, attributed to German physicist and musician Ernst Chladni (1756 - 1827) (not Euler or Gauss this time). The idea of sprinkling sand on a plate and bow it is not modest, but Chladni took our understanding to a new level by coming up with a formula that describe the pattern, and hypothesized a formula that predicts the number of rings on the circular plate. This formula is later proven true by Kirchoff and Rayleigh.

Chladni's Law

Part of the formula states that the vibrating frequency is proportional to four times the number of rings on the plate squared, that is, $f \sim 4n^2$, where $n$ is the number of rings (see in Fig. 4). The full formula that he derived is \begin{equation} f \sim (m+2n)^2, \end{equation} where $m$ here is the number of diameter lines. This is called Chladni's law. We don't see any $m$'s in the figure below, there is a reason for this, and we'll get to that part later.

The general model and the boundary condition

Depends on whether you clamp your edge or not, clamp the edge partially, or simply support the plate by placing some other objects below the plate; the boundary condition may be different. For this project, we will emphasize the free boundary condition. That is, only the wave generator supports the plate.

The free boundary condition means that there are no forces acting on the edges of the plate; that is, expressed mathematically,

which gave us four boundary conditions. The uniqueness of the solution requires our equation to have two more initial conditions. Similar to ordinary differential equations, these conditions are the initial displacement and velocity of the plate at $t=0$. That is,

The Square Plate

Consider the wave equation on a rectangle $\Omega = \{(x,y)\in\R^2\;|\; -a\leq x,y\leq a\}$,

Here, $a$ is half the length of the plate, and also note that we denote the length of the plate to be $L$. A standard way to solve a forced differential equation is to solve its homogenous case first, i.e. when $Q=0$.There are many ways to solve this initial-boundary value problem, and one common, and easy-to-understand technique is via separation of variables. If we assume that $u$ is separable and let $u(x,y,t) = X(x)Y(y)T(t)$, the equation becomes

Divide both sides by $XYT$, we obtain

Notice that $\frac{X''}{X}+\frac{Y''}{Y}=-\lambda$ is constant if and only if $X''/X=-\mu$ and $Y''/Y = -\nu$, where $\mu$ and $\nu$ are themselves constants and $\lambda = \mu + \nu$. The two constants above yield two eigenvalue problems with Neumann boundary conditions:

and the time equation yields the third eigenvalue problem:

The detail solution of these eigenvalue problems is beyond the scope of this article. For the detail solution, again, please refer to my project report. The general idea is we need to solve the eigenvalue problems (\ref{eq:X1}) and (\ref{eq:Y1}) first, and from there we solve the third eigenvalue problem, since $\lambda = \mu + \nu$. After carrying out all the hard work, the eigenvalues for equations (\ref{eq:X1}) and (\ref{eq:Y1}) are

and the solution to the homogeneous equation is the double Fourier series

where,

The coefficients $a_n$ and $\overline{a_n}$ are defined as $a_n = \frac{1 + (-1)^n}{2},\overline{a_n} = \frac{1 - (-1)^n}{2}$, and,

The forcing term $Q$ doesn't change the overall shape of the patterns, but it does change the displacement of the plate overtime, hence the time function $T_{nm}(t)$ becomes

Here, $q_{nm}$ is a function defined on the disk of radius $r$ — the radius of the tip of the wave generator,

The Circular Plate

Consider the polar wave equation on a disk $\Omega = \{(r,\theta)\in\R^2\;|\; 0 \leq r\leq a, 0 \leq \theta\leq 2\pi\},$

Let $u(r, \theta, t) = R(r)\Theta(\theta)T(t)$. Skipping all the hard work., we may arrive at the final solution to be

Algorithm for the Square Plate

The reason for all those derivations above is to plot our theoretical patterns. Recall that our final solution is derived as

Before we design our algorithm, there are a few things to note:

• In realistic experimentation, there is no way to tell what the initial conditions $f(x,y)$ and $g(x,y)$ are. However, we can assume that the initial velocity, $u(\mathbf{x},0) = g(x,y)$ is zero. Moreover, the forcing term is non-zero only when $n=m$, and it approaches zero as $\omega\to\infty$. This means that it does not matter that much for the overall shape of the plot, but it only affects the insignificant microscopic height of the plate. The solution is now simplified to
  \begin{equation} u(\mathbf{x}, t) = \sum_{n\geq 0}\sum_{m\geq 0}X_n(x)Y_m(y)A_{nm}\cos(\omega t). \end{equation}
• For the initial condition $f(x,y)$, we can pick a function $f(x,y)$ such that $|A_{ij}|\neq |A_{nm}|$, for any $(i, j)\neq (n,m)$. The reason why we can do this is that the realistic vertical displacement of the plate is extremely small, so the choice of $f$ doesn't matter that much. We simply choose the function $f$ that satisfies $|A_{ij}|\neq |A_{nm}|$ for the clearest patterns.
• What we are really interested in is the shape of the nodal lines. Therefore, for the sake of simplicity, we can choose certain time points such that $\cos(\omega t) = 1$, that is, $t = 0,\; 2\pi/\omega,\; 4\pi/\omega,\ldots$. In this way, if we set $u(\mathbf{x}, t) = 0$, the coefficients $A_{nm}$ will be cancelled out, leaving us with the contour equation
  \begin{equation} u(\mathbf{x}, t) = \sum_{n\geq 0}\sum_{m\geq 0}X_n(x)Y_m(y) = 0. \end{equation}

The final step is to notice the angular frequency $\omega$. Our solution $u(\mathbf{x}, t)$ is derived as the superposition sum of all other possible solutions. However, this is a double infinite series; how are we going to plot this function? This is a lesson I learned throughout my time as an engineering student: Occasionally, we must reinterpret abstract mathematics to make it meaningful in the context of the physical world..

• Recall that we derived our angular frequency to be
  \begin{equation} \omega = \frac{\pi c}{L}\sqrt{n^2 + m^2}. \end{equation}
  What we are really interested in is not the general super-position solution, but solutions at certain frequencies $\omega$. After all, each frequencies must only correspond to some certain shape of the plate. This means that each $\omega$ must correspond so some combination of integer pairs $n$ and $m$.
• For example, if $\omega = \frac{\pi c}{L}\sqrt{5}$, then all possible combinations of $n$ and $m$ must be $(1, 2)$ and $(2, 1)$ (since $1^2 + 2^2 = 5$). Notice previously that we use the coefficients $|A_{nm}|$ in absolute value, so negative coefficients are possible. The two equations that describe the Chladni pattern found at this frequency are
  \begin{equation} X_1(x)Y_2(y) \pm X_2(x)Y_1(y)=0, \end{equation}
  For equations with more sum of squares combinations, there will be more combinations for a single mode. For example, for $\omega = \frac{\pi c}{L}\sqrt{65}$, all the possible implicit functions are
  \begin{equation} X_1(x)Y_8(y) \pm X_4(x)Y_7(y) \pm X_7(x)Y_4(y) \pm X_8(x)Y_1(y) =0\ \end{equation}
  We need not consider the first term due to symmetry. In general, if there are $N$ terms in the equation, there will be a total of $2^{N-1}$ possible plots.
• Special Case: If $n^2 + m^2$ is a perfect square ($1, 4, 9, 16, 25,\ldots$), then we also need to consider the pairs with 0. For example, if $n^2 + m^2 = 25$, then besides the standard pairs of $(3, 4)$ and $(4, 3)$, we also need to consider the pairs $(0, 5)$ and $(5, 0)$. This is because since we don't know what $f(x,y)$ is, we also don't know whether $A_{00}$ is 0. If it's zero, then we just plot the regular pairs $(3, 4)$ and $(4, 3)$,
  \begin{equation} X_3(x)Y_4(y) + X_4(x)Y_3(y)=0,\quad\text{and}\quad X_3(x)Y_4(y) - X_4(x)Y_3(y)=0. \end{equation}
  If $A_{00}$ is not zero, meaning that the summation actually starts at index $(n,m)=(0,0)$, then we also need to consider eight other patterns
  \begin{equation} X_0(x)Y_5(y) \pm X_3(x)Y_4(y) \pm X_4(x)Y_3(y) \pm X_5(x)Y_0(y) =0. \end{equation}
  In total, for the case of $S = 25$, there will be a total of $2+8=10$ plots.

In summary, the algorithm to plot the patterns of the square plate is just an algorithm to search for pairs of sum of squares $(n, m)$.

    # ALGORITHM 1: PLOT SQUARE CHLADNI PATTERNS
    
    # generate different sign combinations for X using an array of 1's and -1's.
    # always keep the first element positive (1)
    function generate_sign_combinations(pairs) 
        
    # define the basis function
    # we don't use Y(y), or m here because it's basically the exact same function as X
    function X(x, n)
        return mod(n + 1, 2) * cos(n * pi * x/ 2.0) + mod(n, 2) * sin(n * pi * x/ 2.0)

    # plot the full wave function
    # input: x, y, pair (n, m) and sign array, [1, -1, -1] for instance
    function u(x, y, pairs, signs)
        return u

    # search for sum-of-squares pairs
    # find all pairs of (n, m) such that n^2 + m^2 = number
    # return value: a list of tuples (n, m) that satisfies the requirement
    function find_pairs(number) 

    # Check if a given number can be represented as sums of squares pairs
    # return true if it is a sum of square
    function exist_sos(number) 

    # search for sum of squares
    # plot the sum of squares
    function plot_sum_of_squares(start_index, end_index)
        # compute total_plots as the difference between end_index and start_index
        # initialize empty lists: equation, sums_of_squares, and description (to store the title of each plot)
        # iterate over numbers and check for sum of squares:
        #   for each integer k, check if it can be written as a sum of squares by exist_sos(k)
        #   if true, find all pairs of k using find_pairs(k)
        #   for each pair (n, m) generate all possible sign combinations
        #   for each sign combination:
        #       if current plot is in between start_index and end_index, compute u(x, y, pairs, signs)
        #       store u in the list equation, and k in sums_of_squares
        #   stop when the total_plots is reached.
        # for i in range of total_plots:
        #   assign colors for each unique pairs (n. m).
        #   plot the contour plot at z = 0
    

Algorithm for the Circular Plate

The circular plate is a bit easier, since we don't need to design an algorithm to search for sum of squares and perfect squares exceptions. Recall that for a circular plate,

By applying a similar technique as the square plate, we can eliminate $c_{nm}, d_{nm}$, and the forcing term. The circular plate is easier to implement because we don't have to search for any sum of squares, so each Chladni pattern is simply a pair of arbitrary $n$ and $m$. The equations that we need to plot are simply

    # ALGORITHM 2: PLOT CIRCULAR CHLADNI PATTERNS
    
    # Function to get the m-th zero of the derivative of the Bessel function of order n
    function bessel_derivative_zero(n, m):
        # redefine the first zero of J_0 prime to be 0
        # default is not zero for some reason
        if n == 0 then
            zeros = list(jnp_zeros(n, m + 1))
            zeros.insert(0, 0)
            return zeros[m]
        else
            # Get the m-th zero, treating first zero as the 0th zero
            return jnp_zeros(n, m + 1)[-1]  

    # Initial function f(r, theta)
    function f(r, theta):
        return r * np.cos(theta)

    # Two coefficients
    # use built-in integration tools to obtain these
    # input: pair (n, m), radius a
    function a_nm(n, m, a)   
    function b_nm(n, m, a)

    # Define the function J_n(z_nm * r / a) * (a_nm * cos(n*theta) + b_nm * sin(n*theta))
    def u(n, m, r, theta, a)
    

The next part of this article will be dedicated to implementing the code in python and compare the theoretical results to the experimental result. By then, we may know whether the wave speed is constant or not. You can read Part 2 by pressing the directory below.