The Mathematica Journal: Volume 10, Issue 1: Approximating Solutions of Linear Ordinary Differential Equations with Periodic Coefficients by Exact Picard Iterates

Approximating Solutions of Linear Ordinary Differential Equations with Periodic Coefficients by Exact Picard Iterates

5. Application to Linearized Systems of Driven Pendula

Consider a system of

upside-down physical pendula (homogeneous rods, in fact) such as in the following graphic, where we depict the case

. It is assumed that each rod may rotate without friction around its bearing and that gravity points down.

It is a fact of everyday experience that the upside-down equilibrium position in which all rods are fixed at the vertical position

is highly unstable.

Curiously, however, if the support point of the first rod is driven by an external force such that it executes a vertical sinusoidal motion of amplitude

and angular frequency

, and if

and

are chosen from a suitable range, then the upside-down equilibrium becomes stable. This fact seems to have been first predicted in 1908 by Stephenson [15] for a single rod and subsequently extended by him to the case of two and three rods in 1909 [16]. Interest in such systems was aroused more recently by Acheson, who proved in [7] that if the nonlinear equations of motion for the rods are linearized in the neighborhood of the equilibrium solution, then, for any value of

there exists a stability region for this system of linearized upside-down rods. In other words, there exists a range in

and

where the upside-down equilibrium for the linearized equations is stable. If the values of the angles

and their derivatives are small, one may hope that the solution of the linearized equations is close to the solution of the nonlinear equations, but a mathematical proof of stability for nonlinear rods is much more difficult [9]. In fact, even an experimental proof of stability for systems of one, two, and three rods has merited a good deal of attention recently [8]. Since then, the number of papers in this field has been very large and even outside of the physics and mathematics communities. Stability of upside-down rods (or pendula) have attracted interest in areas ranging from control engineering to biomechanics. A good reference from the mathematical physics point of view with many further references is [9].

The Lagrangian

In this subsection we use the Lagrangian formalism of classical mechanics [17] to deduce the equations of motion for undamped systems of periodically driven physical pendula and then linearize them around the upside-down equilibrium position. The Euler-Lagrange equations of motion are

where the Lagrangian function

is the difference between kinetic and potential energies.

Consider a system such as the one in the previous graphic, where the

th pendulum has length

and mass

. Suppose that the mass of the

th pendulum is distributed such that its center of mass is located at a distance

from its bottom extremity and its inertia moment through an axis passing through the bottom extremity and orthogonal to the plane of motion is

, that is, the gyration radius is

. Suppose also that the support point of the lower pendulum is externally driven so that its position vector is

. We will denote with

the velocity vector of the bottom extremity of the

th pendulum with respect to the inertial lab frame. The velocity vector of the top of pendulum

with respect to the noninertial frame of its bottom will be denoted with

. Using well-known expressions for

and

, the kinetic and gravitational potential energies of a rigid body, and

for the acceleration of gravity, we obtain as follows the Lagrangian

for the system.

A Single Rod

For a rod with uniform mass distribution

and sinusoidal (in fact, cosinusoidal) vertical driving, the equation of motion is

It is possible to put the equation of motion in a simpler form, depending on a smaller number of parameters, if we rescale the time variable defining

So, the question of stability of an upside-down vertically driven single rod amounts to stability of the solutions of the Mathieu equation. This has been studied for a long time [18].

The Floquet theory [18] for linear differential equations with periodic coefficients examines, among other things, the stability of equilibrium solutions. It is proved that the equilibrium solution

for a linear system with periodic coefficients is stable if and only if all eigenvalues

of the Floquet transition matrix (FTM) are such that

In the particular case of Mathieu equations, one proves that the transition from stability to instability occurs when both eigenvalues are equal to 1 or both equal to

. In other words, the curves in the

plane where solutions to the Mathieu equation change from stable to unstable or vice versa are given by the condition that

Before we use this condition to produce a plot of the stability regions for a single inverted rod, let us make a numerical comparison of our results with exact ones. We will use the approximation ftmmathieu to the FTM calculated previously with 20 iterates and choose two values of

and

, in order to find the values of

corresponding to the stability boundaries at these values of

. We do that as follows.

For the sake of comparison, the exact values corresponding to the approximate ones, are related to the so-called Mathieu characteristic numbers implemented in Mathematica.

Instead of plotting the curves where equation (15) holds in variables

and

, it is interesting to use physical variables related to the pendula. So we define

(dimensionless amplitude) and

(dimensionless frequency) as

The relationship between these variables and

and

. In order to get a good graphic, we need to use a large value for the option PlotPoints of ContourPlot. It takes a reasonable time to run, but the resulting graph is as follows.

We claim that the inverted rod is stable if the frequency and amplitude of the driving are such that the corresponding point lies in the region between the two curves. To see that, we choose an arbitrary point in this region and check stability by computing the corresponding eigenvalues of the FTM. For example, the point

coordinates lies in the mentioned region. The absolute values of the eigenvalues of the FTM are thus.

This confirms our claim! We see that for reasonably small amplitudes, such as

, the inverted rod is stable (in the linear approximation) if the frequency is high enough. For larger amplitudes, the inverted rod is stable only in a precise range of frequencies.

Two Rods

The left-hand sides of the Euler-Lagrange equations of motion for a system of two pendula, restricted to the simpler case of equal masses, equal lengths, and uniform mass distributions follow.

Notice the presence of several nonlinear terms. Again, the first step in linearizing the equations around the equilibrium

is to replace sines and cosines by their Taylor expansions, retaining terms of at most order 1. But now, this does not eliminate all nonlinearities. In order to finish the linearization procedure, we first create a function to detect the degree of a monomial in a given variable and a Boolean function that tests for monomials of degree not larger than 1.

Then we use it to select only the terms of degree not larger than 1 in all variables

We would like to write the previous linear expression in matrix form. In order to obtain the matrices, we use the command LinearExpressionToMatrix in the Developer` context (or, equivalently, transform the expression into a system of equations and then use the LinearEquationstoMatrices command in the standard package LinearAlgebra`MatrixManipulation`).

Notice that the remainder vanishes in the case of vertical driving, that is,

. Using the cosine vertical driving

, the term proportional to

becomes

we can rewrite it into the form of equation (4) by making

, with the matrix

in equation (4) being

and the dimensionless variables

and

are the same as defined in equations (16) and (17).

In order to produce a picture of the stability boundaries for the double rod in the dimensionless variables, in analogy with what was done in the case of one rod, we must locate the points where (in absolute value) the maximum eigenvalue of the Floquet matrix changes from greater than one to equal to one. In the case of the single rod, this was conveniently accomplished by the trace condition of equation (15), but that does not hold anymore.

As we do not know any condition similar to equation (15) to use instead, the best we can do is locate some points in the stability boundaries. We may do that by using an idea similar to the bisection method for solving equations. Whenever the largest eigenvalue (in absolute value) of the FTM is greater than one at some point

and equal to one at another point

, then there must exist a point between

and

through which a stability boundary passes. We may approximate this point by taking the point

halfway between

and

(hence the name bisection) and calculating at

the maximum eigenvalue of the FTM in absolute value. If this number is equal to one, then the stability boundary must pass between

and

, otherwise it must pass between

and

. By repeating the bisection procedure a sufficient number of times we may approximate a point in the stability boundary between

and

with an error smaller than any prescribed tolerance.

By experience we know that the stability boundary lies in the region

. We then produce a grid of points in that region and calculate the largest eigenvalue of the FTM in absolute value at each point in the grid.

In order to select those pairs of points in the grid between which we know that the stability boundary passes, we create the following Boolean function and use it.

Notice that because of the particular way pointgrid was produced, we have selected pairs of points in the grid with the same horizontal coordinate, that is, points on vertical lines. Later we will repeat the procedure for selecting points lying on horizontal lines.

We can now create the function bisect that takes each pair of points in the output of the previous line and calculates the midpoint of the pair, the largest eigenvalue in absolute value of the FTM at the midpoint, and then selects from among the three points the two between which the stability boundary must lie.

If we want to locate the point in the stability boundary between each pair of selected points, it is necessary to iterate bisect the right number of times. If the tolerance tol is specified, that can be done by the following function.

So, this finds with tolerance 0.1 (in the vertical direction) some points in the stability boundaries.

We may now repeat the procedure using Transpose[pointgrid] instead of pointgrid. That will produce pairs of points lying on horizontal lines between which some point in the stability boundaries must exist. Of course, as the scale of the horizontal axes in the picture is smaller, we must accordingly reduce the tolerance tol.

An alternative exists for generating this graph: we may use our Mathieu equation results together with the introduction of normal coordinates that decouple equation (18) into two separate Mathieu equations. These normal coordinates, which we now briefly explain, were used by Acheson in the proof of his pendulum theorem [7].

Notice that as the matrix

has two distinct real eigenvalues, then it must be diagonalizable over the real numbers. Let

be the orthogonal matrix that diagonalizes

. We define the normal coordinates

. Substituting

in equation (18) and then left multiplying it by

, we obtain

where

and

are the eigenvalues of

. These are two Mathieu equations (7) with

and

As the original variables

are linear combinations of the normal variables

, then the equilibrium solution of equation (18) will be stable if and only if the equilibrium solution for both equations (21),

is stable. We may then determine the stability region for the linearized double rod in equation (18) by intersecting the two stability regions for Mathieu equations. This can be done by using the already calculated approximate FTM for the Mathieu equation. In the following result, the red curve is the stability picture for the mode

, the cyan curve corresponds to the mode

, and the points over the curves are the ones calculated previously. Notice that both ways of producing the picture lead to indistinguishable results!

In order to check the numerical accuracy of one of the points in the picture, we may use the built-in Mathieu characteristic numbers together with the normal coordinates method.

Given that the first coordinate is 0.06, here is the exact value for the second coordinate of the chosen point.

This is in perfect agreement with the value obtained by our method. Remember that the second coordinate of that point was calculated by the bisection method with a tolerance of 0.1.

It can be shown that the method of normal coordinates works for any number of rods, not only for two. The procedure for producing the previous graph for three or more rods is the same and we do not reproduce it here, leaving it as an exercise for the reader.