1D Normal Modes¶

Once again I've been inspired by my friend Rhett Allain who has been working on the 2 mass/3 spring eigenvalue problem. Check out his most recent video.

What I want to do is something I've blogged a little about before, but in 1D it ought to be a little easier.

The main approach is to set up a system and then shake it a tiny bit at either a single frequency (which we'll adjust) or a bunch of frequencies at once (that makes the analysis harder). What I'll do is imagine Rhett's system where the so-called fixed red dot on the left will actually be driven wiht a tiny amplitude.

Ok, let's get started. I'll do the usual Lagrangian approach and we'll drive the left edge at a single sinusoidal frequency to start. Note that x1[t] and x2[t] will be measured from their natural starting position:

In [1]:

k1=1;
k2=2;
k3=3;
m1=1;
m2=2;
l1=l2=l3=1;
drive=1;
left[t_]:=0.1 Sin[drive t];
KE=1/2 m1 x1'[t]^2+1/2 m2 x2'[t]^2;
PE=1/2 k1 (x1[t]-left[t])^2+1/2 k2 (x2[t]-x1[t])^2+1/2 k3 (x2[t])^2;
L=KE-PE;
el[a_]:=D[L, a[t]]-D[L, a'[t],t]==0;
sol=First[NDSolve[{el/@{x1,x2},
        x1[0]==0,
        x2[0]==0,
        x1'[0]==0,
        x2'[0]==0},{x1,x2}, {t,0,tmax=100}]];

Here's a plot of the two variables (x1 in blue, x2 in orange):

In [14]:

Plot[{x1[t]/.sol, x2[t]/.sol}, {t,0,tmax}]

Out[14]:

Here's a movie of that:

In [15]:

frames=Table[Graphics[{PointSize[0.05], 
                        Red,
                        Point[{left[t],0}/.sol],
                        Black,
                        Point[{l1+x1[t],0}/.sol],
                        Point[{l1+l2+x2[t],0}/.sol],
                        Red,
                        Point[{l1+l2+l3,0}],
                        Green,
                        Line[{ {left[t],0},{l1+l2+l3,0}}]}, PlotRange->{ {-0.2, l1+l2+l3+0.2}, {-1,1}}],
                    {t,0,tmax,tmax/100}];
ListAnimate[frames]

Out[15]:

Ok, that's a little fast, but you get the point.

So, to find the resonances, we could do a couple of different things:

1 - Try a bunch and see who causes the biggest response¶

I'll turn my ODE solver code into a function that returns the highest amplitue of x1 and x2 (maybe summed?) so that I can hunt for big ones:

In [34]:

Clear[singleFrequency]
singleFrequency[{k1_,k2_,k3_}, {m1_,m2_}, {l1_,l2_,l3_}][drive_?NumericQ]:=(
    left[t_]:=0.1 Sin[drive t];
    KE=1/2 m1 x1'[t]^2+1/2 m2 x2'[t]^2;
    PE=1/2 k1 (x1[t]-left[t])^2+1/2 k2 (x2[t]-x1[t])^2+1/2 k3 (x2[t])^2;
    L=KE-PE;
    el[a_]:=D[L, a[t]]-D[L, a'[t],t]==0;
    {sol, {steps}}=Reap[First[NDSolve[{el/@{x1,x2},
            x1[0]==0,
            x2[0]==0,
            x1'[0]==0,
            x2'[0]==0},{x1,x2}, {t,0,tmax=20},
            StepMonitor :> Sow[{Abs[x1[t]], Abs[x2[t]]}]]]];
    Max[steps[[All,1]]]+Max[steps[[All,2]]]
    )

In [33]:

singleFrequency[{1,2,3}, {1,2},{1,1,1}][1]

Out[33]:

0.439052

In [29]:

Plot[singleFrequency[{1,2,3}, {1,2},{1,1,1}][omega], {omega,0,10},PlotRange->All]

Out[29]:

Aha! You can see two clear resonances! So this works pretty well, though it's a little slow (that cell took about a minute to run). Let's see what the two peaks look like:

In [39]:

first=FindMaximum[singleFrequency[{1,2,3}, {1,2},{1,1,1}][omega],{omega,1.1}]
second=FindMaximum[singleFrequency[{1,2,3}, {1,2},{1,1,1}][omega],{omega,2.0}]

FindMaximum::lstol: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient increase in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

FindMaximum::sdprec: Line search unable to find a sufficient increase in the function value with MachinePrecision digit precision.

Out[39]:

Ok, now we know the omegas that lead to a peak. Let's look at the mass locations for those:

In [41]:

singleFrequency[{1,2,3}, {1,2},{1,1,1}][1.12989]
Plot[{x1[t]/.sol, x2[t]/.sol}, {t,0,20}]

Out[41]:

They're growing together and in phase, which is usually what you expect for the lower resonance.

Now the second one:

In [43]:

singleFrequency[{1,2,3}, {1,2},{1,1,1}][2.04046]
Plot[{x1[t]/.sol, x2[t]/.sol}, {t,0,20}]

Out[43]:

As expected, they're both growing and are 180 degrees out of phase. Cool!

Here's a movie of both of those:

In [82]:

tmp=singleFrequency[{1,2,3}, {1,2},{1,1,1}][1.12989];
sol1=sol;
tmp=singleFrequency[{1,2,3}, {1,2},{1,1,1}][2.04046];
sol2=sol;
singleFrame[sol_,shift_,omega_][t_]:=Graphics[{PointSize[0.05], 
                        Red,
                        Point[{0.1 Sin[omega t],shift}/.sol],
                        Black,
                        Point[{l1+x1[t],shift}/.sol],
                        Point[{l1+l2+x2[t],shift}/.sol],
                        Red,
                        Point[{l1+l2+l3,shift}],
                        Green,
                        Line[{ {0.1 Sin[omega t],shift},{l1+l2+l3,shift}}]}, PlotRange->{ {-0.2, l1+l2+l3+0.2}, {-1,2}}]
frames=Table[Show[singleFrame[sol1,0, 1.12989][t], singleFrame[sol2,1, 2.04046][t]],{t,0,tmax,tmax/100}];
ListAnimate[frames]

Out[82]:

2 - Excite them with all frequencies at once¶

Here I'll move the left red ball quickly in a linear fashion and then quickly return it to zero. Then I'll Fourier Transform the results to see the frequency spectrum

In [184]:

k1=1;
k2=2;
k3=3;
m1=1;
m2=2;
l1=l2=l3=1;
drive=1;
left[t_]:=Piecewise[{ {10 t, t<0.01},{0.1-10(t-0.01), t<0.02}}];
KE=1/2 m1 x1'[t]^2+1/2 m2 x2'[t]^2;
PE=1/2 k1 (x1[t]-left[t])^2+1/2 k2 (x2[t]-x1[t])^2+1/2 k3 (x2[t])^2;
L=KE-PE;
el[a_]:=D[L, a[t]]-D[L, a'[t],t]==0;
sol=First[NDSolve[{el/@{x1,x2},
        x1[0]==0,
        x2[0]==0,
        x1'[0]==0,
        x2'[0]==0},{x1,x2}, {t,0,tmax=100}]];

Here's a plot of the left edge zoomed in to the first tenth of a second of time:

In [121]:

Plot[left[t],{t,0,.1},PlotRange->All]

Out[121]:

and here's a plot of x1 and x2

In [155]:

Plot[{x1[t]/.sol,x2[t]/.sol},{t,0,tmax}]

Out[155]:

They don't have a lot of amplitude, but it's enough for a Fourier Transform:

In [171]:

xdata=Table[x1[t]/.sol,{t,0,tmax,tmax/100}];
ydata=Table[x2[t]/.sol,{t,0,tmax,tmax/100}];
xFT=Fourier[xdata];
yFT=Fourier[ydata];
ListLinePlot[{Abs[xFT]^2, Abs[yFT]^2},PlotRange->All]

Out[171]:

See, two peaks! (I know, I know, it looks like four, but that's because of the dumb way most fft (fast fourier transforms) work). Really you should just show the first half of the data if the original data is composed of real numbers, which this data set is:

In [183]:

ListLinePlot[{Abs[xFT]^2, Abs[yFT]^2},PlotRange->{ {0,50}, All}]

Out[183]:

There, see, 2 peaks!

So how do you get the properly scaled frequencies from a discrete fft? Well, every point represents a frequency, with the first one being zero and then they're all spaced by $\frac{2\pi}{\Delta t}$ where $\Delta t$ it to total time that the data spans.

So first we have to find the peak locations:

In [177]:

FindPeaks[Abs[xFT]^2]

Out[177]:

                -6                 -6                 -6                  -6
{{19, 1.76219 10  }, {34, 1.9945 10  }, {69, 1.9945 10  }, {84, 1.76219 10  }}

So we have a peak in the 19th and 34th position. That means that the first one is 18 of the frequency spacings:

In [181]:

2 Pi (19-1)*1./(tmax)

Out[181]:

1.13097

and the second on is at 33 of the frequency spacings:

In [182]:

2 Pi (34-1)*1./(tmax)

Out[182]:

2.07345

They're both pretty darn close to what we found before, and we did by just doing the ode solver once!

Your thoughts?¶

I'd love to hear your thoughts. Here are some starters for you:

1. This is great! I especially like . . .
2. This is dumb! What's truly dumb is . . .
3. I saw Rhett posted about a 3-mass system on twitter. Why are you so slow?
4. I don't understand why the Fourier command duplicates the data for real values
5. Wait, the Wolfram Engines starts counting lists at one and not zero? That's crazy!
6. What's with those warnings after the FindMaximum commands?
7. Why do you do option 2 with a triangle pulse? Would something else work better?
8. I don't think summing the max amplitudes of x1 and x2 is the right way to measure the resonances. A better way would be . . .