\[ \newcommand{\lbrac}{\left(} \newcommand{\rbrac}{\right)} \]
Topics: Mathematics/Physics/Standing Waves/Python

The Quest to Finding Chladni Patterns, Part 2: Codes and Results

“That there are such limits to our knowledge of nature, must be borne with patience by every sound mind whether he be a scientist or a workman.”
- Gustav Robert Kirchhoff

Review of Part 1

The previous part of this article talked about one of the theoretical models behind Chladni patterns: the wave equation model. We used this model to construct our contour equation and established two algorithms to plot the square model and the circular model. The square wave equation algorithm requires us to check if a number can be represented as a sum of squares. If it can be represented as sum(s) of squares, then we must search for all possible combinations of $n$ and $m$ such that

\begin{equation} n^2 + m^2 = \mathrm{number} \end{equation}

including pairs with zero. The circular wave equation model is simpler to implement, as we only need to utilize built-in integration methods, and do a double for loop to search over 0 to $n$ and 0 to $m$. Right below is the draft of the two algorithms that we have designed previously.

    # 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)
    # return a list of list
    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, for example, [1, -1, -1]
    function u(x, y, pairs, signs)

    # 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)

In the next part, we will generate these patterns in Python, but before proceeding, remember to install Python (select install for all users) and download related dependencies by using the pip package manager:

        -m pip install matplotlib
        -m pip install numpy
        -m pip install itertools
        -m pip install scipy

In this article, we will use matplotlib to plot our coordinate systems, numpy for the mathematical functions and constants, itertools for convenient sign alternation tool, scipy for integration and built-in special functions.

Generate Patterns in Python

Code for the Square Plate

Most of the code sketch that from the section above and the previous part is already pretty straight forward, but I want to make a few remarks about some of the functions used in this code snippet.

def generate_sign_combinations(pairs):
    sign_combinations = list(itertools.product([1, -1], repeat=(len(pairs) - 1)))
    return [[1] + list(comb) for comb in sign_combinations]  # Keep the first sign positive        
        

The function .product() works similar to a cartesian product, or equivalent to a nested for-loop in the context of computer science. For example, product('ABCD', repeat=2) will output

  AA AB AC AD BA BB BC BD CA CB CC CD DA DB DC DD 

For our case, .product([1, -1], repeat = (len(pairs) - 1)) will generate all possible combinations of 1 and -1 for a length of len(pairs) - 1. We need to use len(pairs) - 1 because we want to keep the first sign positive, as seen in the return statement [[1] + list(comb)]. If pairs = [(0, 5), (3, 4), (4, 3), (5, 0)], then len(pairs) = 4. We want to keep the first term to be positive, always, so .product([1, -1], repeat = (len(pairs) - 1)) will generate

    [1, 1, 1] 
    [1, 1, -1]
    [1, -1, 1] 
    [1, -1, -1]    

The statement sign_combinations = list(itertools.product([1, -1], repeat=(len(pairs) - 1))) assigns a list of [1, -1] into sign_combinations up to len(pairs) - 1. The return statement return [[1] + list(comb) for comb in sign_combinations] returns a list of list, consisting of all combinations of 1 and -1. For the example above, if len(pairs) is 4, then the function will return

[ [1, 1, 1], [1, 1, -1], [1, 1, -1], [1, -1, -1] ]

The next function, $u(x, y, t)$, defined by

def u(x, y, pairs, signs):
    equation = np.zeros_like(x)  # Return an array of 0's with the same type
    for (n, m), sign in zip(pairs, signs):
        equation += sign * X(x, n) * X(y, m)
    return equation        
        

takes in $x$, $y$, the pairs which we need to plots, and the signs list from the generate_sign_combinations function. The output is obviously the z-value of the function $u$. For example, if the pairs of $n$ and $m$ is given as (n, m) = [[0, 5], [3, 4], [4, 3], [5, 0]], and the sign combinations is [1, -1, -1], then the contour equation $u(x, y, t)$ must be

\begin{equation} u(x, y, t_0) = X(x, 0)X(y, 5) - X(x, 3)X(y, 4) - X(x, 4)X(y, 3) - X(x, 5)X(y, 0) = 0 \end{equation}

The next function, exist_sos(number) is very straight-forward. It simply return true if a number can be represented as a sum of squares. Similarly, if a number is found, then the function find_pairs(number) find all such pairs, and then return all of them in a list.

def find_pairs(number):
    return [(n, m) for n in range(0, int(np.sqrt(number)) + 1)
            for m in range(0, int(np.sqrt(number)) + 1) if n**2 + m**2 == number]
    return equation        
        

For example, if number == 145, then the pairs are [(1, 12), (8, 9), (9, 8), (12, 1)]. Most sum of squares numbers will have two or three pairs. Some have four like 145 or 25, but some number can have a quite large number of pairs. For example, if number == 725, then the pairs are [(7, 26), (10, 25), (14, 23), (23, 14), (25, 10), (26, 7)], or if number == 9425, then there are up to 12 pairs, [(4, 97), (20, 95), (31, 92), (41, 88), (55, 80), (64, 73), (73, 64), (80, 55), (88, 41), (92, 31), (95, 20), (97, 4)]. This function is $O(n)$, which is pretty efficient. I think there is a way to make it even more efficient, but it lies beyond the scope of this article.

The final, and longest function is the plot_sum_of_squares(start_index, end_index). This function input the starting and ending pattern number, iterate through all patterns in between, and output them to the plots. We first focus on the first half of the function's code block and start by creating three lists: equations, description, and sum_of_squares.

def plot_sum_of_squares(start_index, end_index):
    # ... something up here        
    
    equations = []
    descriptions = []
    sums_of_squares = []
        

The list equations stores all equations to plot, while the other two list are simply for documentation purposes. The list sum_of_squares is to store the sum-of-squares numbers, call it k, to output to the plot, while description store all pairs that made up that number and the signs combination to output to the terminal.

The next part of the code is the while loop to iterate over all k to search for sum of squares. I want the user to be able to determine which index they want to start and end with to avoid overcrowding the screen with too much plot. Therefore, the len(equations) < total_plots is the breaking condition.

    # ...something up here
    k = 0
    found_plots = 0
    while len(equations) < total_plots:
        if exist_sos(k):
            # find all pairs from k
            pairs = find_pairs(k)                                   

            # find non zero pairs using slicing operator, copying element from second index of list, [1] 
            # to second to last element index of list, [-1]
            non_zero_pairs = pairs[1:-1]                            

            # generate a list of different sign combinations from pairs
            sign_combinations = generate_sign_combinations(pairs)  
            
            # iterates over each combination of signs produced in the previous step. 
            # for each combination, it will attempt to create an equation.
            for signs in sign_combinations:
                found_plots += 1

                # if the current found_plots is between start_index and end_index,
                # then append elements to lists
                # not all generated plots will be kept—only those within this specific range.
                if start_index <= found_plots <= end_index:
                    equations.append((pairs, signs, u(x, y, pairs, signs)))
                    descriptions.append(f"{pairs} {signs}")
                    sums_of_squares.append(k)
                if found_plots >= end_index:
                    break

            # check if non_zero_pairs is not empty
            if non_zero_pairs:
                sign_combinations_0 = generate_sign_combinations(non_zero_pairs)
                for signs in sign_combinations_0:
                    found_plots += 1
                    if start_index <= found_plots <= end_index:
                        equations.append((non_zero_pairs, signs, u(x, y, non_zero_pairs, signs)))
                        descriptions.append(f"{non_zero_pairs} {signs}")
                        sums_of_squares.append(k)
                    if found_plots >= end_index:
                        break
            if found_plots >= end_index:
                break
        k += 1
        

The later half of the function is only to plot all graphs in the list. One technique for visibility of the plot is to assign colors to plots of the same modes, but different signs (different sign combinations of [(1, 12), (8, 9), (9, 8), (12, 1)] will have the same color, for example). Here, pair_color_mapping is a dictionary that maps each unique pair of values (pairs_tuple) to a color.

    colors = ['darkblue', 'red', 'darkgreen', 'deeppink', 'darkviolet', 'blue', 'darkgoldenrod',
             'teal', 'darkred', 'darkcyan']
    pair_color_mapping = {}
    color_index = 0
        

And, there it is, our code for generating Chladni patterns on the square plate using the wave equations model and eigenvalues is completed. The full code is available right below. I added a few options for you to adjust based on your own preference.

import matplotlib.pyplot as plt
import numpy as np
import itertools

# Settings here
title = True
font_title = 9.2        # Title font size
color_plots = True      # All black or color plots
thickness = 1.2         # Plot line thickness
start = 1               # Start at pattern #
end = 72                # End at pattern #

# Range and step size
delta = 0.005
xrange = np.arange(-1.0, 1.0, delta)
yrange = np.arange(-1.0, 1.0, delta)
x, y = np.meshgrid(xrange, yrange)

# Define the basis function
def X(x, n):
    return np.mod(n + 1, 2) * np.cos(n * np.pi * x / 2.0) + np.mod(n, 2) * np.sin(n * np.pi * x / 2.0)

# Sign combination generator
def generate_sign_combinations(pairs):
    sign_combinations = list(itertools.product([1, -1], repeat=(len(pairs) - 1)))
    return [[1] + list(comb) for comb in sign_combinations]  # Keep the first sign positive

# The wave function with contour cut z = 0
def u(x, y, pairs, signs):
    equation = np.zeros_like(x)  # Return an array of 0's with the same type
    for (n, m), sign in zip(pairs, signs):
        equation += sign * X(x, n) * X(y, m)
    return equation

# Check if a given number has sums of squares pairs
def exist_sos(number):
    return any(n**2 + m**2 == number for n in range(0, int(np.sqrt(number)) + 1)
                for m in range(0, int(np.sqrt(number)) + 1))

# Searching combinations of sum of squares
def find_pairs(number):
    return [(n, m) for n in range(0, int(np.sqrt(number)) + 1)
            for m in range(0, int(np.sqrt(number)) + 1) if n**2 + m**2 == number]
    
plt.rcParams['mathtext.fontset'] = 'stix'
plt.rcParams['font.family'] = 'serif'

# Plot the sum of squares patterns from start_index to end_index (1-based indexing)
def plot_sum_of_squares(start_index, end_index):
    total_plots = end_index - start_index + 1
    
    # Error messages
    if end_index <= start_index:
        print("Error: End_index must be greater than start_index.\n")
        return
    if total_plots > 100:
        print("Error: Cannot plot more than 100 plots at a time.\n")
        return

    equations = []
    descriptions = []
    sums_of_squares = []
    
    # Iterating over all possible combinations
    k = 0
    found_plots = 0
    while len(equations) < total_plots:
        if exist_sos(k):
            pairs = find_pairs(k)
            non_zero_pairs = pairs[1:-1]
            sign_combinations = generate_sign_combinations(pairs)
            for signs in sign_combinations:
                found_plots += 1
                if start_index <= found_plots <= end_index:
                    equations.append((pairs, signs, u(x, y, pairs, signs)))
                    descriptions.append(f"{pairs} {signs}")
                    sums_of_squares.append(k)
                if found_plots >= end_index:
                    break
            if non_zero_pairs:
                sign_combinations_0 = generate_sign_combinations(non_zero_pairs)
                for signs in sign_combinations_0:
                    found_plots += 1
                    if start_index <= found_plots <= end_index:
                        equations.append((non_zero_pairs, signs, u(x, y, non_zero_pairs, signs)))
                        descriptions.append(f"{non_zero_pairs} {signs}")
                        sums_of_squares.append(k)
                    if found_plots >= end_index:
                        break
            if found_plots >= end_index:
                break
        k += 1

    ncols = 9  # Fixed number of columns
    nrows = (total_plots + ncols - 1) // ncols  # Calculate the number of rows (// = Floor division)
    
    # Define fig and axis
    fig, axs = plt.subplots(nrows, ncols, figsize=(1.2 * ncols, 4 * nrows))
    
    # Colors
    colors = ['darkblue', 'red', 'darkgreen', 'deeppink', 'darkviolet', 'blue', 'darkgoldenrod', 'teal', 'darkred', 'darkcyan']
    pair_color_mapping = {}
    color_index = 0
    
    # Plotting
    for i in range(total_plots):
        # Print all pairs and signs to the console
        print(f"N = {start + i} | S = {sums_of_squares[i]}: {descriptions[i]}")
        
        # pairs = equations[i][0], signs = equations[i][1], equation = equations[i][2]
        pairs, signs, equation = equations[i]
        pairs_tuple = tuple(pairs)
        if pairs_tuple not in pair_color_mapping:
            pair_color_mapping[pairs_tuple] = colors[color_index % len(colors)]
            color_index += 1
        
        mode_color = pair_color_mapping[pairs_tuple] if color_plots else 'black'
        row, col = divmod(i, ncols)
        plot_title = rf"$N = {start + i},\; S = {sums_of_squares[i]}$"
        
        ax = axs[row, col] if nrows > 1 else axs[col]
        ax.contour(x, y, equation, levels=[0], colors=mode_color, linewidths=thickness)
        
        if title:
            ax.set_title(plot_title, fontsize=font_title, ha='center', rotation='vertical', x=-0.1, y=0)
        ax.grid(True)
        ax.set_aspect('equal')
        ax.set_xticks([])
        ax.set_yticks([])
    
    plt.tight_layout()
    fig.subplots_adjust(hspace=0.2, wspace=0.1)
    plt.show()

# Plotting the sum of squares patterns for all perfect squares
plot_sum_of_squares(start, end)
        

Code for the Circular Plate

As mentioned earlier in this article, the circular plate is somewhat easier to implement, as it doesn't require searching for any sum-of-squares numbers. Each pair of $(n, m)$ is a pattern by itself, so we only need to make a double for loop to iterate over every single of of them. The first thing to do is to of course set up a function to find the m-th zero of the n-th derivative of the bessel function of the first kind, $J_n(x)$. For this, we must import jn, jnp_zeros from scipy.special

# Function to get the m-th zero of the derivative of the Bessel function of order n
def bessel_derivative_zero(n, m):
    if n == 0: #Redefine the first zero of J_0 prime to be 0
        zeros = list(jnp_zeros(n, m + 1))
        zeros.insert(0, 0)
        return zeros[m]
    else:
        return jnp_zeros(n, m + 1)[-1]  # Get the m-th zero, treating first zero as the 0th zero
        

One of the weird things about the jnp_zeroes is that it only output zeroes that are greater than zero, but zero itself is one of the zeroes of the derivative of the zeroth-order Bessel function, $J_0'(x)$. You can see this right in the figure below:

First few Bessel functions of the 1st kind.

As seen, $J_0(x)$ attains a maximum at $x=0$, hence $x=0$ is one of the zeroes of $J_0'(x)$. Therefore, we must add it to the list, and shift the rest of the zeroes by one unit. The next step in the algorithm is to calculate the coefficients $a_{nm}$ and $b_{nm}$,

\begin{align} a_{nm} &= \frac{\langle J_0\cos(n\theta), f \rangle_w}{\langle J_0, J_0 \rangle_w} = \frac{\int_{0}^{a}\int_{0}^{2\pi}J_n\lbrac kr\rbrac\cos(n\theta) f(r,\theta)r\;drd\theta}{2\pi\int_{0}^{a}J_n^2\lbrac kr\rbrac rdr},\quad n, m = 0, 1, \ldots\\ b_{nm} &= \frac{\langle J_0\sin(n\theta), f \rangle_w}{\langle J_0, J_0 \rangle_w} = \frac{\int_{0}^{a}\int_{0}^{2\pi}J_n\lbrac kr\rbrac\sin(n\theta) f(r,\theta)r\;drd\theta}{2\pi\int_{0}^{a}J_n^2\lbrac kr\rbrac rdr},\quad n, m = 0, 1, \ldots \end{align}

We can do this using the built-in methods for double integration: scipy.integrate, dblquad:

# Define the double integrals for coefficients a_nm and b_nm
def a_nm(n, m, a):
    z_nm = bessel_derivative_zero(n, m)
    numerator = integrate.dblquad(lambda theta, r: jn(n, z_nm * r * 1.0/a) * np.cos(n * theta) * 
                        f(r, theta) * r, 0, a, lambda r: 0, lambda r: 2 * np.pi)[0]
    denominator = integrate.quad(lambda r: jn(n, z_nm * r * 1.0/a) ** 2 * r, 0, a)[0]
    return numerator / (2 * np.pi * denominator)

def b_nm(n, m, a):
    z_nm = bessel_derivative_zero(n, m)
    numerator = integrate.dblquad(lambda theta, r: jn(n, z_nm * r * 1.0/a) * np.sin(n * theta) * 
                        f(r, theta) * r, 0, a, lambda r: 0, lambda r: 2 * np.pi)[0]
    denominator = integrate.quad(lambda r: jn(n, z_nm * r * 1.0/a) ** 2 * r, 0, a)[0]
    return numerator / (2 * np.pi * denominator)
        

After calculating the coefficients, the final step is straight forward. Simply iterate through all $n$ and $m$ using a double for loops and plot the figure.

# Plotting function in polar coordinates
def plot_polar_contour(n_max, m_max, a, resolution=100):
    r = np.linspace(0, a, resolution)
    theta = np.linspace(0, 2 * np.pi, resolution)
    R, Theta = np.meshgrid(r, theta)
    
    if n_max < 1 or m_max < 1:
        print("Columns or lines cannot be less than 1\n")
        return

    # Define fig and axis
    fig, axs = plt.subplots(m_max, n_max, figsize=(1.4* n_max, 1.2 * m_max), subplot_kw=
            {'projection': 'polar'})
    
    for m in range(m_max):
        for n in range(n_max):
            Z = u(n, m, R, Theta, a)
            ax = axs[m, n] if m_max > 1 else axs[n]
            ax.contour(Theta, R, Z, levels=[0], colors=color, linewidths=thickness)
            ax.set_xticks([])
            ax.set_yticks([])
            ax.grid(False)
            
            plot_title = rf"$(n,m) = {n, m}$"
            if title:
                ax.set_title(plot_title, fontsize = 9, ha='center', rotation='vertical', 
                x=-0.1, y=0)

    plt.tight_layout()
    plt.show()

plot_polar_contour(columns, rows, a)
        

Put together all the pieces, and we may arrive at the final code for the circular plate, using the 2d wave equation to be:

import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as integrate
from scipy.special import jn, jnp_zeros

a = 1.0                 # Radius of the plate
thickness = 1           # Line thickness
color = 'darkblue'      # Color
title = True            # Title?
columns = 7
rows = 5

# Function to get the m-th zero of the derivative of the Bessel function of order n
def bessel_derivative_zero(n, m):
    if n == 0: #Redefine the first zero of J_0 prime to be 0
        zeros = list(jnp_zeros(n, m + 1))
        zeros.insert(0, 0)
        return zeros[m]
    else:
        return jnp_zeros(n, m + 1)[-1]  # Get the m-th zero, treating first zero as the 0th zero


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

# Define the double integrals for coefficients a_nm and b_nm
def a_nm(n, m, a):
    z_nm = bessel_derivative_zero(n, m)
    numerator = integrate.dblquad(lambda theta, r: jn(n, z_nm * r * 1.0/a) * np.cos(n * theta) * 
                        f(r, theta) * r, 0, a, lambda r: 0, lambda r: 2 * np.pi)[0]
    denominator = integrate.quad(lambda r: jn(n, z_nm * r * 1.0/a) ** 2 * r, 0, a)[0]
    return numerator / (2 * np.pi * denominator)

def b_nm(n, m, a):
    z_nm = bessel_derivative_zero(n, m)
    numerator = integrate.dblquad(lambda theta, r: jn(n, z_nm * r * 1.0/a) * np.sin(n * theta) * 
                        f(r, theta) * r, 0, a, lambda r: 0, lambda r: 2 * np.pi)[0]
    denominator = integrate.quad(lambda r: jn(n, z_nm * r * 1.0/a) ** 2 * r, 0, a)[0]
    return numerator / (2 * np.pi * denominator)

# 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):
    z_nm = bessel_derivative_zero(n, m)
    A_nm = a_nm(n, m, a)
    B_nm = b_nm(n, m, a)
    return jn(n, z_nm * r * 1.0/a) * (A_nm * np.cos(n * theta) + B_nm * np.sin(n * theta))

plt.rcParams['mathtext.fontset'] = 'stix'
plt.rcParams['font.family'] = 'serif'

# Plotting function in polar coordinates
def plot_polar_contour(n_max, m_max, a, resolution=100):
    r = np.linspace(0, a, resolution)
    theta = np.linspace(0, 2 * np.pi, resolution)
    R, Theta = np.meshgrid(r, theta)
    
    if n_max < 1 or m_max < 1:
        print("Columns or lines cannot be less than 1\n")
        return

    # Define fig and axis
    fig, axs = plt.subplots(m_max, n_max, figsize=(1.4* n_max, 1.2 * m_max), subplot_kw=
            {'projection': 'polar'})
    
    for m in range(m_max):
        for n in range(n_max):
            Z = u(n, m, R, Theta, a)
            ax = axs[m, n] if m_max > 1 else axs[n]
            ax.contour(Theta, R, Z, levels=[0], colors=color, linewidths=thickness)
            ax.set_xticks([])
            ax.set_yticks([])
            ax.grid(False)
            
            plot_title = rf"$(n,m) = {n, m}$"
            if title:
                ax.set_title(plot_title, fontsize = 9, ha='center', rotation='vertical', 
                x=-0.1, y=0)

    plt.tight_layout()
    plt.show()

plot_polar_contour(columns, rows, a)
        

Theoretical, Experimental Results and Data Analysis

Theoretical Results for the Square Plate

In this section, we will compare our theoretical results from the codes to the experimental results from my experiment. The designed code above should be good for up and running. Before we proceed any further, there are a few things to note. The patterns displayed below is arranged so that the frequency increases from left to right, top to bottom. Each set of patterns of the same color correspond to a different node, with different sign combinations.

The information about the modes and sign combinations is printed to the terminal. On each pattern's caption, there will be two numbers: $N=1$, which is simply the number of the pattern, and $S$, which is the sum of squares which that pattern is generated from (if $S=50$, then the modes are $(1, 7)$, $(5, 5)$ and $(7, 1)$, since $1^2+7^2=5^2 + 5^2=50$). By setting start = 1 and end = 72, we shall obtain the patterns right below:

Pattern #1 to #72.

Furthermore, if you want more information about each pattern, open the terminal and the console of this output should generate something like this:

    N = 1 | S = 0: [(0, 0)] [1]
    N = 2 | S = 1: [(0, 1), (1, 0)] [1, 1]
    N = 3 | S = 1: [(0, 1), (1, 0)] [1, -1]
    N = 4 | S = 2: [(1, 1)] [1]
    N = 5 | S = 4: [(0, 2), (2, 0)] [1, 1] 
    N = 6 | S = 4: [(0, 2), (2, 0)] [1, -1]
    N = 7 | S = 5: [(1, 2), (2, 1)] [1, 1] 
    N = 8 | S = 5: [(1, 2), (2, 1)] [1, -1]
    N = 9 | S = 8: [(2, 2)] [1]
    N = 10 | S = 9: [(0, 3), (3, 0)] [1, 1] 
    N = 11 | S = 9: [(0, 3), (3, 0)] [1, -1]
    N = 12 | S = 10: [(1, 3), (3, 1)] [1, 1] 
    N = 13 | S = 10: [(1, 3), (3, 1)] [1, -1]
    N = 14 | S = 13: [(2, 3), (3, 2)] [1, 1]
    N = 15 | S = 13: [(2, 3), (3, 2)] [1, -1]
    ... and so on

If we continue to set start = 73 and end = 144, we will obtain even more mesmerizing patterns like these:

Pattern #73 to #144.

and more, if we go further

Pattern #145 to #216.

and so much more,

Pattern #217 to #288..

There is no limit in how far you can go in a theoretical sense. But in reality, there will be no more new patterns after around $f= 5000$ Hz due to the amplitude of the plate being too small. After we are done with our theoretical result, we can compare these patterns with our actual experimental results.

For this project, I obtained the necessary equipment from my college's engineering stock room. Out of all the choices that we had, this experiment is probably the simplest to perform but quite challenging to explain and derive. If you are planning to conduct this experiment yourself, before conducting the experiment, remember to make sure that the surface on which we place our generator must be perfectly flat, or else the sand will be bouncing in a chaos manner. Wear ear plugs as well, since the plate will output some very loud noise if we are not at resonance frequency.

Turn on the wave generator, put on the ear plugs, and these patterns will arise as we increases the frequency $f$.

Experimental Results of the Square Plate.

Data Analysis of Square's Plate

The figure right below contains selected patterns that (relatively) matches with the patterns generated by our program. Note that some patterns are similar, while some other are not so much. I just pick that ones that are closest to our actual results.

Matching patterns on the Square Plate.
(Somewhat) similar patterns.

The modes and predicted wave speed of the patterns are listed in the table right below. Here, recall from the previous part that the wave speed is calculated by the formula

\begin{equation} c = \frac{2Lf}{\sqrt{n^2 + m^2}}. \end{equation}

For this experiment, the length of our plate is $24.1\pm 0.1$ cm, and we are assuming that the thickness is negligible. For better organization, I made a table below, comparing the similar patterns obtained from the actual experiment and the patterns generated by the program. Each pattern has each own modes and sign-combination, so I also listed them in the table.

Wave Speed of the Square Plate
Pattern No. $(N)$ Sum $(S)$ Modes $(n,m)$ Sign Found at Wave Speed
$21$ $20$ $(2, 4)\; (4, 2)\; $ $[+, +]$ $245\pm 50$ Hz $26.41\pm 5.39$
$17$ $16$ $(0, 4)\; (4, 0)\; $ $[+, -]$ $453 \pm 50$ Hz $54.59\pm 6.03$
$33$ $26$ $(1, 5)\; (5, 1)\; $ $[+, -]$ $56 \pm 75$ Hz $53.41\pm 7.11$
$22$ $21$ $(2, 4)\; (4, 2)\; $ $[+, -]$ $647 \pm 75$ Hz $69.73\pm 8.09$
$62$ $58$ $(3, 7)\; (7, 3)\; $ $[+, -]$ $934 \pm 100$ Hz $59.11\pm 6.33$
$52$ $50$ $(1, 7)\; (5, 5)\; (7, 1) $ $[+, +, +]$ $2160 \pm 500$ Hz $147.24\pm 34.09$
$223$ $200$ $(2, 14)\;(10, 10)\; (14, 2)\; $ $[+, +, +]$ $3163 \pm 800$ Hz $107.80\pm 27.27$

Based on the wave equation model that we used, the wave speed should remain a constant. This is not what we see here. The wave speed seems to increases as the frequency $f$ increases, and some results seem to fall off the increasing trend as well. Let's see what the circular plate has to offer.

Theoretical Results for the Circular Plate

I did not set any start and end for the circular plate, simply because they all follow the same set of rules: the number of diameter lines increases as we further move to the right, and the number of rings increases as we move downward. Set columns = 7 and rows = 5 and our output will be generated as below. The computer-generated patterns for the circular plate a bit less interesting, and is not as versatile and mesmerizing as the square patterns.

First 35 patterns of the Circular Plate.

On the other hand, our experimental results for the circular yield much more interesting outcome. I have yet to seen any diameter lines on any of my experiment, but our actual experimentation generates more patterns than just rings like our theoretical analysis.

Experimental Results of the Circular Plate.

Data Analysis of the Circular Plate

Similar to the square plate, the circular plate also has similar experimental patterns.

Matching patterns on the Circular Plate.
Similar Patterns.

From the angular frequency of the plate from the previous part, we can easily derive the wave speed equation for the circular plate:

\begin{equation} c=\frac{2\pi f r}{z_{nm}} \end{equation}

For our experiment, the plate's radius is $r = 12.1\pm 0.1$ cm.

Wave Speed of the Circular Plate
Order $(n)$ Zero ($m$) Zeroes of $J'_n(r)$ Found at Wave Speed ($m/s$)
$0$ 0 $0.000000$ $0.000$ Hz N/A
$0$ $1$ $3.831706$ $131 \pm 50$ Hz $25.885\pm 9.88$
$0$ $2$ $7.0155867$ $429 \pm 75$ Hz $46.298\pm 8.10$
$0$ $3$ $10.1734681$ $1079 \pm 100$ v $80.301\pm 14.05$
$0$ $4$ $13.3236919$ $1806 \pm 500$ Hz $102.627\pm 7.47$
$0$ $5$ $16.4706301$ $3320 \pm 750$ Hz $152.614\pm 34.50$
$0$ $6$ $19.6158585$ $4661 \pm 800$ Hz $172.184\pm 29.59$

The wave speed of the circular plate also increases as the frequency $f$ increases, which is far from our expectation. What is going on here? Well, it turns out that the wave equation model might not be sufficient after all. The problem seems to be much more complicated than I expected, and we might need a much more precise model in order to accurately describe these patterns.

Next

References

  1. Linqi. Shao. Modeling a square vibrating plate. University of Waterloo, pages 1-58, 2018.
  2. Lord Rayleigh. Vibration of Plates, pages 352-361. Princeton University Press, New York, NY, 1894.
  3. John M. Davis. Introduction to Applied Partial Differential Equation. W. H. Freeman, New York City, 1st edition, 2012.
  4. Richard. Haberman. Applied Partial Differential Equations with Fourier Series and Boundary Value Problems. Pearson, New York City, 5th edition, 2018.
  5. Rudolph. Szilard. Theories and Applications of Plates Analysis. John Wiley and Sons, Inc., 111 River Street Hoboken, NJ 07030-5774, 2004
  6. Dispersion of flexural waves. https://www.acs.psu.edu/drussell/Demos/Dispersion/Flexural.html, 2004. Accessed: 2024-05-26
  7. Thomas D. Rossing and Neville H. Fletcher. Two-Dimensional Systems: Membranes and Plates, pages 65-94. Springer New York, New York, 2004.
  8. Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio. Stable and accurate numerical methods for generalized kirchhoff-love plates, 2020. Accessed: 2024-05-26.

More Articles