String Modes

Columbia University

This notebook continues our study of waves and oscillations using our numerical model of a uniform string. The goal of this notebook is to understand the relationship between the structure of oscillations on a string and the frequency of oscillations on a string. This relationship is called a dispersion relation.

We will necessarily discuss the important concept of eigensystems, and we will use the built-in functions Eigensystem and Fourier to analyze the normal modes of our model string.

Our Model String

Last week, we constructed a model of a "string" consisting of a collection of "balls" of mass, m, connected by "springs" with Hook constants, k. For N balls, we wrote 2N simplified equations for 2N unknown functions of time: the ball positions, y[i][t], and vertical velocities, vy[i][t].

The equations of motion are:
	(∂_t) y[i][t]  =  vy[i][t]  (* for all i *)
	(∂_t) vy[1][t]  =  - 2(ω^2) y[1][t] +  (ω^2)  y[2][t]
	(∂_t) vy[i][t]  =    (ω^2)  y[i-1][t] - 2(ω^2) y[i][t] +  (ω^2)  y[i+1][t]
	(∂_t) vy[N][t]  =  (ω^2)  y[N-1][t] - 2(ω^2) y[N][t] 
In these equations, the characteristic frequency, ω, is the frequency of a single ball connected to fixed neighbors with two springs. For simplicity, we normalize time to 2π/ω which is equivalent to assigning ω-> 1.

Equations of Motion (N balls)

nBalls = 30;

p = Table[{y[i][t],vy[i][t]}, {i,1,nBalls}] // Flatten;

dp = {vy[1][t],-2 y[1][t] + y[2][t]}~ Join ~  
Table[{vy[i][t],y[i-1][t]-2 y[i][t] + y[i+1][t]}, {i, 2, nBalls-1}] ~Join~
{vy[nBalls][t], y[nBalls-1][t]-2 y[nBalls][t]} // Flatten;

A particular initial condition

Defining an initial condition for the string...

i0 = 10;
di = 4;
initial[i_] := Exp[-(i-i0)^2/(di^2)]

ListPlot[Table[initial[i], {i,1,30}], 	PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None]

p0 = Table[{y[i][0] == initial[i], vy[i][0]==0}, {i,1,nBalls}] // Flatten;

Example solution

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

ballN = NDSolve[eqs , p, {t, 0, 100}] // Flatten;

Plot[Evaluate[y[3][t] /. ballN], {t, 0, 100}]

Animation

ListPlot3D[Table[y[i][t] /. ballN,{i,1,nBalls}, {t,0,100,1}]]

ListAnimate[Table[ListPlot[Evaluate[Table[y[i][t],{i,1,nBalls}]/. ballN],
			PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None,
			DisplayFunction->Identity],{t,0,100,0.5}], AnimationRunning->False]

Exciting Normal Modes

The frequency of oscillations of a string are very dependent on the initial deflection of the string. In particular, the frequencies are inversely proportional to the "wavelength" of the oscillation. In this section, we experiment with our numerical string and determine the relationship between the wavelength and the frequency of oscillation.

Sine[π m x] modes

As you might guess, the natural modes of a uniform string are sine waves. They can be excited by choosing an initial condition proportional to Sin[π m i/(nBalls+1)]; where m = 1, 2, 3, 4, 5, ...  Notice, these modes form a discrete set of oscillations whose wavelength decreases as m increases. The wavelength increases in particular integral steps because the oscillations must satisfy boundary conditions on both ends of the string!

m = 1

m = 1;
initial[i_] := Sin[π m i/(nBalls+1)]

ListPlot[Table[initial[i], {i,1,30}], 	PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None]

p0 = Table[{y[i][0] == initial[i], vy[i][0]==0}, {i,1,nBalls}] // Flatten;

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

ballN = NDSolve[eqs , p, {t, 0, 100}] // Flatten;

Plot[Evaluate[y[3][t] /. ballN], {t, 0, 100}]

m = 2

m = 2;
initial[i_] := Sin[π m i/(nBalls+1)]

ListPlot[Table[initial[i], {i,1,30}], 	PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None]

p0 = Table[{y[i][0] == initial[i], vy[i][0]==0}, {i,1,nBalls}] // Flatten;

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

ballN = NDSolve[eqs , p, {t, 0, 100}] // Flatten;

Plot[Evaluate[y[3][t] /. ballN], {t, 0, 100}]

m =5

m = 5;
initial[i_] := Sin[π m i/(nBalls+1)]

ListPlot[Table[initial[i], {i,1,30}], 	PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None]

p0 = Table[{y[i][0] == initial[i], vy[i][0]==0}, {i,1,nBalls}] // Flatten;

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

ballN = NDSolve[eqs , p, {t, 0, 100}] // Flatten;

Plot[Evaluate[y[3][t] /. ballN], {t, 0, 100}]

Animation

ListAnimate[Table[ListPlot[Evaluate[Table[y[i][t],{i,1,nBalls}]/. ballN],
			PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None,
			DisplayFunction->Identity],{t,0,100,0.5}], AnimationRunning->False]

Mode Combinations

If the wavelength is slightly different than π m/(nBalls+1), then the string oscillation will not be uniform. Multiple modes are excited.

Almost m = 1

m = 1;
initial[i_] := Sin[0.9π m i/(nBalls+1)]

ListPlot[Table[initial[i], {i,1,30}], 	PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None]

p0 = Table[{y[i][0] == initial[i], vy[i][0]==0}, {i,1,nBalls}] // Flatten;

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

ballN = NDSolve[eqs , p, {t, 0, 100}] // Flatten;

Plot[Evaluate[y[3][t] /. ballN], {t, 0, 100}]

Animation

ListAnimate[Table[ListPlot[Evaluate[Table[y[i][t],{i,1,nBalls}]/. ballN],
			PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None,
			DisplayFunction->Identity],{t,0,100,0.5}], AnimationRunning->False]

Combining m=1 and m=4

m1 = 1;m2=4;
initial[i_] := 0.7 Sin[π m1 i/(nBalls+1)]+0.3 Sin[π m2 i/(nBalls+1)]

ListPlot[Table[initial[i], {i,1,30}], 	PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None]

p0 = Table[{y[i][0] == initial[i], vy[i][0]==0}, {i,1,nBalls}] // Flatten;

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

ballN = NDSolve[eqs , p, {t, 0, 100}] // Flatten;

Plot[Evaluate[y[3][t] /. ballN], {t, 0, 100}]

Animation

ListAnimate[Table[ListPlot[Evaluate[Table[y[i][t],{i,1,nBalls}]/. ballN],
			PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None,
			DisplayFunction->Identity],{t,0,100,0.5}], AnimationRunning->False]

String EigenSystem: Wavelength & Frequency

As can be seen from the examples in the previous section, the frequency of a string oscillation depends upon the spatial structure or wavelength  of the oscillation. This is a fundamental property of waves, and it is called wave dispersion. The relationship between the frequency and wavelength of oscillations is called the dispersion relation. For light, f = c/λ, where c is the speed of light. For sound, f = (Subscript[c, s])/λ, where (Subscript[c, s]) is the speed of sound. These linear dispersion relations with constant wave propagation speeds are the most simple found in nature. In general, the relationship between frequency and wavelength is very complex.

Mathematica has several built-in algebraic tools to determine the dispersion relation (i.e. the relationship between frequency and wavelength) of our string of balls and springs. 
The general mathematical subject required is call eigensystem analysis–an essential part of linear algebra.
We can express all 2N coupled equations of motion as a single matrix equation. Mathematica can then determine all of the frequencies for each independent mode of oscillation of the string. The frequencies are called eigenvalues and the spatial structures of each mode are called eigenvalues.

Writing Equations of Motion as a Matrix

How do we write our equations of motion as a matrix equation?
The answer is simple: we've already done so.

The equations for the positions of the balls are written in matrix form as:
	(∂_t)y = v
where y = {y[1][t], y[2][t], y[3][t],...y[nBalls][t]}and v = {vy[1][t], vy[2][t], vy[3][t],...vy[nBalls][t]}.
The equations for the rate of change of the velocities of the balls are written as:
	(∂_t)v = - (ω^2) M . y
where M is a simple matrix which expresses how each ball is connected by two springs to its neighboring balls.
Since (ω^2) and M are constants, both equations can be combined:
	(∂_t) ∂Nully = -  (ω^2) M . y

We can write down M very easily in Mathematica:

nBalls = 30;

mMatrix = Table[Which[j == i, 2, j == i+1, -1, j == i-1, -1, True, 0], {i,1,nBalls}, {j, 1, nBalls}];

mMatrix  // Shallow

Using Eigensystem

Eigensystem solves the matrix equation:
	(α^2) v = M . v
where (α^2) is an eigenvalue and v is an eigenvector. The matrix, M, must be square, meaning that the number of rows must equal the number of columns. If the matrix is N × N, then there will be N eigenvalues and N corresponding eigenvectors. These related to the frequencies and structures (i.e.  wavelength) of the modes on the string. For our problem, α is proportional to the normalized frequency of the eigenmode. This can be seen by substitution
	(∂_t) (∂_t) y = -  (ω^2) M . y = -  (((ωα))^2) y
which has simple harmonic solutions proportional to y[i][t] ∝ Sin[ωα t].

In general, it is very time-consuming to evaluate the eigenvalues and eigenvectors of a large matrix. To help Mathematica, be sure to use approximate numerical computations.

{values, vectors} = Eigensystem[mMatrix // N];

Mathematica orders the modes from the highest frequency to the lowest:

ListPlot[values]

The lowest, or fundamental, frequency is much lower than the natural frequency of a single ball (which we have normalized to unity):

Sqrt[values[[30]]]

The corresponding mode structure is

ListPlot[vectors[[30]],
			PlotRange->{{0,nBalls+1},{-0.5, 0.5}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None]

The first six modes are

modes = Table[ListPlot[vectors[[i]],
			PlotRange->{{0,nBalls+1},{-0.5, 0.5}},
			PlotStyle->{Hue[(30-i)/5],PointSize[0.02]},
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None,
	DisplayFunction-> Identity],{i,30,25,-1}]

Show[modes, DisplayFunction->$DisplayFunction]

Finding the Wavelength

We can use the built-in function, Fourier, we can find the wavelength of each eigenvector. Remember, Fourier produces a list of the "frequencies" contained in its argument. The wavelength is analogous to the temporal period of an oscillation. Fourier will produce a number proportional to the inverse of the wavelength.

A function to find the dominate wavelength is:

wavelength[i_Integer] := Module[{absFourier},
		
		(* Find the Fourier transform of a double-length eigenvector *)
		absFourier = Abs[Fourier[vectors[[i]]~Join~(If[OddQ[i],1,-1] vectors[[i]])]];
		
		(* Find the average wavelength = inverse of the average mode index from Fourier *)
		(Plus @@ absFourier[[Range[1,nBalls]]])/(Plus @@ Table[(j-1) absFourier[[j]], {j,1,nBalls}])
	
	]

This function is interesting. First, it joins together two copies of the ith mode. When the mode index is even, it reverses the sign of the copy to be joined so that the oscillation structure is continuous.
For example, the two lowest modes are shown below after joining together a proper copy:

ListPlot[vectors[[29]]~Join~(If[OddQ[29],1,-1] vectors[[29]])]

ListPlot[Take[Abs[Fourier[Join[vectors[[29]],If[OddQ[29],1,-1] vectors[[29]]]]],{1,30}],Joined->True,PlotRange->{0,1}]

ListPlot[vectors[[30]]~Join~(If[OddQ[30],1,-1] vectors[[30]])]

ListPlot[Take[Abs[Fourier[Join[vectors[[30]],If[OddQ[30],1,-1] vectors[[30]]]]],{1,30}],Joined->True,PlotRange->{0,1}]

After forming the copy, Fourier is used to compute the spectrum of the spatial oscillations. We use the Plus function to find the average spatial oscillation rate produced by Fourier. The wavelength is the inverse of this rate.
The first few wavelengths are:

Table[wavelength[i],{i,30,25,-1}]

Finding the Dispersion Relation

If you recall, for simple oscillations, the frequency and wavelength are related by a propagation speed,
	f = c/λ
This is exactly what we see...

ListPlot[Table[{wavelength[i],(Sqrt[values[[i]]])},{i,30,15,-1}],Joined->True,AxesLabel->{"Wavelength","Frequency"}]

The wavenumber is proportional to the inverse of the wavelength:

ListPlot[Table[{(1/(wavelength[i])),(Sqrt[values[[i]]])},{i,30,15,-1}],Joined->True,AxesLabel->{"Wavenumber","Frequency"}]

This is a very important figure. It shows that oscillations of all wavelengths propagate at the same speed!

This is just what we have already seen in our numerical experiments. When we excite perturbations, no matter what their shape, they move along the string without distortion. However, when they reach the ends of the string, there is a reflection, and oscillations interfere with each other as they bounce back and forth along the string.