Exploring the stability of inverted multi-linked pendulums

Worldwide Mathematica Conference 1998

Peder Thusgaard Ruhoff

The Maersk Mc-Kinney Moller Institute for Production Technology

Odense University, Denmark

Initialization

PendulumSystem = 

      Show[Graphics[{

            Text[StyleForm["q[1][t]", FontFamily -> "Courier",  

                FontWeight -> "Bold"], {(4/3) Sin[(π/12)], 

                ((-4)/3) Cos[(π/12)]}, {0, 0}], 

            Text[StyleForm["l[1]", FontFamily -> "Courier",  

                FontWeight -> "Bold"], {Sin[(π/4)], 

                -Cos[(π/4)]}, {0, 0}], 

            Text[StyleForm["m[1]", FontFamily -> "Courier",  

                FontWeight -> "Bold"], {(7/3) Sin[(π/5)], 

                ((-7)/3) Cos[(π/5)]}, {0, 0}], {Thickness[0.015], 

              Line[{{0, 0}, {2  Sin[(π/6)], -2 Cos[(π/6)]}}]}, 

            Disk[{2  Sin[(π/6)], -2 Cos[(π/6)]}, 0.1], {

              Dashing[{0.03, 0.03}], 

              Circle[{0, 0}, 

                2, {((3  π)/2) - (π/6), 

                  ((3  π)/2) + ((2  π)/6)}]}, 

            Circle[{0, 0}, 

              1, {((3  π)/2), ((3  π)/2) + (π/6)}]}, 

          AspectRatio -> Automatic,  Axes -> True, 

          Ticks -> None]]

\(MultiPendulumSystem = 

      Show[Graphics[{

            Text[StyleForm["q[1][t]", FontFamily -> "Courier", \ 

                FontWeight -> "Bold"], {\(4\/3\) Sin[π\/12], 

                \(\(-4\)\/3\) Cos[π\/12]}, {0, 0}], 

            Text[StyleForm["l[1]", FontFamily -> "Courier", \ 

                FontWeight -> "Bold"], {Sin[π\/4], 

                \(-Cos[π\/4]\)}, {0, 0}], 

            Text[StyleForm["m[1]", FontFamily -> "Courier", \ 

                FontWeight -> "Bold"], {\(7\/3\) Sin[π\/5], 

                \(\(-6\)\/3\) Cos[π\/5]}, {0, 0}], {Thickness[0.015], 

              Line[{{0, 0}, {2  Sin[π\/6], \(-2\) Cos[π\/6]}}]}, 

            Disk[{2  Sin[π\/6], \(-2\) Cos[π\/6]}, 0.1], 

            Circle[{0, 0}, 1, {\(3  π\)\/2, \(3  π\)\/2 + π\/6}], 

            
			Disk[{2  Sin[π\/6], \(-2\) Cos[π\/6]} + {

                  2  Sin[π\/3], \(-2\) Cos[π\/3]}, 0.1], 
				

            Line[{{2  Sin[π\/6], \(-2\) Cos[π\/6]}, {

                  2  Sin[π\/6], \(-2\) Cos[π\/6] - 1.2}}], 
				

            Circle[{2  Sin[π\/6], \(-2\) Cos[π\/6]}, 

              1, {\(3  π\)\/2, \(3  π\)\/2 + π\/3}], 
				

            Text[StyleForm["q[2][t]", FontFamily -> "Courier", \ 

                FontWeight -> "Bold"], {1.6, \(-2.8\)}, {0, 0}], 


            Text[StyleForm["l[2]", FontFamily -> "Courier", \ 

                FontWeight -> "Bold"], {2.0, \(-2.0\)}, {0, 0}], 

            Text[StyleForm["m[2]", FontFamily -> "Courier", \ 

                FontWeight -> "Bold"], {2.8, \(-2.3\)}, {0, 0}], 


            	{Thickness[0.015], 

              Line[{{2  Sin[π\/6], 

                    \(-2\) Cos[π\/6]}, {2  Sin[π\/6], 

                      \(-2\) Cos[π\/6]} + {2  Sin[π\/3], 

                      \(-2\) Cos[π\/3]}}]}		}, 

          AspectRatio -> Automatic, \ Axes -> True, 

          Ticks -> None]]\)

Introduction

The purpose of this talk is not so much to demonstrate new scientific 
discovery. Rather, I would like to demonstrate how I prefer to work with Mathematica.

First, however, I will take a look back at "my life with Mathematica"

1989: First encounter with Mathematica. My physics professor asked me to take a look at this new program for 
doing mathematics by computer. I was hooked right away. Problem: Only user at the university - no one to discuss with.

1991: Completed my master thesis in Physics and Computer Science  (Mathematica version 1.2). The subtitle of the thesis was applications of Mathematica in molecular physics. Used Mathematica for Van Vleck perturbation theory, operator algebra, numerical solution 
of the Schroedinger equation, etc.

1995: I was responsible for the introduction of Mathematica in the curriculum for the education in applied mathematics and computer 
systems enginnering at our university.

A simple pendulum

Anyone taking a course in classical mechanics have been introduced to the 
simple pendulum. In figure 1 we show such a system.

Note that the pendulum can only swing in the vertical plane; the system is 
two-dimensional.

The Lagrangian for the system can be written

The equation of motion for the pendulum is derived in the usual manner

eqm[1]  = D[D[L[1], q[1]'[t]], t] - D[L[1], q[1][t]] == 0

We can go ahead and solve this second order differential equation 
numerically using the built-in function NDSolve. In order to do so we must  insert a numerical values for the constants 
and specify a set of initial conditions.

sol[1]  =  

      NDSolve[{eqm[1]  /. {m[1] -> 1, l[1] -> 1, g -> 9.82}, 

          q[1][0] == 0, q[1]'[0] == 3}, q[1][t], {t, 0, 20}]

We can plot the solution in various ways.

Plot[Evaluate[q[1][t] /. First[sol[1]]], {t, 0, 20}, 

      Frame -> True]

ParametricPlot[

      Evaluate[{q[1][t] /. sol[1][[1]], 

          D[q[1][t] /. sol[1][[1]], t]}], {t, 0, 20}, 

      PlotPoints -> 200, AspectRatio -> Automatic]

ParametricPlot3D[

      Evaluate[{t/10, Sin[q[1][t]], -Cos[q[1][t]]} /. 

          sol[1][[1]]], {t, 0, 20}, PlotPoints -> 200, 

      ViewPoint -> {3.032,  -1.481,  0.245}]

A multi-linked pendulum

We can easily generalize the simple pendulum. In the general case of an n-linked pendulum we can write the Lagrangian function as

where

Again, the motion can only take place in a vertical plane. The system is 
shown below in the case of the double pendulum.

The equations of motion can be derived like we did for the simple pendulum

Let us consider a particular example; namely, a 4-linked pendulum. For 
simplicity, we will assume that all masses are equal to 1 Kilogram and all 
lengths are equal to 1Meter.

sol[2] = 

      NDSolve[Flatten[{

            eqm[4] /. {m[_] -> 1. , l[_] -> 1. , g -> 9.82},  

            Table[{q[i][0] == π - i 0.01, q[i]'[0] == 0}, {i, 1, 

                4}]}], Table[q[i][t], {i, 1, 4}], {t, 0, 5}, 

        MaxSteps -> 10000]

We can plot the result

Plot[

      Evaluate[{q[1][t], q[2][t], q[3][t], q[4][t]} /.  

          First[sol[2]]], {t, 0, 5}, PlotStyle -> ColorList[4], 

      Frame -> True]

Alternatively, one can generate an animation of the system.

Table[

      Show[Pendulum[{1, 1, 1, 1}, 

          Evaluate[{q[1][t], q[2][t], q[3][t], q[4][t]} /.  

              First[sol[2]]], {0, 0}], AspectRatio -> Automatic, 

        PlotRange -> {{-4, 4}, {-4, 4}}, Frame -> True], {

        t, 0, 5, 0.10}]

A driven multi-linked pendulum

The systems considered sofar are charaterized by having a fixed suspension 
point. However, it is straight forward to include the situation where the 
suspension is allowed to move in some specific manner. Assuming the 
suspension point follows the curve (x(t), y(t)) = (ξ(t), ψ(t)), the Lagrangian function can be written

where T is as before and

The equation of motion is readily found from the Lagrangian function, we 
get

In order to numerically integrate the equation of motion, we must insert 
numerical values for the constants, give explicit expressions for ξ[t] and ψ[t], and specify a set of initial conditions.

sol[3] = 

      NDSolve[Flatten[{

            deqm[4] /. {m[_] -> 1. , l[_] -> 1. , g -> 9.82}, 

             Table[{q[i][0] == 0, q[i]'[0] == 0}, {i, 1, 4}]}], 

        Table[q[i][t], {i, 1, 4}], {t, 0, 2  π}, 

        MaxSteps -> 10000]

Plot[

      Evaluate[{q[1][t], q[2][t], q[3][t], q[4][t]} /.  

          First[sol[3]]], {t, 0, 2  π}, PlotStyle -> ColorList[4], 

      Frame -> True]

Table[

      Show[{Pendulum[{1, 1, 1, 1}, 

            Evaluate[{q[1][t], q[2][t], q[3][t], q[4][t]} /. 

                 First[sol[3]]], {ξ[t], ψ[t]}], 

          Graphics[

            Line[Table[{ξ[t], ψ[t]}, {t, 0, π, π/50}]]]}, 

        AspectRatio -> Automatic, 

        PlotRange -> {{-4, 4}, {-5, 3}}, Frame -> True], {

        t, 0, 2  π, 0.10}]

Acheson´s general pendulum theorem

In 1993 Acheson proved the following theorem on the stability of inverted 
pendulum systems.

Let N pendulums hang down, one from another, under gravity g, each having one 
degree of freedom, the uppermost being suspended from a pivot point O. Let (Subscript[ω, max]) og (Subscript[ω, min]) denote the largest and the smallest of the natural frequncies of small 
oscillation about this equillibrium state.

Now turn the whole system upside-down. The resulting configuration of the 
pendulums can be stabilized (according to linear theory, at least) if we 
subject the pivot point O to vertical oscillations of suitable amplitude ϵ and frequency (Subscript[ω, 0]). When \(TraditionalForm

      \`ω\_0\%2 ≫ \ ω\_max\%2\)the stability criterion is

\(TraditionalForm

    \`\(\(\@2\) g\)\/\(\(ω\_0\) ω\_min\) < ϵ < 

      \(0.450  g\)\/ω\_max\%2\)

The proof of the theorem is simple. It is based on the fact that the 
equations of motion for the system (in a linear approximation) can be written 
as system of N uncoupled Mathieu equations.

Natural frequencies

It is straight forward to determine the natural frequencies of the system 
in its downward equillibrium configuration (see, e.g Landau and Lifshitz: Mechanics). For the present work we only need to know the smallest and largest of 
the frequencies.

First we introduce a function for computing the charateristic polynomial

Then the natural frequencies are computed using the following

An example

NaturalFrequencies[{1, 1, 2}, {1, 1, 2}]

We now make a function that returns the stability criterion for a given 
multi-linked pendulum

We can now make a small table

TableForm[

      Table[{i, StabilityCriterion[Table[1, {i}], Table[1, {i}]]}, {i, 1, 4}]]

Testing the stability criterion

Finally, let us try out a few cases.

First we note that the suspension point must now follow a curve of the 
form \(TraditionalForm

      \`\((x(t), y(t))\) = \((0, \ ϵ\ \(cos(ω\_0)\))\)\), that is

For the double pendulum with masses m[1]=m[2]=1 and lengths l[1]=l[2]=1, the stability criterion is

StabilityCriterion[{1, 1}, {1, 1}]

Thus

should satisfy the criterion. We can, of course, test this easily

StabilityCriterion[{1, 1}, {1, 1}] /.  val[1]

We now solve the equations of motion

sol[5] = 

      NDSolve[Flatten[{

            deqm[2] /. {m[_] -> 1. , l[_] -> 1. , g -> 9.82}

                   /. val[1],  

            Table[{q[i][0] == π - 0.2, q[i]'[0] == 0}, {i, 1, 

                2}]}], Table[q[i], {i, 1, 2}], {t, 0, 10}, 

        MaxSteps -> 10000]

Plot[

      Evaluate[{q[1][t], q[2][t]} /.  First[sol[5]]], {t, 0, 10}, 

      PlotStyle -> ColorList[2], Frame -> True, 

      GridLines -> {None, {π}}]

Now, let us move outside the region of stability

Indeed the values does not satisfy the criterion

StabilityCriterion[{1, 1}, {1, 1}] /.  val[2]

We now solve the equations of motion

sol[6] = 

      NDSolve[Flatten[{

            deqm[2] /. {m[_] -> 1. , l[_] -> 1. , g -> 9.82}

                   /. val[2],  

            Table[{q[i][0] == π - 0.2, q[i]'[0] == 0}, {i, 1, 

                2}]}], Table[q[i], {i, 1, 2}], {t, 0, 10}, 

        MaxSteps -> 10000]

Plot[

      Evaluate[{q[1][t], q[2][t]} /.  First[sol[6]]], {t, 0, 10}, 

      PlotStyle -> ColorList[2], Frame -> True, 

      GridLines -> {None, {π}}]

Let us try and animate the two results

Table[

      Show[{Graphics[{Line[{{- .1, 0}, { .1, 0}}], 

              Line[{{0, -0.1}, {0, 0.1}}]}], 

          Pendulum[{1, 1}, 

            Evaluate[{q[1][t], q[2][t]} /.  First[sol[5]]], {

                ξ[t], ψ[t]} /. val[1]], 

          Pendulum[{1, 1}, 

            Evaluate[{q[1][t], q[2][t]} /.  First[sol[6]]], {

                ξ[t], ψ[t]} /. val[2]]}, 

        AspectRatio -> Automatic, 

        PlotRange -> {{-3, 3}, {-3, 3}}, Frame -> True], {

        t, 0, 10, 0.10}]

Summary

We have demonstrated how Mathematica can be used to investigate the stability of inverted multi-linked 
pendulums in an easy manner. We have shown how one working with Mathematica can refine or generalize a model but adding new features.