Modeling the chaotic motion apparatus

Peter N. Saeta
20 February 2011

Mathematica makes it possible to model the behavior of a mechanical system, once you figure out the appropriate equations of motion. This file shows how to do it, using a different mechanical system to show you how to do it. You will need to edit this document to change the model to match the equations of motion you derive for the chaotic motion apparatus.
If you have never used Mathematica before, please take a few moments to watch a video introduction at http://www.wolfram.com/broadcast/screencasts/handsonstart/

Constants

We define here the parameters of the system. The semicolons keep Mathematica from printing output. Several definitions may be combined in a single cell, in which case all will be evaluated when you press Shift-Enter.

g = 9.8;     (* acceleration of gravity, m/(s^2) *)
M = 0.375;     (* mass of the disk, kg *)
R = 0.125;     (* radius of the disk, m *)
m = 0.01;     (* mass of the pin, kg *)

Equations of motion

Define the equations of motion as a list of equations (surrounded by curly braces). The biggest gotcha to avoid here is using a single "=" (assign to) instead of a double "==" (is equal to). Unfortunately, this is a huge mistake, because it can be difficult to get Mathematica to forget the definition. Try Remove[ θ ] to completely reset any information about the symbol θ. If that doesn’t work, you’ll have to quit Mathematica and start again to remove the unwanted assignment. You can type θ by pressing Esc-q-Esc. Other Greek letters make more sense: Esc-s-Esc gives you σ, Esc-t-Esc gives τ, etc. You can also use the Palettes menu to display a GUI for entering symbols.

EoM = { (1/2)M (R^2) θ''[t]== - g m R Sin[θ[t]] - β θ'[t],
θ'[0] == 0,
θ[0] == θ0};

Note that the prime (apostrophe) indicates differentiation with respect to the function argument. θ'[t] may also be written D[θ[t], t] or Derivative[1][θ][t], but the prime notation is the easiest to understand for functions that depend on a single variable.

Solving the equations of motion

We solve the equations of motion numerically for a specified period of time. Note that there are two as-yet undefined parameters: the damping parameter β and the initial angular position θ0. We will supply these as arguments to a function that runs the solver to produce an interpolating polynomial, which represents the solution. The underscores after function arguments (such as b_) are crucial to tell Mathematica that these symbols are function arguments, not pre-defined parameters. The expression  /. {β -> b, θ0 -> q0} means to substitute for β the parameter value b, and θ0 gets the value of q0.

solv[b_, q0_, duration_] := First[ NDSolve[ EoM /. {β->b,θ0->q0}, θ,{t, 0, duration}]]

Let's see what happens when we call this function:

s1 = solv[ 0.001, π/4, 10 ]

Not particularly informative, but this interpolating function, which we have named s1, makes it easy to evaluate the solution at any time between 0 and 10. We just need to apply s1:

θ[2.3] /. s1

We can also make a plot

Plot[ Evaluate[ θ[t] /. s1],{t,0,10},{PlotRange->All, Frame->True,LabelStyle->Larger}]

We can even make a phase plot. We'll color the line from red through purple to help show the progression in the plot.

ParametricPlot[ Evaluate[ { θ'[t], θ[t]} /. s1], {t,0,10}, Frame->True, LabelStyle->Larger,FrameLabel->{θ[t],θ'[t]},
ColorFunction->Function[{x,y,t},Hue[t]],
PlotStyle->Thick,
AspectRatio->1 ]

Exploring the solution using Manipulate[]

To explore ranges of parameters, you can create an interactive panel with which to vary parameters to observe the changes in a plot nearly in realtime.

Manipulate[
ParametricPlot[ Evaluate[ {θ'[t],θ[t]} /. solv[β, θ0, 100]],{t,0,100},
Frame->True, AspectRatio->1,PlotStyle->Thick,
ColorFunction-> Function[{x,y,t},Hue[t]],
LabelStyle->Larger,FrameLabel->{θ[t],θ'[t]},
PlotRange->All],
{β,0,0.001},
{θ0, 0.1, π}]