Preface

This tutorial was made solely for the purpose of education and it was designed for students taking Applied Math 0330. It is primarily for students who have very little experience or have never used Mathematica before and would like to learn more of the basics for this computer algebra system. As a friendly reminder, don't forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to the main page for the course APMA0340
Return to the main page for the course APMA0330
Return to Part IV of the course APMA0330

Pendulum

 

I. Pendulum Equation

Christiaan Huygens.

There seems no doubt that the first person who investigated and established the mathematical theory and properties of the pendulum was a prominent Dutch mathematician and scientist Christiaan Huygens (this spelling of his name is taken from the title of his 1658 book Horologium Oscillatorium). The main point of Huygens (1629--1695) discovery was that the curve in which a bob hang by a string of insensible weight must move in order to be isochronous in its vibrations, is not a circle, but a cycloid.

The pendulum consists of a particle (or bob) of mass m attached to the end of a light inextensible rod of length \( \ell \) with the motion taking place in a vertical plane. Such pendulum is usually referred to as a simple pendulum. Let θ be the angular coordinate of m measured counterclockwise from the down position. The kinetic and potential energies are

\[ \mbox{K} = \frac{1}{2}\, I_0 \dot{\theta}^2 , \qquad\mbox{and} \qquad \Pi = mgy = mg\ell \left( 1- \cos \theta \right) , \]
where \( I_0 = m\ell^2 \) is the moment of inertia of the pendulum. Calculating partial derivatives, we get
\[ \frac{\partial \mbox{K}}{\partial \dot{\theta}} = I_0 \dot{\theta} = m\ell^2 \dot{\theta} , \qquad \frac{\partial \Pi}{\partial \theta} = mgy = mg\ell \, \sin \theta . \]
Using these expressions, we obtain from the Euler--Lagrange equation \( \displaystyle \frac{\text d}{{\text d} t} \left( \frac{\partial {\cal L}}{\partial \dot{\theta}} \right) = \frac{\partial {\cal L}}{\partial \theta} , \) for the Lagrangian \( {\cal L} = \mbox{K} - \Pi , \) the pendulum equation in a vacuum:
\[ \ddot{\theta} + \left( g/\ell \right) \sin \theta =0 \qquad \mbox{or} \qquad \ddot{\theta} + \omega_0^2 \sin \theta =0 , \]
where \( \ddot{\theta} = {\text d}^2 \theta / {\text d}t^2 \) , \( \omega_0 = \sqrt{g/\ell} >0 ,\) and g is gravitational acceleration. For small oscilaltions, we can replace sine function by its linear approximation, which leads to a linear differential equation
\[ \ddot{\theta} + \left( g/\ell \right) \theta =0 \qquad \Longleftrightarrow \qquad \ddot{\theta} + \omega_0^2 \theta =0. \]
We assume that the oscillations of the pendulum are subjected to the initial conditions \( \theta(0) = \theta_0 , \qquad \dot{\theta}(0) =0, \) where θ0 is the initial amplitude of the oscillation. The system oscillates between symmetric limits \( [- \theta_0 , \theta_0 ] . \) The periodic motion exhibited by a simple pendulum is harmonic only for small angle oscillations. beyond this limit, the equation of motion is nonlinear: the simple harmonic motion is unsatisfactory to model the pendulum motion for large amplitudes and in such cases the period depends on amplitude. The periodic solution \( \theta (t) \) and the angular frequency ω (also with the period \( T = 2\pi /\omega \) ) depends on the initial amplitude \( \theta_0 . \)

Simple pendulum equation \( \ddot{\theta} + \omega_0^2 \sin \theta =0 , \) although straightforward in appearance, is in fact rather difficult to solve because of the nonlinearity of the term \( \sin \theta . \) In order to obtain the exact solution, this equation is multiplied by the integrating factor \( \dot{\theta} = {\text d}\theta / {\text d} t, \) so that it becomes

\[ \dot{\theta} \,\ddot{\theta} + \dot{\theta}\, \sin\theta =0 , \]
which can be rewritten as
\[ \frac{{\text d}}{{\text d}t} \left[ \frac{1}{2} \left( \frac{{\text d} \theta}{{\text d}t} \right)^2 - \omega_0^2\, \cos\theta \right] =0 . \]
This equation, which corresponds to the conservation of the mechanical energy, is immediately integrable, taking into account the initial conditions:
\[ \left( \frac{{\text d} \theta}{{\text d}t} \right)^2 = 2\,\omega_0^2\left( \cos\theta - \cos \theta_0 \right) , \]
which can be written as follows:
\[ \left( \frac{{\text d} \theta}{{\text d}t} \right)^2 = 4\,\omega_0^2\left( \sin^2 \left( \frac{\theta_0}{2} \right) - \sin^2 \left( \frac{\theta}{2} \right) \right] , \]
where we used the trigonometric identity \( \cos \theta = 1 - 2\,\sin^2 \frac{\theta}{2} . \) Now let
\[ y = \sin \left( \frac{\theta}{2} \right) \qquad\mbox{and} \qquad k= \sin \left( \frac{\theta_0}{2} \right) . \]
Then \( y(0) = \sqrt{k} . \) It is easy to obtain \( {\text d}\theta / {\text d}t \) as a function of \( {\text d}y / {\text d}t \) as follows:
\[ \frac{{\text d}y}{{\text d}t} = \frac{{\text d}y}{{\text d}\theta} \, \frac{{\text d}\theta}{{\text d}t} = \frac{1}{2}\, \frac{{\text d}\theta}{{\text d}t} \,\cos \left( \frac{\theta}{2} \right) . \]
Then square it, we have
\[ \left( \frac{{\text d}y}{{\text d}t} \right)^2 = \frac{1}{4} \, \cos^2 \left( \frac{\theta}{2} \right) \left( \frac{{\text d}\theta}{{\text d}t} \right)^2 = \frac{1}{4} \left[ 1 - \sin^2 \left( \frac{\theta}{2} \right) \right] \left( \frac{{\text d}\theta}{{\text d}t} \right)^2 = \frac{1}{4} \left( 1 - y^2 \right) \left( \frac{{\text d}\theta}{{\text d}t} \right)^2 . \]
Its square derivative will be
\[ \left( \frac{{\text d}\theta}{{\text d}t} \right)^2 = \frac{4}{1-y^2} \left( \frac{{\text d}y}{{\text d}t} \right)^2 . \]
Substituting this relation into the integrated pendulum equation, we obtain
\[ \frac{4}{1-y^2} \left( \frac{{\text d}y}{{\text d}t} \right)^2 = 4\,\omega_0^2 \left( k - y^2 \right) , \]
which can be rewritten as
\[ \left( \frac{{\text d}y}{{\text d}t} \right)^2 = \omega_0^2 k \left( 1 - y^2 \right) \left( 1 - \frac{y^2}{k} \right) . \]
To simplify the above expression, we define new variables
\[ \tau = \omega_0 t \qquad\mbox{and}\qquad z = \frac{y}{\sqrt{k}} . \]
This leads to
\[ \left( \frac{{\text d}z}{{\text d}\tau} \right)^2 = \left( 1 - z^2 \right) \left( 1 - k\,z^2 \right) , \qquad z(0) = 1, \quad \left. \frac{{\text d}z}{{\text d}\tau} \right\vert_{\tau =0} =0. \]
Separation of variables yields
\[ {\text d}\tau = \pm \frac{{\text d}z}{\sqrt{\left( 1 - z^2 \right) \left( 1 - k\,z^2 \right)}} . \]
The time τ to go from point (1,0) to point \( \left( z, {\text d}z / {\text d} \tau \right) \) in the lower half-plane of the graph of \( {\text d}z / {\text d} \tau \) as a function of z is
\[ \tau = -\int_1^z \frac{{\text d}\zeta}{\sqrt{\left( 1 - \zeta^2 \right) \left( 1 - k\,\zeta^2 \right)}} . \]
We rewrite the integral in the right hand side as
\[ \tau = \int_0^1 \frac{{\text d}\zeta}{\sqrt{\left( 1 - \zeta^2 \right) \left( 1 - k\,\zeta^2 \right)}} - \int_0^z \frac{{\text d}\zeta}{\sqrt{\left( 1 - \zeta^2 \right) \left( 1 - k\,\zeta^2 \right)}} , \]
which allows us to obtain τ as a function of z and k as
\[ \tau (z) = K(k) - F(\arcsin z, k ) , \]
where K(m) and F(φ , m) are the complete and the incomplete elliptic integral of the first kind defined as follows
\begin{eqnarray*} K(m) &=& \int_0^1 \frac{{\text d}\zeta}{\sqrt{\left( 1 - \zeta^2 \right) \left( 1 - m\,\zeta^2 \right)}} \\ F(\phi , m) &=& \int_0^{\phi} \frac{{\text d}\zeta}{\sqrt{\left( 1 - \zeta^2 \right) \left( 1 - m\,\zeta^2 \right)}} , \qquad z= \sin \phi . \end{eqnarray*}
The period of oscillation T is four times the time taken by the pendulum to swing from \( \theta =0 \ (z=0) \) to \( \theta = \theta_0 \ (z=1) . \) Therefore,
\[ T = 4\,t(0) = \frac{4\,\tau (0)}{\omega_0} = \frac{4}{\omega_0} \,K(k) = \frac{2}{\pi}\, T_0 \,K(k) = \sqrt{\frac{\ell}{g}} \,\int_0^{\pi /2} \frac{{\text d} \beta}{\sqrt{1- k^2 \sin^2 \beta}} , \]
where
\[ T_0 = \frac{2\pi}{\omega_0} = 2\pi \,\sqrt{\frac{\ell}{g}} \]
is the period of the pendulum for small oscillations. The solution formula can be rewritten as
\[ F \left( \arcsin z , k \right) = K(k) - \tau , \]
which can be written in terms of the Jacobi elliptic function sn(u;m):
\[ z = \mbox{sn} \left( K(k) - \tau ; k \right) , \]
where Mathematica has special commands for their evaluations: EllipticK[m] for complete elliptic integral of the first kind, and EllipticF[φ,m] for incomplete elliptic integral of the first kind. Mathematica also uses a special command JacobiSN[u,m] for the Jacobi elliptic function sn(u;m). These functions are discussed in the second part of the course.

Now we return to the original variables
\[ \sin \left( \frac{\theta}{2} \right) = \sin \frac{\theta_0}{2} \, \mbox{sn} \left[ K\left( \sin^2 \frac{\theta_0}{2} \right) - \omega_0 t ; \sin^2 \frac{\theta_0}{2} \right] , \]
which allows us to express θ as a function of t:
\[ \theta (t) = 2\,\arcsin \left\{ \sin\frac{\theta_0}{2} \, \mbox{sn} \left[ K\left( \sin^2 \frac{\theta_0}{2} \right) - \omega_0 t ; \sin^2 \frac{\theta_0}{2} \right] \right\} . \]
Using the relations \( \dot{\bf e}_{r} = \dot{\theta} {\bf e}_{\theta} \) and \( \dot{\bf e}_{\theta} = -\dot{\theta} {\bf e}_{r} , \) we obtain the velocity, acceleration, and jerk:
\begin{eqnarray*} {\bf v} (t) &=& \ell \dot{\theta} {\bf e}_{\theta} = \sqrt{2\ell g} \left( \cos \theta - \cos \theta_0 \right)^{1/2} {\bf e}_{\theta} , \\ {\bf a} (t) &=& \dot{\bf v} (t) = -\ell \dot{\theta}^2 {\bf e}_{r} + \ell \ddot{\theta} {\bf e}_{\theta} = - 2g \left( \cos \theta - \cos \theta_0 \right) {\bf e}_{r} - g \,\sin \theta \, {\bf e}_{\theta} , \\ {\bf j} (t) &=& \dot{\bf a} (t) = 3g\,\sin\theta \, \dot{\theta} \, {\bf e}_r - \dot{\theta} \left[ 2\ell \dot{\theta}^2 + g\,\sin \theta \right] {\bf e}_{\theta} . \end{eqnarray*}
Manipulate[
Module[{run, sol, z, app, a},
If[run == 10, run = 0];
sol = DSolve[{y''[x] + (9.8/l)*y[x] == 0, y[0] == int, y'[0] == 0}, y, x];
z[x_] := y[x] /. sol;
app = NDSolve[{\[Theta]1''[t] + (9.8/l)*Sin[\[Theta]1[t]] == 0, \[Theta]1[0] == int, \[Theta]1'[0] == 0}, \[Theta]1, {t, 0, 10}];
a[t_] := \[Theta]1[t] /. First[app];
Grid[{{
Graphics[{
{Thick, Line[{{0, 0}, {l*Sin[a[r]], -l*Cos[a[r]]}}]}, {Darker[Green], Disk[{0, 0}, .03]}, {Purple,
Disk[{l*Sin[a[r]], -l*Cos[a[r]]}, 0.1 Power[(30 m)/(4 \[Pi]), (3)^-1]]}},
PlotRange -> {{-3.4, 3.4}, {.5, -3.4}},
ImageSize -> {250, 250}],
Plot[{a[t], z[t]}, {t, 0, 10},
PlotRange -> {{0, 10}, {-1.6, 1.6}},
PlotStyle -> {{Thick, Blue}, {Red}}, ImageSize -> {250, 250},
AspectRatio -> 1.2, Frame -> True,
PlotLegends`PlotLegend -> {"numerical", "sin(\[Theta]) = \[Theta]"},
FrameLabel -> {"time", "displacement"},
Epilog -> {Blue, PointSize[.06], Point[{r, a[r]}]}]}}]],
{{int, \[Pi]/3., "initial angle \[Theta](0)"}, \[Pi]/16., \[Pi]/2., Appearance -> "Labeled"},
Delimiter,
{{m, 4, "mass m"}, 1, 10, .01, Appearance -> "Labeled"},
{{l, 2.25, "length L"}, 1, 3, .01,
Appearance -> "Labeled"}, {{r, 0, "release system"}, 0, 10, .01,
ControlType -> Trigger, AnimationRate -> .5},
TrackedSymbols :> {int, m, l, r},
Initialization :> Get["PlotLegends`"]]

 

III. Pendulum Equation with Resistance

The resistance of the air acts on both the pendulum ball and the pendulum wire. It causes the amplitude to decrease with time and increases the period of oscillation slightly. The drag force is proportional to the velocity for values of the Reynolds number of the order 1 or less. For values of the Reynolds number of order \( 10^3 \) to \( 10^5 , \) the force is proportional to the square of the velocity. The maximum Reynolds number based on diameter for the ball is 1100, where the quadratic force law should apply, while the maximum value based on diameter of wire is about 6, where the linear force law should prevail.

Since the damping force is neither linear nor quadratic, but rather a combination of the two, it makes sense to establish a damping function which contains both effects simultaneously. Therefore, the pendulum equation becomes

\[ \ddot{\theta} + \alpha \,\dot{\theta} + \epsilon \,\mbox{sign} \left( \dot{\theta} \right) + \beta \left\vert \dot{\theta} \right\vert \dot{\theta} + \omega_0^2 \sin\theta =0 , \]

where the linear damping coefficient \( \alpha = \kappa /m , \) describes air resistance due to wire, the quadratic damping coefficient β is due to air drag on the pendulum bob, the third damping term \( \quad \epsilon = c\ell /(2m) \) is necessary for taking into count the bearing, which is subject to dry friction (Coulomb damping); here \( \quad \omega_0 = \left( g/\ell \right)^{1/2} ,\) and \( \mbox{sign}(x) = \begin{cases} 1, & \ x>0, \\ -1, & \ x<0. \end{cases} \)

Analysis of the experimental data gives the following values of the parameters:

\begin{eqnarray*} \alpha &=& 0.012 \pm 0.005 \,\mbox{s}^{-1} , \\ \beta &=& 0.0291 \pm 0.0064 , \\ \epsilon &=& 0.0017 \pm 0.0009 \,\mbox{s}^{-1} . \end{eqnarray*}

It is significant that the effect of air drag is apparent in both the parameters α and β. If air drag could be described purely in terms of the velocity squared, it should affect only the parameter β. If one attemps to fit the data with a damping law based on a single power of the velocity, the power would have to be somewhat less than 2.

In order to find out the relationship between the angle of the pendulum and time t, a good way is to plot θ against t. In order to make the plot, first solve the differential equation with the command NDSolve.

This can be done by assigning constant values to α, ε and ω when using the NDSolve command, get three different functions in respect to three different α values, and use the Plot command to plot the three functions in one graph

numsol1 =
  NDSolve[{y''[x] == -10*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsol2 =
  NDSolve[{y''[x] == -1*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsol3 =
  NDSolve[{y''[x] == -0.1*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
Plot[Evaluate[y[x] /. {numsol1, numsol2, numsol3}], {x, 0, 100}, PlotRange -> All, PlotStyle -> {GrayLevel[0.1], Dashed, Thick}, AxesLabel -> {"t", "\[Theta]"}, PlotLegends -> {"\[Alpha]=10", "\[Alpha]=1", "\[Alpha]=0.1"}]

Another method is to set up a function that takes θ0, ω0, α, ε and tmax as input and out put a function of θ

DampedDrivenPendulumFunction2[{\[Theta]0_, \[Omega]0_}, \[Alpha]_, \ \[Epsilon]_, tmax_] :=
Module[{\[Theta]}, \[Theta] /. First[ NDSolve[{ \[Theta]''[t] == -Sin[\[Theta][t]] - 2 \[Alpha] \[Theta]'[ t] - \[Epsilon] Sign[\[Theta]'[t]] (\[Theta]'[t])^2, \[Theta][0] == \[Theta]0, \[Theta]'[0] == \[Omega]0}, \[Theta], {t, 0, tmax}
Plot[{DampedDrivenPendulumFunction2[{0.7, 0}, 0.1, 0, 100][x], DampedDrivenPendulumFunction2[{0.7, 0}, 1, 0, 100][x], DampedDrivenPendulumFunction2[{0.7, 0}, 10, 0, 100][x]}, {x, 0, 100}, PlotLegends -> "Expressions", PlotRange -> Full, PlotStyle -> {Thick, Automatic, Dashed}, AspectRatio -> 1.0]
The third method presented here is very different from the first and the second: we create a loop to add solutions into a list.
\[Epsilon] = 0; g = 9.8; \[Omega]0 = 1;
listofsolutions = {};
listofvaluesforalpha = {.1, 1, 10};
Numerically solve the corresponding second order differential equation, leaving alpha as a variable:
pendulumeq = \[Theta]''[t] + 2*\[Alpha] \[Theta]'[ t] + \[Epsilon] Sign[\[Theta]'[ t]]*((\[Theta]'[t])^2) + (\[Omega]0^2)*Sin[\[Theta][t]] == 0;
Get distinct equations for the different values of \[Alpha] and add them to the listofsolutions:
counter = 1;
While[counter <= Length[listofvaluesforalpha], AppendTo[listofsolutions, NDSolve[{Evaluate[ pendulumeq /. {\[Alpha] -> listofvaluesforalpha[[counter]]}], \[Theta][ 0] == .7, \[Theta]'[0] == 0}, \[Theta][t], {t, 0, 100}]]; counter++]
Plot[Evaluate[theta[t] /. listofsolutions], {t, 0, 50}, PlotRange -> All]

The three different methods all get the same result, from which we can see that the largest damping coefficient does not stop the pendulum the fastest. So the qualitative difference among the three different coefficients is that when α is 0.1, the pendulum oscillates (cross 0). In other cases, the pendulum doesn't oscillate (doesn't cross 0).

The critical damping value is the damping coefficient that stop the pendulum the fastest.

In order to find out the critical damping value of alpha, first we need to define the meaning of stop, because the pendulum will end up always be in very small oscillation.

The first method defined - 0.035 < θ < 0.035 as stop. By plotting different functions of θ with different alpha values (0.6, 0.7, 0.8, 0.9, 1.0) in one graph, you can see that α = 0.7 starts to be in the range [-0.035, 0.035] the fastest.

numsola1 =
  NDSolve[{y''[x] == -0.6*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolb1 =
  NDSolve[{y''[x] == -0.7*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolc1 =
  NDSolve[{y''[x] == -0.8*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsold1 =
  NDSolve[{y''[x] == -0.9*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolf1 =
  NDSolve[{y''[x] == -1.0*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
Plot[Evaluate[ y[x] /. {numsola1, numsolb1, numsolc1, numsold1, numsolf1}], {x, 3, 9}, PlotStyle -> {GrayLevel[0.2], Dashed, Thick}, PlotRange -> All, AxesLabel -> {"t", "\[Theta]"}, PlotLegends -> {"\[Alpha]=0.6", "\[Alpha]=0.7", "\[Alpha]=0.8", "\[Alpha]=0.9", "\[Alpha]=1.0"}]
To find out more accurate value of α, plot several different α values close to 0.7:
numsola =
  NDSolve[{y''[x] == -0.74*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolb =
  NDSolve[{y''[x] == -0.73*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolc =
  NDSolve[{y''[x] == -0.72*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsold =
  NDSolve[{y''[x] == -0.71*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolf =
  NDSolve[{y''[x] == -0.70*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolg =
  NDSolve[{y''[x] == -0.69*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolh =
  NDSolve[{y''[x] == -0.68*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsoli =
  NDSolve[{y''[x] == -0.67*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
Plot[Evaluate[ y[x] /. {numsola, numsolb, numsolc, numsold, numsolf, numsolg, numsolh, numsoli}], {x, 3, 8}, PlotStyle -> {GrayLevel[0.1], Dashed, Thick}, PlotRange -> All, AxesLabel -> {"t", "\[Theta]"}, PlotLegends -> {"\[Alpha]=0.74", "\[Alpha]=0.73", "\[Alpha]=0.72", "\[Alpha]=0.71", "\[Alpha]=0.70", "\[Alpha]=0.69", "\[Alpha]=0.68", "\[Alpha]=0.67"}]
Then α = 0.69 should roughly be the critical damping value.
Pi/2*0.005
Out[25]= 0.00785398
findmid[x_, y_] := (x + y)/2
Plot[{DampedDrivenPendulumFunction2[{0.7, 0}, 1, 0, 100][x], DampedDrivenPendulumFunction2[{0.7, 0}, 0.7, 0, 100][x]}, {x, 0, 30}, PlotLegends -> "Expressions", PlotRange -> {-0.007853981633974483`, 0.007853981633974483`}, PlotStyle -> {Thick, Automatic, Dashed}, AspectRatio -> 1.0]
From this graph, when α = 1, the pendulum settles in the stop range faster than when α = 0.7. This means that 0.7 is too small, and 1 may be too large. Then try the mid value 0.75.
findmid[1, 0.7]
Plot[{DampedDrivenPendulumFunction2[{0.7, 0}, 0.85, 0, 100][x], DampedDrivenPendulumFunction2[{0.7, 0}, 1, 2, 100][x], DampedDrivenPendulumFunction2[{0.7, 0}, 0.7, 0, 100][x]}, {x, 0, 30}, PlotLegends -> "Expressions", PlotRange -> {-0.007853981633974483`, 0.007853981633974483`}, PlotStyle -> {Thick, Automatic, Dashed}, AspectRatio -> 1.0]

When α = 0.85, the pendulum settles in the range faster than when α = 1. So the critical value may be smaller than 0.85. Try the mid value of 0.85 and 0.7. Finally, the answer got from this method is 0.81.

The third method defines stop as not passing 0 in the time span 0 < t \ < 500, uses a function "piano" to check if the pendulum oscillates.

\[CurlyEpsilon] = 0; g = 9.8;
\[CapitalOmega]0 = 1;
listofsolutions = {};
listofvaluesfor\[Alpha] = Range[0, 2, .0005];
pendulumeq = \[Theta]''[t] + 2*\[Alpha]*\[Theta]'[t] + \[CurlyEpsilon]* Sign[\[Theta]'[t]]*((\[Theta]'[t])^2) + (\[CapitalOmega]0^2)* Sin[\[Theta][t]] == 0;
counter = 1;
While[counter <= Length[listofvaluesfor\[Alpha]], AppendTo[listofsolutions, NDSolve[{Evaluate[ pendulumeq /. {\[Alpha] -> listofvaluesfor\[Alpha][[counter]]}], \[Theta][ 0] == .7, \[Theta]'[0] == 0}, \[Theta][t], {t, 0, 1000}]]; counter++];
evaluatedlistofsolutions = Evaluate[\[Theta][t] /. listofsolutions];
Piano is a function that checks to see if the pendulum ever swings past θ = 0 for each of the functions in the list. If the answer is no, then it means that the pendulum no longer oscillates for that value of α.
piano[key_] :=
  Catch[If[FindMinValue[ evaluatedlistofsolutions[[key]], {t, 0, 500}] > 0, Throw[key], Throw[906]]];
counter2 = 1;
While[piano[counter2] == 906, counter2++];
Print["The Value of \[Alpha] at which the function \[Theta][t] stops \ oscillating is ", listofvaluesfor\[Alpha][[counter2]]]
The Value of α at which the function θ[t] stops oscillating is 0.8985
Plot[evaluatedlistofsolutions[[piano[counter2]]], {t, 0, 20}, PlotRange -> All]
Now we assume that alpha is not zero. .
DampedDrivenPendulumFunction2[{\[Theta]0_, \[Omega]0_}, \[Alpha]_, \ \[Epsilon]_, tmax_] :=
  Module[{\[Theta]}, \[Theta] /. First[ NDSolve[{ \[Theta]''[t] == -Sin[\[Theta][t]] - 2 \[Alpha] \[Theta]'[ t] - \[Epsilon] Sign[\[Theta]'[t]] (\[Theta]'[t])^2, \[Theta][0] == \[Theta]0, \[Theta]'[0] == \[Omega]0}, \[Theta], {t, 0, tmax}] ]]
g5 = DampedDrivenPendulumFunction2[{0.7, 0}, 0.1, 1, 100]
g6 = DampedDrivenPendulumFunction2[{0.7, 0}, 10, 1, 100]
g7 = DampedDrivenPendulumFunction2[{0.7, 0}, 1, 1, 100]
Plot[{g5[x], g6[x], g7[x]}, {x, 0, 100}, PlotLegends -> "Expressions", PlotRange -> Full, PlotStyle -> {Thick, Automatic, Dashed}]
When ε is 1, the qualitative difference among different alpha when apha is 1, 0.1, or 10 is the same. But when alpha is 0.1, it oscillates with less times compared to the case when epsilon is 0.

 

III. Driven Pendulum

Driven pendulum may exhibit a chaotic motion. It consists of a mass m fixed at a distance \( \ell \) from a pivot which is subject to a vertical oscillation \( y = A\, \cos \left( \omega \, t \right) . \) Let θ be the angular coordinate of m measured counterclockwise from the down position, and \( \phi = \pi - \theta \) the complementary displacement measured clockwise from the up orientation. The kinetic and potential energies are

 

IV. Rocking Rigid Pendulum

plus = Polygon[{{-1/2, 0}, {-1/2, 1/4}, {-4, 1/4}, {-4, 0}}]
minus = Polygon[{{1/2, 0}, {1/2, 1/4}, {4, 1/4}, {4, 0}}]
a = Show[Graphics[{Pink, minus}], Graphics[{Pink, plus}]]
line1 = Graphics[Line[{{-2, -3.5}, {0, 0.6}}], PlotRange -> {{-2, 2}, {-4, 1}}]
line2 = Graphics[Line[{{-0.85, 1.0}, {0.8, 0.25}}], PlotRange -> {{-2, 2}, {-4, 1}}]
makeArrowPlot[g_Graphics, ah_: 0.05, dx_: 1*^-6, dy_: 1*^-6] :=
Module[{pr = PlotRange /. Options[g, PlotRange], gg, lhs, rhs},
gg = g /. GraphicsComplex -> (Normal[GraphicsComplex[##]] &);
lhs := Or @@
Flatten[{Thread[Abs[#[[1, 1, 1]] - pr[[1]]] < dx], Thread[Abs[#[[1, 1, 2]] - pr[[2]]] < dy]}] &;
rhs := Or @@ Flatten[{Thread[Abs[#[[1, -1, 1]] - pr[[1]]] < dx], Thread[Abs[#[[1, -1, 2]] - pr[[2]]] < dy]}] &;
gg = gg /.
x_Line?(lhs[#] && rhs[#] &) :> {Arrowheads[{-ah, ah}], Arrow @@ x};
gg = gg /. x_Line?lhs :> {Arrowheads[{-ah, 0}], Arrow @@ x};
gg = gg /. x_Line?rhs :> {Arrowheads[{0, ah}], Arrow @@ x};
gg]
curve = Plot[{-Sqrt[9 - x^2]}, {x, -2.2, 2.2}, PlotStyle -> {Thick, Dashed}, Axes -> False] // makeArrowPlot
Show[a, line1, line2, curve]

Figure 1: Rocking rigid pendulum.

A rigid pendulum, such as a long rod, has a short cross-bar rigidly attached to it. As it oscillates in vertical plane, it is supported on ends of the cross-bar as they alternately contact a fixed horizontal surface:
plus = Polygon[{{-1/2, 0}, {-1/2, 1/4}, {-4, 1/4}, {-4, 0}}]
minus = Polygon[{{1/2, 0}, {1/2, 1/4}, {4, 1/4}, {4, 0}}]
a = Show[Graphics[{Pink, minus}], Graphics[{Pink, plus}]]
line1 = Graphics[Line[{{-2, -3.5}, {0, 0.6}}], PlotRange -> {{-2, 2}, {-4, 1}}]
line2 = Graphics[Line[{{-0.85, 1.0}, {0.8, 0.25}}], PlotRange -> {{-2, 2}, {-4, 1}}]
p2 = Graphics[{Dashed, Arrow[{{-0.66, -0.8}, {-0.66, -3.8}}]}]
point = Graphics[{PointSize[Large], Green, Point[{-0.66, -0.8
t1 = Graphics[ Text[Style["\[Theta]", FontSize -> 14, Red], {-1.0, -2.4}]]
t2 = Graphics[Text[Style["G", FontSize -> 14, Blue], {-1.0, -0.8}]]
a1 = Graphics[Text[Style["a", FontSize -> 14, Blue], {0.5, 0.7}]]
a2 = Graphics[Text[Style["a", FontSize -> 14, Blue], {-0.2, 1.0}]]
Q = Graphics[Text[Style["Q", FontSize -> 14, Blue], {1.0, 0.45}]]
b = Graphics[Text[Style["b", FontSize -> 14, Blue], {0.0, -0.2}]]
mg = Graphics[Text[Style["mg", FontSize -> 14, Black], {-0.2, -3.6}]]
Show[a, line1, line2, point, p2, t1, t2, a1, a2, Q, b, mg]
plus = Polygon[{{-1/2, 0}, {-1/2, 1/4}, {-4, 1/4}, {-4, 0}}]
minus = Polygon[{{1/2, 0}, {1/2, 1/4}, {4, 1/4}, {4, 0}}]
a = Show[Graphics[{Pink, minus}], Graphics[{Pink, plus}]]
line1 = Graphics[Line[{{2, -3.5}, {0, 0.6}}], PlotRange -> {{-2, 2}, {-4, 1}}]
line2 = Graphics[Line[{{0.85, 1.0}, {-0.8, 0.25}}], PlotRange -> {{-2, 2}, {-4, 1}}]
p2 = Graphics[{Dashed, Arrow[{{0.66, -0.8}, {0.66, -3.8}}]}]
point = Graphics[{PointSize[Large], Green, Point[{0.66, -0.8}]}]
t1 = Graphics[ Text[Style["\[Theta]", FontSize -> 14, Red], {1.0, -2.4}]]
t2 = Graphics[Text[Style["G", FontSize -> 14, Blue], {0.3, -0.8}]]
P = Graphics[Text[Style["P", FontSize -> 14, Blue], {-0.9, 0.45}]]
b = Graphics[Text[Style["b", FontSize -> 14, Blue], {0.0, -0.2}]]
mg = Graphics[Text[Style["mg", FontSize -> 14, Black], {0.22, -3.6}]]
Show[a, line1, line2, point, p2, t1, t2, P, b, mg]
plus = Polygon[{{-1/2, 0}, {-1/2, 1/4}, {-4, 1/4}, {-4, 0}}]
minus = Polygon[{{1/2, 0}, {1/2, 1/4}, {4, 1/4}, {4, 0}}]
a = Show[Graphics[{Pink, minus}], Graphics[{Pink, plus}]]
line1 = Graphics[Line[{{0, -3.5}, {0, 0.6}}], PlotRange -> {{-2, 2}, {-4, 1}}]
line2 = Graphics[Line[{{0.85, 0.27}, {-0.85, 0.27}}], PlotRange -> {{-2, 2}, {-4, 1}}]
P = Graphics[Text[Style["P", FontSize -> 14, Blue], {0.9, 0.45}]]
Q = Graphics[Text[Style["Q", FontSize -> 14, Blue], {-0.9, 0.45}]]
mg = Graphics[Text[Style["mg", FontSize -> 14, Black], {0.4, -3.6}]]
point = Graphics[{PointSize[Large], Green, Point[{0, -0.8}]}]
t2 = Graphics[Text[Style["G", FontSize -> 14, Blue], {-0.3, -0.8}]]
p3 = Graphics[{Dashed, Line[{{0, -0.8}, {0.85, 0.27}}]}]
b = Graphics[Text[Style["b", FontSize -> 14, Blue], {-0.3, -0.2}]]
a1 = Graphics[Text[Style["a", FontSize -> 14, Blue], {0.3, 0.5}]]
a2 = Graphics[Text[Style["a", FontSize -> 14, Blue], {-0.3, 0.5}]]
ell = ToExpression["\ell", TeXForm]
t3 = Graphics[Text[Style[ell, FontSize -> 14, Black], {0.6, -0.4}]]
Show[a, line1, line2, point, p3, t2, a1, a2, P, Q, b, mg, t3]
     

Let the mass of pendulum be m, and the length of attached bar be 2 a. Since the construction of rigid pendulum is symmetric, the center of gyration, which we denote by G, is along the main rod. Let k be the radius of gyration of pendulum about G, so that the square of its radius of gyration about P and Q is \( k^2 + \ell^2 . \)

 

The first half-cycle of the motion

Suppose that the pendulum is set in motion from the central position with initial conditions

\[ \theta =0 \quad \mbox{and} \quad \dot{\theta} = \omega_1 \qquad\mbox{at} \quad t=0. \]
Since energy is assumed to be conservative, the total energy is a constant:
\[ \frac{1}{2}\, m \left( k^2 + \ell^2 \right) \dot{\theta}^2 + mg \left( a\,\sin \theta -b\,\cos \theta \right) = \mbox{constant} . \]
For small oscillations, powers of θ higher than the first may be neglected, yielding
\[ \left( k^2 + \ell^2 \right) \dot{\theta}^2 + 2ag \,\theta = \mbox{constant} \qquad\mbox{or} \qquad \dot{\theta}^2 + c^2 \theta = \mbox{constant} , \]
where \( c^2 = 2ga \left( k^2 + \ell^2 \right) . \) The initial conditions enable this to be written as
\[ \dot{\theta}^2 + c^2 \theta = \omega_1^2 , \]
When the pendulum first comes to rest, let
\[ t= \frac{1}{2}\,\tau_1 , \quad \theta = \theta_1 , \quad \dot{\theta} =0. \]
Then
\[ \dot{\theta}^2 + c^2 \theta = \omega_1^2 \qquad\Longrightarrow \qquad \theta_1 = \omega_1^2 / c^2 ., \]
Now
\begin{align*} \frac{1}{2}\,\tau_1 &= \int_0^{\theta_1} \, \frac{{\text d} \theta}{\dot{\theta}} = \int_0^{\theta_1} \, \frac{{\text d} \theta}{\sqrt{\omega_1^2 - c^2 \theta}} = \frac{2}{c^2} \left[ \sqrt{\omega_1^2 - c^2 \theta} \right]_{0}^{\theta_1} \\ &= \frac{2\omega_1}{c^2} = \frac{2\sqrt{\theta_1}}{c} . \end{align*}
Clearly, the time taken for θ to reach zero again will also be \( \frac{1}{2}\,\tau_1 , \) and the angular speed will then be ω1 again. Thus, the first "half period" of oscillation will be \( \tau_1 = 4 \sqrt{\theta_1} /c . \)

 

The first transition through the central position

 

 

 

Return to Mathematica page

Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)