Circular membrane vibration simulation

Question

I'm new in Mathematica and I'm trying to simulate the vibration of a circular membrane for math project but I don't even know how to start.

The wave equation describes the displacement of the membrane $(z)$ as a function of its position $(r,\theta)$ and time $(t)$.

$$ \frac{\partial^2 z}{\partial t^2}=c^2 \nabla^2 z $$

For an ordinary circular drumhead the imposed boundary conditions are that $z$ is $2π$ periodic, finite within $0 ≤ r ≤ a$ and at $r = a$ (the edge of the drum) the displacement is zero: $z(a) = 0$. The wave equation is solved through a method of separation of variables and produces the solution:

$$ z=A_{nj}J_n\left(\frac{\omega r}{c}\right)\sin(\omega t+\phi)\sin(n\theta+\varphi) $$

How can I simulate/plot this solution in Mathematica?

Answer 1

The separation-of-variables solution you quoted has two indices appearing in it: n and j (the subscripts of the coefficient $A_{nj}$). Here, n is azimuthal mode order, i.e. it counts the number of nodes along the direction in which the polar-angle $\theta$ varies (divided by 2).

The index j is needed because the wave is supposed to satisfy the boundary condition of being zero at the radius a. Here, it's best to employ a as the unit of length, so that r=1 is the radius of the circle, and the boundary condition becomes

$$J_{n}(\frac{\omega_j}{c}) = 0$$

Here, I added the index j to the frequency, because only a discrete set of $\omega_j$ can satisfy the above equation.

To determine these allowed frequencies, you can use BesselJZero.

For the plot, I'll convert the polar coordinate form of the separated solution to Cartesian coordinates. This is done by defining the function fXY below. It takes the index n and the wave number $k_j\equiv\omega_j/c$ as inputs. I'll leave out the time dependent factor for now, i.e., consider only the spatial variation of the wave at a given fixed time. Also, I'll choose the phase $\phi$ of the azimuthal solution to be such that I get a cosine instead of a sine:

fXY[n_, k_][x_?NumericQ, y_?NumericQ] := 
 BesselJ[n, k Sqrt[x^2 + y^2]] Cos[n ArcTan[y, x]]

Now do the plot, given a pair of indices {n, j}:

wavePattern[n_, j_] := Module[
  {k0 = N[BesselJZero[n, j]]},
  DensityPlot[fXY[n, k0][x, y],
   {x, -1.1, 1.1}, {y, -1.2, 1.1},
   RegionFunction -> Function[{x, y}, x^2 + y^2 < 1],
   ColorFunction -> "BlueGreenYellow",
   PlotPoints -> 100,
   MaxRecursion -> 0, 
   Epilog -> 
    Inset[Grid[{{"n", n}, {"j", j}}, Frame -> All], {-.9, .9}],
   BaseStyle -> {FontFamily -> "Arial"}]
  ]

I count the index j starting at 1 because there is always at least one radial node (at the boundary).

As an example, here are some plots:

Show[GraphicsGrid@Table[wavePattern[m, n], {m, 0, 2}, {n, 1, 3}], 
 ImageSize -> 700]

In the DensityPlot inside the function wavePattern, I set the option MaxRecursion -> 0 to speed up the plotting. An additional ingredient in making the plot for this solution is the use of RegionFunction to define the circular domain of the wave.

Since I left out the time dependence, it may be better to apply a color scheme that emphasizes the nodal lines, because these lines are the only thing that stays constant in time for a standing wave such as this. So here is an alternative plotting function:

wavePattern[n_, j_] := 
 Module[{k0 = N[BesselJZero[n, j]]}, 
  DensityPlot[fXY[n, k0][x, y], {x, -1.1, 1.1}, {y, -1.2, 1.1}, 
   RegionFunction -> Function[{x, y}, x^2 + y^2 < 1], 
   ColorFunction -> 
    Function[{x}, 
     Blend[{White, Darker@Brown, White}, 2 ArcTan[10 x]/Pi + .5]], 
   ColorFunctionScaling -> False, PlotPoints -> 100, 
   MaxRecursion -> 0, 
   Epilog -> 
    Inset[Grid[{{"n", n}, {"j", j}}, Frame -> All], {-.9, .9}], 
   BaseStyle -> {FontFamily -> "Arial"}]]

Show[GraphicsGrid@Table[wavePattern[m, n], {m, 0, 2}, {n, 1, 3}], 
 ImageSize -> 700]

This is pretty close to the kind of pattern you would actually observe if you took a vibrating plate and spread sand on it: the grains will move to the places in the standing wave where the amplitude remains zero at all times, i.e., the nodes. That's what is shown in brown above.