Volume 10, Issue 1 Articles Trott's Corner New Products New Publications Calendar News Bulletins New Resources Classifieds Download This Issue Editorial Policy Staff and Contributors Submissions Subscriptions Advertising Back Issues Contact Information 
Approximating Solutions of Linear Ordinary Differential Equations with Periodic Coefficients by Exact Picard Iterates
5. Application to Linearized Systems of Driven PendulaConsider a system of upsidedown 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 upsidedown 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 upsidedown 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 upsidedown rods. In other words, there exists a range in and where the upsidedown 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 upsidedown 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 LagrangianIn 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 upsidedown equilibrium position. The EulerLagrange 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 wellknown 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.
L[n] will produce the Lagrangian for a system of pendula. A Single RodHere is the Lagrangian for a single upsidedown driven rod, for example.
From it we can evaluate the lefthand side of equation (10).
For a rod with uniform mass distribution and sinusoidal (in fact, cosinusoidal) vertical driving, the equation of motion is
Here linearization is substituting by .
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
the new form assumed by the equation is then
which is the Mathieu equation (7) with
and
So, the question of stability of an upsidedown 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 socalled 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
and
The relationship between these variables and and is , . 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 in 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 RodsThe lefthand sides of the EulerLagrange 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`).
The matrix that serves as the coefficient for the vector is obtained as
The coefficient of the vector is
and the remaining terms are
Notice that the remainder vanishes in the case of vertical driving, that is, . Using the cosine vertical driving , the term proportional to becomes
and again rescaling the time as in equation (11), we arrive finally at
where was previously obtained as
and is the diagonal matrix
Multiplying both sides of equation (18) on the left by and realizing that
we can rewrite it into the form of equation (4) by making , with the matrix in equation (4) being
and obeying equation (5) with . At this point, we can calculate the twentieth approximation to the FTM as
where we should have in mind that
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 as . Substituting by in equation (18) and then left multiplying it by , we obtain
which, because is diagonal, decouples into
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 builtin 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.


About Mathematica  Download Mathematica Player © Wolfram Media, Inc. All rights reserved. 