Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part I of the course APMA0340
Introduction to Linear Algebra with Mathematica

Inverted Pendulum


An inverted pendulum is not stable. Nevertheless, as pointed out by Pyotr Kapitza (1894--1984) in 1951, if you oscillate vertically the suspension point at a sufficientlyfast frequency Ω, you can stabilise it. A similar phenomenon occurs in a rotating saddle (see video), and is also used to create radio-frequency traps.

Let us start considering a very familiar one-dimensional system: a planar pendulum made of a massless rod of length ℓ ending with a point mass m. In the familiar swing, the driving occurs in different ways: if you “drive” the swing yourself, you do it by effectively modifying the position of your “center-of-mass”, hence the effective length ℓ(t) of the “pendulum”. If you are pushed by someone else, then you have a pendulum with a “periodic external force”. We use the generalize coordinate q = θ that denotes the angle formed with the vertical (θ = 0 being the downward position), and y0(t) denotes the position of its suspension point, we can derive the equations of motion from the Lagrangian formalism. In a short while we will assume that y0(t) = Acos Ωt, where A is the amplitude of the driving and Ω the driving frequency, but for the time being, let us proceed by keeping y0(t to be general. In a system of reference with the y-axis oriented upwards and the x-axis horizontally, the position x(t) and y(t) of the mass m is:

\[ \begin{cases} x(t) &= \ell \,\sin\theta (t) , \\ y(t) &= y_0 (t) - \ell\,\cos \theta (t) , \end{cases} \qquad \Longrightarrow \qquad \begin{cases} \dot{x}(t) &= \ell \dot{\theta} \cos\theta , \\ \dot{y}(t) &= \dot{y}_0 + \ell \dot{\theta} \sin\theta . \end{cases} \]
The Lagrangean is given by
\begin{align*} {\cal L}\left( \theta , \dot{\theta} , t \right) &= \frac{m}{2} \left( \dot{x}^2 + \dot{y}^2 \right) - mgy \\ &= \frac{m}{2} \,\ell^2 \dot{\theta}^2 + m\ell \dot{y}_0 \dot{\theta} \sin\theta + mg\ell\,\cos\theta , \end{align*}
where we dropped the term \( \frac{m}{2}\,\dot{y}_0^2 - mgy_0 \) because it would not enter in the Euler--Lagrange equations due to pure dependence on time t. The associated momentum is given by:
\[ p_{\theta} = \frac{\partial {\cal L}}{\partial \theta} = m\ell^2 \dot{\theta} + m\ell \dot{y}_0 \sin\theta \qquad \Longrightarrow \qquad \dot{\theta} = \frac{p_{\theta}}{m\ell^2} - \frac{\dot{y}_0}{\ell} \,\sin\theta . \]
The Hamilton’s equations read:
\[ \begin{cases} \dot{\theta} &= \frac{\partial H}{\partial p_{\theta}} = \frac{p_{\theta}}{m\ell^2} , \\ \dot{p}_{\theta} &= - \frac{\partial H}{\partial \theta} = -m \left[ g - A \Omega^2 \cos\Omega t \right] \sin \theta . \end{cases} \]
Hence, transforming it into a second-order equation:
\begin{equation} \label{EqI.1} \frac{{\text d}^2 \theta}{{\text d} t^2} + \frac{1}{\ell} \left[ g - A\Omega^2 \cos (\Omega t) \right] \sin \theta = 0 , \end{equation}
where \( \omega_0^2 = g/\ell \) is the square of the frequency of the unperturbed pendulum in the linear regime, Ω is the driving frequency, A is the amplitude of the driving, which we model as y0(t) = Acos(Ωt).

Example 1: If a rod or a pencil is held upright on a table, and then released, it will fall onto the ground. We consider a simple case when a rigid rod of length 2ℓ is fixed to a horizontal table by a smooth hinge at one end while another end is free to fall.
      Falling pencil
rod = Graphics[{LightGray, Polygon[{{0.1, 0}, {0, 0.1}, {2.0, 2.1}, {2.1, 2.0}}]}];
ar = Graphics[{Black, Thickness[0.01], Arrowheads[0.08], Arrow[{{-0.2, 0}, {2.5, 0}}]}];
ar2 = Graphics[{Black, Dashed, Thickness[0.01], Arrowheads[0.08], Arrow[{{0, -0.2}, {0, 2.2}}]}];
ar3 = Graphics[{Black, Dashed, Thickness[0.01], Arrowheads[0.08], Arrow[{{1, 1}, {1, 0.4}}]}];
tl = Graphics[{Black, Text[Style["\[ScriptL]", 18, FontFamily -> "Mathematica1"], {1.5, 1.75}]}];
tx = Graphics[{Black, Text[Style["x", 18], {2.4, 0.2}]}];
t1 = Graphics[{Black, Text[Style["C", 18], {1.0, 1.2}]}];
t2 = Graphics[{Black, Text[Style["\[Theta]", 18], {0.15, 0.4}]}];
t3 = Graphics[{Black, Text[Style["mg", 18], {1.05, 0.3}]}];
Show[rod, ar, ar2, ar3, tx, tl, t1, t2, t3]
       Falling pencil.            Mathematica code

Since the bottom end is hinged, there is no horizontal move of it. The kinetic energy, K, consists of the kinetic energy of the motion of the center of mass plus the kinetic energy of rotation about the bottom end of the pencil:

\[ \mbox{K} = \frac{m}{2}\,v_{cm}^2 + \frac{1}{2}\, I_0 \dot{\theta}^2 \]

Let us consider a particular example of inverted pendulum that initially was displaced by 1° ·

\[ \ddot{\theta} = 100\,\sin (\theta ) , \qquad \theta (0) = 0.0174533, \quad \dot{\theta} (0) = 0.2. \]
We solve this problem numerically:
sol = NDSolve[{y''[t] == 100*Sin[y[t]], y[0] == Pi/180, y'[0] == 9.2}, y, {t, 0, 20}];
Plot[Evaluate[y[t] /. sol], {t, 0, 2}, PlotStyle -> Thickness[0.01]]
Then we plot θ(t) for three different initial velocities.
           
Initial velocity is 0.2;       Initial velocity is 1.2;       Initial velocity is 9,2.
      For a curious reader, we present the graph of the following initial value problem (IVP):
\[ \ddot{\theta} = 100\,\cos (\theta ) , \qquad \theta (0) = 0.0174533, \quad \dot{\theta} (0) = 0.2. \]
sol2 = NDSolve[{y''[t] == 100*Cos[y[t]], y[0] == Pi/180, y'[0] == 9.2}, y, {t, 0, 20}];
Plot[Evaluate[y[t] /. sol2], {t, 0, 2}, PlotStyle -> Thickness[0.01]]
       Solution of the initial value problem d²θ/dt² = 100 cos(θ).            Mathematica code

   ■

Example 2: Let us consider a stick or pencil or any elongated object of length ℓ is at at rest in upright position standing on frictionless table. Suppose that initial at time t = 0 it was pushed in the horizontal direction either by slight move or initial velocity or both. This will cause the pencil to fall, with bottom end sliding along the ground surface. Since there is no horizontal force acting on the object, its center of mass G moves only vertically. One generalized coordinate is enough to determine the state of the system. Let it be θ, the angle between the vertical axis and the pencil.

      Falling pencil on ficktionless table:
pencil = Graphics[{Yellow, Polygon[{{0.3, 0.2}, {0, 0}, {0.2, 0.3}, {2.0, 2.1}, {2.1, 2.0}}]}]; ar = Graphics[{Black, Thickness[0.01], Arrowheads[0.08], Arrow[{{-0.2, 0}, {2.5, 0}}]}]; ar2 = Graphics[{Black, Dashed, Thickness[0.01], Arrowheads[0.08], Arrow[{{0, -0.2}, {0, 2.2}}]}]; ar3 = Graphics[{Black, Dashed, Thickness[0.01], Arrowheads[0.08], Arrow[{{1.4, 1.4}, {1.4, 0.8}}]}]; ar4 = Graphics[{Black, Thickness[0.01], Arrowheads[0.08], Arrow[{{0, 0}, {0, 1.1}}]}]; tl = Graphics[{Black, Text[Style["\[ScriptL]", 18, FontFamily -> "Mathematica1"], {1.8, 2.05}]}]; tx = Graphics[{Black, Text[Style["x", 18], {2.4, 0.2}]}]; t1 = Graphics[{Black, Text[Style["G", 18], {1.4, 1.5}]}]; t2 = Graphics[{Black, Text[Style["\[Theta]", 18], {0.15, 0.4}]}]; t3 = Graphics[{Black, Text[Style["mg", 18], {1.45, 0.5}]}]; t4 = Graphics[{Black, Text[Style["R", 18], {0.75, 0.85}]}]; t5 = Graphics[{Black, Text[Style["C", 18], {-0.1, 1.4}]}]; t6 = Graphics[{Black, Text[Style["N", Bold, 18], {0.15, 1.1}]}]; t7 = Graphics[{Black, Text[Style["O", 18], {-0.13, -0.13}]}]; t8 = Graphics[{Black, Text[Style["y", 18], {0.13, 2.13}]}]; t9 = Graphics[{Black, Text[Style["table", 18], {1.1, -0.15}]}]; line = Graphics[{Black, Dashed, Thick, Line[{{1.4, 1.4}, {0, 1.4}}]}]; Show[pencil, ar, ar2, ar3, ar4, tx, tl, t1, t2, t3, t4, t5, t6, t7, t8, t9, line]
       Fig. 2.1:     Falling pencil on fricktionless table.            Mathematica code

To derive the equation of motion, we will denote by ω = dθ/dt the derivative of θ(t) with respect to time, and by α = dω/dt its acceleration. Let O denote a reference point of contact of the rod with the table, which moves in the horizontal direction. So its vector ro has only one coordinate in the x-direction; the vector rg describing the center of mass has only one coordinate in the y-direction. Let Ao denote the angular momentum of the system with respect to O as follows:

\[ {\bf A}_o = \sum m_i \left( {\bf r}_i - {\bf r}_o \right) \times \left( \dot{\bf r}_i - \dot{\bf r}_o \right) . \]
Let τo be the sum of the external torques:
\[ {\bf \tau}_o = \sum m_i \left( {\bf r}_i - {\bf r}_o \right) \times {\bf F}_i , \]
where Fi is the external force on the i-th particle. The torque equation is
\[ \frac{\text d}{{\text d}t} \left( {\bf A}_0 \right) + M \left( {\bf r}_g - {\bf r}_0 \right) \times \ddot{\bf }_o = {\bf \tau}_o , \tag{2.1} \]
whereM is the total mass of the system. For our problem of falling pencil, we have
\[ \dot{\bf r}_o = - \frac{\text d}{{\text d}t} \left( R\,\sin \theta \right) {\bf i} = - \left( R\,\cos\theta \,\dot{\theta} \right) {\bf i} , \]
\[ \ddot{\bf r}_o = R \left( \dot{\theta}^2 \sin\theta - \ddot{\theta}\,\cos\theta \right) {\bf i} , \]
and
\[ {\bf r}_g - {\bf r}_o = R\,\sin\theta\,{\bf i} + R\,\cos\theta\,{\bf j} \]
Since the motion is restricted to the xy-plane, we use the z-component of Eq.(1):
\[ {\bf k}: \ \left( {\bf r}_g - {\bf r}_o \right) \times \ddot{\bf r}_o = R^2 \left( \ddot{\theta}\,\cos^2 \theta -\dot{\theta}^2 \sin\theta\,\cos\theta \right) . \tag{2.2} \]
The z-component of the derivative of the angular momentum is
\[ \dot{\bf A}_{oz} = -\frac{4}{3}\,m r^2 \ddot{\theta} . \tag{2.3} \]
The torque with respect to the point of contact of the pencil with the ground is
\[ {\bf tau}_o = - mgR\,\sin\theta \,{\bf k} . \tag{2.4} \]
We substitute Eqs.(2.2), (2.3), and (2.4) into the equation (2.1) and obtain
\[ \left( \frac{R}{g} \right) \ddot{\theta} = \frac{\sin\theta - \left( R/g \right) \dot{\theta}^2 \cos\theta\,\sin\theta}{1/3 + \sin^2 \theta} . \tag{2.5} \]
      We plot the direction field for Eq.(5):
Table[VectorPlot[{y, (0.5*Sin[x] - 1/0.5*y^2*Cos[x]*Sin[x])/(1/3 + Sin[x]^2)}, {x, -2, 2}, {y, -2, 2}, StreamPoints -> Coarse, StreamScale -> {Full, All, 0.05}, StreamStyle -> "Arrow"]]
       Fig. 2.2:     Direction field of the falling pencil on frickless table.            Mathematica code

Now we rederive Eq.(2.5) using another point, C, as an instantaneous center of rotation. The corresponding moment of inertia is

\[ J_c = \frac{m}{3}\,R^2 + m\,R^2 \sin^2 \theta \qquad \Longrightarrow \qquad \dot{J}_c = 2m r^2 \dot{\theta}\,\cos\theta\,\sin\theta . \]
We rewrite the equation of motion (1) for the case of rotational rigid body in plane
\[ \frac{\text d}{{\text d}t} \left( \frac{\omega^2}{2}\, J_c \right) = \frac{\omega^2}{2}\,\dot{ J}_c + \omega\,\dot{\omega} \,J_c = \tau_c \omega , \qquad \omega = \dot{\theta} , \]
where C denotes the instantaneous center, and τc is the torque with respect to this point. This equation can be simplified as follows
\[ \tau_c = \frac{\omega}{2}\,\dot{ J}_c + \ddot{\theta} \,J_c , \qquad \omega = \dot{\theta} . \tag{2.6} \]
Since the torque of the pencil with respect to point C is
\[ \tau_c = mg\,R\, \sin\theta , \]
we obtain from Eq.(2.6) that
\[ mg\,R\, \sin\theta = m R^2 \left( \frac{1}{3} + \sin^2 \theta \right) \ddot{\theta} + mR^2 \dot{\theta}^2 \sin\theta \,\cos\theta , \]
which is equivalent to Eq.(2.5). The first integral (which is the sum of the kinetic and potential energies) of this equation becomes
\[ mR^2 \left( \frac{1}{3} + \sin^2 \theta \right) \frac{1}{2}\,\dot{\theta}^2 + mgR\,\cos \theta = E (= \mbox{K} + \Pi ). \tag{2.7} \]
      We plot the phase trajectories using Mathematica:
ContourPlot[(1 + 3*(Sin[x])^2)*y^2 /2 + 5*Cos[x], {x, 0, 15}, {y, 0, 3}, ContourStyle -> Black]
       Fig. 2.3:     Phase portrait of the falling pencil on ficktionless table.            Mathematica code

   ■

Example 3:    ■

 

  1. Abalmassov, V.A., Maljutin, D.A., Where will a pen fall to?, School Teacher in Physics
  2. Bundy, F.P., Stresses in freely falling chimneys columns, Journal of Applied Physics, 1940, 11, No. 2, pp. 112--123.
  3. Crawford, F.S.,Moments to remember, American Journal of Physics, 1989, 57, Issue 2, p. 177. https://doi.org/10.1119/1.16117
  4. Cross, R., The fall and bounce of pencil and other elongated objects, American Journal of Physics, 2006, 74, Issue 1, pp. 26--30. https://doi.org/10.1119/1.2121752
  5. Cross, R., A falling rod on wheels, Physics Education, 2021, vol. 56, September, 5pp.
  6. Härtel, H., The falling stick with 𝑎 > g, The Physics Teacher, 2000, 38, January, pp. 54--55. https://doi.org/10.1119/1.880430
  7. Madson, E.L., Theory of the chimney breaking while falling, American Journal of Physics, 1977, 45, No. 2, pp. 182--184. https://doi.org/10.1119/1.10651
  8. Oliveira, V., Experiments with a falling rod, American Journal of Physics, 2016, 84, Issue 2, pp. 113--117. https://doi.org/10.1119/1.4934950
  9. Theron, W., The "faster than gravity" demonstration revisited, American Journal of Physics, 1988, 56, Issue 8, pp. 736--739. https://doi.org/10.1119/1.15513
  10. Tiersten, M.S., Moments not to forget---The conditions for equating torque and rate of change of angular momentum around the instantaneous center, American Journal of Physics, 1991, 59, pp. 733--738. https://doi.org/10.1119/1.16752
  11. Tiersten, M.S., Erratum: Moments not to forget---The conditions for equating torque and rate of change of angular momentum around the instantaneous center, American Journal of Physics, 1992, 60, No. 2, p. 187.
  12. Turner, L., Pratt, J.L., Does a falling pencil levitate?, Quantum, 1998, 8, No. 4, pp. 22-25. Online at National Science Teachers Association: https://www.nsta.org/quantum-magazine-math-and-science
  13. Varieschi, G, and Kamiya, K., Toy models for the falling chimney, American Journal of Physics, 2003, 71, No, 9, pp. 1025--1031. https://doi.org/10.1119/1.1576403
  14. Zypman, F.R., Moments to remember---The conditions for equating torgue and rate of change of angular momentum, American Journal of Physics, 1990, 58, No. 1, pp.41--43.

 

  1. Bukov, Martin, D’Alessio, Luca, and Polkovnikov, Anatoli, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to floquet engineering, Advances in Physics, 64(2):139–226, 2015.
  2. Kapitza, P. L., “Dynamic stability of the pendulum with vibrating suspension point, Sov. Phys. – JETP 21, 588–597 (1951) (in Russian). See also Collected papers of P. L. Kapitza, edited by D. Ter Haar ( Pergamon, London, 1965), Vol. 2, pp. 714–726.
  3. F. M. Phelps III and J. H. Hunter, Jr., An analytical solution of the inverted pendulum, American Journal of Physycs, 1965, 33, 285–295 (1965); https://doi.org/10.1119/1.1971474
  4. F. M. Phelps, III and J. H. Hunter, Jr., Am. J. Phys. 34, 533–535 (1966). https://doi.org/10.1119/1.1973087,

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions