Departments of Civil Engineering and Mathematics
CE 109: Computing for Engineering

Mathematica Session 2:
Movies

1. Vibrating strings

In Session 1, you were briefly introduced to some of Mathematica's capabilities for generating mathematical graphics: the functions Plot, ParametricPlot, Plot3D and so on. These functions can, of course, be used to generate single, static diagrams, and that's what you saw in Session 1. But Mathematica also has the capability to generate and display moving, animated graphics. These are called "movies", and indeed they work on the same principle as the cinema: the user creates a series of still pictures (frames) and uses Mathematica to display them in rapid succession, creating the illusion of motion.

We're going to use this feature of Mathematica to study the engineering topic of standing vibrations: in rod- or string-like structures to begin with, then in disc- or drum-like ones. In a rod or string, many different types of vibration are possible. Think of the different sounds made by a guitar (in which the string is usually plucked), a piano (in which it is struck sharply) and a violin (in which it is subjected to a continual driving vibration from the bow): the differences in sound quality among these instruments are partly due to the different types of vibration response these three forms of stimulus produce in strings. But it turns out, as you'll discover in the Second Year mathematics course, that all these types of vibration, indeed every type of vibration, can be broken down into a sum of simple sinusoidal components, called the normal modes.

In a sinusoidal standing vibration, with no damping, the displacement y m of a point x m along the string at time t s is given by

\(TraditionalForm
    \`\(y = a\ \(sin(λ\ x)\)\ \(cos(c\ λ\ t + α)\), 
    \)\)

where a, λ , c and α are constants. For simplicity, let's assume that a and c are both 1 and that α is zero, giving

y =  sin(λ x) cos(λ t).

We assume, too, that the string is of length 1 m and that the ends of the string are fixed; in other words, that y is equal to zero at x = 0 and x = 1. It follows from this that λ must be a multiple of π (why?).

In Mathematica, we can set up the simplest normal mode, in which λ=π, by typing the following:

wave1[x_, t_]  :=  Sin[π x] Cos[π t]

(use the palette on the right of the screen for Greek letters, or type <Esc>pi<Esc>).

Notice the use of the combined "colon-equals" sign, :=, and (on the left-hand-side)  the underscore character, _. This is an example of a new kind of Mathematica command, called a function definition: it doesn't produce any output itself, but instead "teaches" Mathematica a new function, called wave1, to add to the ones it already knows like Sqrt, Sin and so on. We can test this new function by trying commands like

wave1[1/3,  2/3]

wave1[1,  1/4]wave1[0.5,  1.75]

etc.

We can get a look at the wave at time t = 0 like this:

Plot[Evaluate[wave1[x,  0]],  {x,  0,  1}]

We can examine the wave for other values of t, for example:

Plot[Evaluate[wave1[x,  0.3]],  {x,  0,  1}]

The problem, as you can see, is that Mathematica's automatic choice of plotting range makes this diagram look very similar to the previous one. We can get round that by overriding this automatic choice:

Plot[Evaluate[wave1[x,  0.3]],  {x,  0,  1},  
      PlotRange  ->  {-1, 1}]

Now for the animation. The following command generates 24 frames:

Table[ 
      Plot[Evaluate[wave1[x,  t]],  {x,  0,  1},  
        PlotRange  ->  {-1, 1}], {t,  0,  23/12,  1/12}];

Note that one complete cycle, with these choices for our constants, takes 2 seconds, since

\(TraditionalForm
    \`y(x, 0) = 
      \(\(sin(π\ x)\)\ \(cos(0)\)
	\ \ \ \ \ \ \ \ \ \ \ \ \  = 
        \(sin(π\ x)
	\ \ \ \ \ \ \ \ \ \ \ \ \  = 
          \(\(sin(π\ x)\)\ \(cos(2\ π)\)
			\ \ \ \ \ \ \ \ \  = 
            \(y(x, 2).\)\)\)\)\)

(Of course, this doesn't mean that the animation will take 2 seconds: we can play the "movie" at any speed we like.) You see, though, that in creating the animation frames we've only gone up, in steps of 1/12 s, as far as t = 23/12. Why do you think we haven't gone all the way to t = 2?

To animate these frames, the first thing you have to to is double-click on the cell bracket that contains them all (on the right-hand edge of the Mathematica window). They should all disappear apart from the first frame in the sequence (consult a tutor if this isn't working for you). Then, either select Animate Selected Graphics from the Cell menu or (quicker) hold down <Ctrl> and type Y.

You can change the speed or direction of the movie using the buttons at the bottom left of the screen. To stop, click the mouse button.

(i) What if we hadn't set α equal to zero, but to some other value? Generate and view an animation of 
wave1a[x, t], where, say,

wave1a[x_,  t_]  :=  Sin[π x] Cos[π t  +  1.9]

(ii) Generate and view an animation of the next simplest normal mode, beginning this time with

wave2[x_,  t_]  :=  Sin[2 π x] Cos[2 π t]

What's different about this one? How long does a single cycle of this vibration take? Look at wave3, at wave4, and so on.

(iii) Experiment with adding multiples of your normal modes together, for example

combinedWave[x_,  t_]  :=  
      wave1[x, t]  +  wave2[x, t]/2  +  wave3[x, t]/3

2. Vibrating drums (optional)

The mathematics in this section is slightly more advanced than that in the other two. It's worth it for the brilliant 3D animations you get, though. However, if you're running short of time you may want to skip this section and perhaps come back to it later. Animations use up a lot of memory: to clear the memory you've used in Section 1, it may be worth quitting Mathematica and relaunching from scratch.

If we study a drum instead of a string, our animated graphics will, of course, be three-dimensional, so we'll have to use Plot3D instead of Plot. The mathematics of three-dimensional vibration is a bit more advanced. The normal modes of a circular drum aren't sinusoidal, but instead involve a special mathematical function (called a Bessel function) that can't be expressed in terms of the familiar functions (sin x, cos x, (e^x) etc.) you've already met. The equation that describes symmetrical normal-mode standing vibrations in a drum is

\(TraditionalForm
    \`z = a\ \(\(J\_0\)(λ\ r)\)\ 
        \(cos(c\ λ\ t\  + \ α)\)\)

where (Subscript[J, 0]) is the Bessel function mentioned above, and r = (Sqrt[(x^2) + (y^2)]). Making the same simplifying assumptions as we did for the string, this becomes

\(TraditionalForm
    \`z = \(\(J\_0\)(λ\ r)\)\ \(cos(λ\ t)\)\)

If we assume that the drum has radius 1 m, then it turns out that the simplest normal mode corresponds to λ = 2.4048. We set this up in Mathematica as follows:

drumWave1[r_,  t_]  :=  BesselJ[0,  2.4048 r] Cos[2.4048 t]

We can see how things look at time t = 0 by typing:

Plot3D[Evaluate[drumWave1[(Sqrt[(x^2) + (y^2)]),  0]],  {x, -1, 1}, 
       {y, -1, 1}]

(i) With λ = 2.4048, how long is one cycle of this vibration? How long is one twenty-fourth of one cycle? Generate a 24-frame animation of the vibration.

(ii) The next two normal modes of the drum correspond to λ = 5.5201 and λ = 8.6537. Animate each of these normal modes. Study some combinations of them, like you did for the string. Why do you think stringed instruments have pure, harmonious tones whereas drums, gongs and cymbals do not?

(iii) (Advanced, optional.) Plot the Bessel function (Subscript[J, 0])(r)  for values of r between 0 and, say, 15. What do you notice about the values of r for which the curve crosses the horizontal axis? Can you find these values? (Check out the Mathematica function FindRoot in the Help System.) Find some more zeroes for this function, and use them to find some more normal modes for the drum. Animate these.

3. Beams

N.B. Animations use up a lot of memory: to clear the memory you've used in Section 1 or Section 2, it may be worth quitting Mathematica and relaunching from scratch.

The following diagram shows a cantilever, supported at the left-hand end, and subjected to a concentrated load of magnitude p N at a point d m along its length. The vertical scale has been exaggerated. By analysing bending moments it's not hard to show that for a point x m from the support, and to the left of the load (i.e.  x < d),   the displacement v(x)  has a cubic profile, given by the equation

v=-((p (d - x/3) (x^2))/(2 EI)),

where p is the load in N and EI is the flexural stiffness in N (m^2). However, this cubic equation is not valid to the right of the load (x > d). Here, the beam is straight in profile, as shown in the diagram (the graph of the cubic equation for this region is shown, for purposes of contrast, as a thin dashed line). The straight line segment that corresponds to the beam in this region has the same gradient as the cubic at the point x = d, and therefore corresponds to the tangent at this point.

(i) Determine the formula for the displacement for x > d. If you wish to use Mathematica for this step, recall from Session 1 that the D command calculates derivatives, and that the "/." operation can be used to substitute values into expressions, as in

-((p (d - x/3) (x^2))/(2 EI)) /. x -> d

Also recall that the equation for a straight line of gradient m through a given point ((Subscript[x, 0]), (Subscript[v, 0])) is v = m(x - (Subscript[x, 0])) + (Subscript[v, 0]).   Try using the command Simplify on your expressions if they seem to contain terms that will simplify.

(ii) You now have two expressions for the displacement v, one that applies to the left of the load (for x < d) and one that applies to the right (for x > d). Combine them into one displacement function by typing, and completing, the following command:

v[p_, d_, x_] := If[x < d, -((p (d - x/3) (x^2))/(2 EI)), ({{... fill in your }, {expression here...}})]

(iii) Try out your displacement function by typing, for example

v[5000,  1.2,  0.8]v[p,  1.1,  2.7]

(iv) Set EI=5000000 and plot the graph of the beam displacement for p = 5000, d = 1.4 for a beam of length 
3 m, by typing

Plot[v[5000,  1.4,  x],  {x,  0,  3}]

Plot some examples of your own, for different values of p and d.

(v) Generate an animation of a 3 m beam under a slowly* increasing load applied at 1.2 m by typing

Table[Plot[v[p,  1.2,  x],  {x,  0, 3}],  {p,  0,  10000,  500}]

If you examine your animation frames carefully, there's something very wrong with them! Change the above command to put this right, then generate the animation again and view it.

(vi) Generate and view a 31-frame animation of a beam under a 5 kN load that moves slowly* from one end of the beam to the other.

*"Slowly" in this context means you can ignore any vibrational effects.