13: Pendulum Dynamics
( \newcommand{\kernel}{\mathrm{null}\,}\)
The undriven pendulum has only a two dimensional space, forming a phase plane where it is easy to visualize the phase portraits. The phase portrait of the small amplitude simple pendulum is particularly easy to draw. The dimensionless form of the simple pendulum equation is
¨θ+θ=0,
and the solutions for θ=θ(t) and u=˙θ(t) are given by
θ(t)=θmcos(t+φ),u(t)=−θmsin(t+φ)
The phase portrait in the θ−u phase plane consists of concentric circles, with clockwise motion along these circles, as seen in the small diameter circles of Fig. 13.1. As the amplitude increases, the approximation sinθ≈θ loses validity, and the relevant dimensionless equation becomes
¨θ+sinθ=0.
Beyond the small diameter circles of Fig. 13.1, one observes that the connected lines become elongated in u as the amplitude increases, implying the pendulum is slowed down when the amplitude becomes large (i.e., the period of the pendulum is lengthened). Eventually, a final closed curve is drawn, called the separatrix, that separates the pendulum’s motion from periodic to rotary, the latter motion corresponding to a pendulum that swings over the top in a circular motion. Exactly on the separatrix trajectory, the pendulum comes to rest at the angle π, or 180∘, which is an unstable fixed point of the pendulum.
The simple pendulum is a conservative system, exhibiting a conservation law for energy, and this implies a conservation of phase space area (or volume). If one evolves over time a given initial area of the phase space of Fig. 13.1, the area of this phase space will be conserved. When the pendulum is damped, however, the area of this phase space will shrink to zero.
Basin of attraction of the undriven pendulum
We consider the damped, undriven pendulum with governing equation
¨θ+1q˙θ+sinθ=0
and the stable fixed point given by (θ,˙θ)=(0,0). How can we determine the basin of attraction of this fixed point? Of course, with ˙θ=0, only the values −π<θ<π will lie in the basin of attraction. But we also need to compute the basin of attraction for nonzero initial values of ˙θ.
As is often the case, to devise a numerical algorithm it is best to appeal directly to the physics. We want to find the borderline of initial conditions between either the attractive fixed points (0,0) and (2π,0), or the attractive fixed points (0,0) and (−2π,0). The initial conditions just on the border result in the pendulum reaching either the unstable fixed points (π,0) or (−π,0). A small

perturbation from (π,0) will result in one of the attractive points (0,0) or (2π,0); from (−π,0), one of the attractive points (0,0) or (−2π,0).
The algorithm then initializes the calculation at either the unstable fixed point (π,0) or (−π,0), with a small perturbation in either the position or velocity that results in the final attracting point being (0,0). We can call these values final conditions because the differential equations are then integrated backward in time to determine all possible initial conditions that can result in these final conditions. You can convince yourself that the four final conditions that should be used are the two pairs (π,−ϵ),(π−ϵ,0), and (−π,ϵ),(−π+ϵ,0), with ϵ very small. The first of each pair corresponds to the immediately previous motion being rotary, and the second of each pair corresponds to the immediately previous motion being oscillatory.
A graph of the basins of attraction of (0,0) for q=4, corresponding to an underdamped pendulum (critical damping is q=1/2 ) is shown in Fig. 13.2. The attracting fixed point is marked by an ’ x ’, and the region inside the two curved lines is the basin of attraction.
Spontaneous symmetry-breaking bifurcation
An interesting supercritical pitchfork bifurcation occurs in the pendulum equations, where at the bifurcation point one stable limit cycle splits into two. The single limit cycle displays a symmetry that is no longer respected individually by the two new limit cycles. This type of pitchfork bifurcation is a manifestation of what is called spontaneous symmetry breaking.

The pendulum equation (11.14), written again here, is given by
¨θ+1q˙θ+sinθ=fcosωt.
Using sin(−θ)=−sinθ and cos(ωt−π)=−cosωt, the pendulum equation can be seen to be invariant under the transformation
θ→−θ,t→t−π/ω
with the physical interpretation that the equations of motion make no distinction between the right and left sides of the vertical, a consequence of the symmetry of both the pendulum and the external force.
Consider again the asymptotic small amplitude solution given by (11.10), which in dimensionless variables is
θ(t)=f√(1−ω2)2+(ω/q)2cos(ωt+ϕ),
with
tanϕ=ω/qω2−1
We can observe that this solution is also invariant under (13.6) : the small amplitude solution obeys θ(t)=−θ(t−π/ω). We say that this solution is symmetric, meaning it obeys the same symmetry as the governing equations; that is, the motion of the small amplitude pendulum is symmetrical about the vertical.

also satisfies (13.5) by the following calculation:
¨θ2(t)+1q˙θ2(t)+sin(θ2(t))=−¨θ1(t−π/ω)−1q˙θ1(t−π/ω)−sin(θ1(t−π/ω))=−fcos(ωt−π)=fcosωt
If the two solutions θ1(t) and θ2(t) are equal, then we say that this solution is symmetric. Here, we are considering asymptotic solutions that are independent of the initial conditions, since the initial conditions themselves can also break the symmetry. However, if θ1(t) and θ2(t) are not equal, we say that these solutions are asymmetric, and that spontaneous symmetry breaking has occurred. Evidently, spontaneous symmetry breaking is a decidedly nonlinear phenomena. After symmetry breaking occurs, the asymmetric solutions must occur in pairs, and the bifurcation point looks like a pitchfork bifurcation, and can be super- or sub-critical. Any subsequent bifurcation that occurs to one of the asymmetric solutions must also be mirrored by the other.
Spontaneous symmetry breaking occurs in the pendulum dynamics at the dimensionless parameter values f=1.5,ω=2/3, and approximately q=1.246. For q just less than 1.246, the single stable asymptotic solution for θ=θ(t) is symmetric, and for q just greater than 1.246, there exists a pair of stable asymptotic solutions that are asymmetric.
By projecting the phase-space trajectory onto the θ−˙θ plane, in Fig. 13.3 we display both the symmetric solution when q=1.24 and the two asymmetric solutions when q=1.30. To explicitly observe the symmetry of the solutions, we mark the value of (θ,θ) at the time t= n2π/ω, with n a positive integer (equivalent to the time t=0), and the value of (−θ,−˙θ) at the time t=n2π/ω+π/ω (equivalent to the time t=π/ω ). For a symmetric solution, these two points mark the same point on the trajectory, and for asymmetric solutions they mark points on different trajectories. Notice that after the occurrence of symmetry breaking, one of the asymptotic solutions undergoes an oscillation centered to the right of the vertical, and the other, centered to the left.


We can graph a bifurcation diagram associated with the symmetry breaking of the solutions. We fix f=1.5 and ω=2/3 and vary q across the bifurcation point q=1.246. We need to distinguish the symmetric from the asymmetric limit cycles, and one method is to compute the average value of θ(t) over one period of oscillation; that is,
⟨θ⟩=1T∫T0θdt
where T=2π/ω. For the symmetric solution, ⟨θ⟩=0, whereas ⟨θ⟩ takes on both positive and negative values after symmetry breaking occurs. In Fig. 13.4, we plot the value of ⟨θ⟩ versus q. At a value of approximately q=1.246, spontaneous symmetry breaking occurs and the stable symmetric limit cycle splits into two asymmetric limit cycles, in what is evidently a supercritical pitchfork bifurcation.
The symmetric limit cycle still exists after the bifurcation point, and though unstable, can also be computed. To compute the unstable cycle, one could determine the values of θ and ˙θ at t=0 that lie on this cycle, and then integrate over one period (over which the instability doesn’t have sufficient time to develop). The problem of determining the correct initial conditions can be cast as a problem in multidimensional root-finding. The key idea is that a symmetric limit cycle satisfies θ(t)=−θ(t−π/ω). We therefore determine the solution vector of initial conditions (θ(0),θ(0)) that satisfies the two equations
θ(0)+θ(π/ω)=0˙θ(0)+˙θ(π/ω)=0
where θ(π/ω) and ˙θ(π/ω) are determined by integrating from t=0 the differential equations using ode 45.m with the initial conditions (θ(0),˙θ(0)). Root-finding can be done using either a two-dimensional version of the secant rule or, more simply, the built-in MATLAB function fsolve.m. Convergence to the roots is robust, and an initial guess for the roots can be taken, for instance, as (θ(0),˙θ(0))=(0,0). A plot of the two stable asymmetric limit cycles, and the unstable symmetric limit cycle when f=1.5,ω=2/3 and q=1.3 is shown in Fig. 13.5.
Period-doubling bifurcations



The period-doubling bifurcations continue with increasing q. The second doubling from period two to period four occurs at approximately q=1.370 and the third doubling from period four to period eight occurs at approximately q=1.375. The θ−˙θ phase-space projections for period four and period eight are shown in Fig. 13.8.
A bifurcation diagram can be plotted that illustrates these period-doubling bifurcations. We use a Poincaré section to plot the value of θ at the times corresponding to 2πn/ω, with n an integer. The control parameter for the bifurcation is q, and in Fig. 13.9, we plot the bifurcation diagram for 1.34<q<1.38. In general, the angle θ should be mapped into the interval −π< θ<π, but for these parameter values there is no rotary motion. Period-doubling bifurcations are observed, and eventually the pendulum becomes aperiodic. Additional windows of periodicity in the aperiodic regions of q are also apparent.
The aperiodic behavior of the pendulum observed in Fig. 13.9 corresponds to a chaotic pendulum. If only the pendulum exhibited the period-doubling route to chaos, then the results shown in Fig. 13.9, though interesting, would be of lesser importance. But in fact many other nonlinear systems also exhibit this route to chaos, and there are some universal features of Fig. 13.9, first discovered by Feigenbaum in 1975. One of these features is now called the Feigenbaum constant, δ

Note that for the pendulum, we can already compute
δ1=1.370−1.3481.375−1.370=4.400
The meaning of δ can be better elucidated by writing
qn+2−qn+1=qn+1−qnδn
and by continuing to iterate this equation, we obtain
qn+2−qn+1=q2−q1δ1δ2⋯δn.
With all the δn ’s approximately equal to δ, we then have the scaling
qn+2−qn+1∝δ−n
A δ larger than one would then insure that the bifurcations occur increasingly closer together, so that an infinite period (and chaos) is eventually attained. The value of δ can be computed to high accuracy and has been found to be
δ=4.669201609102990…,
and our value of δ1=4.4 is a rough approximation.
Period doubling in the logistic map
Feigenbaum originally discovered the period-doubling route to chaos by studying a simple onedimensional map. A one-dimensional map with a single control parameter μ can be written as
xn+1=fμ(xn)
where fμ(x) is some specified function. A one-dimensional map is iterated, starting with some initial value x0, to obtain the sequence x1,x2,x3,…. If the sequence converges to x∗, then x∗ is a stable fixed point of the map.
The specific one-dimensional map we will study here is the logistic map, with
fμ(x)=μx(1−x)
The logistic map is perhaps the simplest nonlinear equation that exhibits the period-doubling route to chaos. To constrain the values of xn to lie between zero and unity, we assume that 0<μ<4 and that 0<x0<1.
A period-1 cycle for the logistic map corresponds to a stable fixed point. Stable fixed points are solutions of the equation x=fμ(x), or
x=μx(1−x).
The two solutions are given by x∗=0 and x∗=1−1/μ. The first fixed point x∗=0 must be stable for 0<μ<1, being the only fixed point lying between zero and one that exists in this range. To determine the stability of the second fixed point, we make use of the linear stability analysis discussed in §12.1.
For a one-dimensional map, x∗ is a stable fixed point of (13.17) if |f′μ(x∗)|<1. For the logistic map given by (13.18), f′μ(0)=μ so that x∗=0 is stable for 0<μ<1 as we have already surmised, and for the second fixed point
f′μ(1−1/μ)=2−μ.
Therefore, we find that x∗=1−1/μ is stable for 1<μ<3.
What happens when μ becomes larger than three? We will now show that the first perioddoubling bifurcation occurs at μ=3. Because of the simplicity of the logistic map, we can determine explicitly the period-2 cycle. Consider the following composite map:
gμ(x)=fμ(fμ(x))
Fixed points of this map will consist of both period-1 cycles and period-2 cycles of the original map (13.18). If x=x∗ is a fixed point of the composite map (13.21), then x∗ satisfies the equation
x=gμ(x)
A period-2 cycle of the map fμ(x) necessarily corresponds to two distinct fixed points of the composite map (13.21). We denote these two fixed points by x0 and x1, which satisfy
x1=fμ(x0),x0=fμ(x1).
We will call x0,x1 the orbit of the period-2 cycle, with the later generalization of calling x0,x1, …,x2n−1 the orbit of the period −2n cycle.
The period-2 cycle can now be determined analytically by solving
x=fμ(fμ(x))=fμ(μx(1−x))=μ(μx(1−x))(1−μx(1−x))
which is a quartic equation for x. Two solutions corresponding to period-1 cycles are known: x∗=0 and x∗=1−1/μ. Factoring out these two solutions, the second by using long division, results in the quadratic equation given by
μ2x2−μ(μ+1)x+(μ+1)=0.
The period-2 cycle, then, corresponds to the two roots of this quadratic equation; that is,
x0=12μ((μ+1)+√(μ+1)(μ−3)),x1=12μ((μ+1)−√(μ+1)(μ−3)).
These roots are valid solutions for μ≥3. Exactly at μ=3, the period-2 cycle satisfies x0=x1= 2/3, which coincides with the value of the period-1 cycle x∗=1−1/μ=2/3. At μ=3, then, we expect the fixed point of the composite map corresponding to a period-1 cycle to go unstable via a supercritical pitchfork bifurcation to a pair of stable fixed points, corresponding to a period-2 cycle.
When does this period-2 cycle become unstable? At μ=3 and x∗=2/3, we have
f′μ(x∗)=μ(1−2x∗)=−1,
so that the period-1 cycle becomes unstable when the derivative of the map function attains the value of −1, and it is reasonable to expect that the period-2 cycle also becomes unstable when
g′μ(x0)=−1,g′μ(x1)=−1.
Now,
g′μ(x0)=f′μ(fμ(x0))f′μ(x0)=f′μ(x1)f′μ(x0),
and
g′μ(x1)=f′μ(fμ(x1))f′μ(x1)=f′μ(x0)f′μ(x1).
The two equations of (13.26) are therefore identical, and the period-2 cycle will go unstable at the value of μ satisfying
f′μ(x0)f′μ(x1)=−1
where x0 and x1 are given by (13.25).
The equation for μ, then, is given by
μ2(1−2x0)(1−2x1)=−1
Or
μ2(1−2(x0+x1)+4x0x1)+1=0
Now, from (13.25),
x0+x1=μ+1μ,x0x1=μ+1μ2.
Substitution of (13.30) into (13.29) results, after simplification, in the quadratic equation
μ2−2μ−5=0,
with the only positive solution given by
μ=1+√6≈3.449490.
We therefore expect the period-2 cycle to bifurcate to a period-4 cycle at μ≈3.449490.
CHAPTER 13. PENDULUM DYNAMICS
n | μn |
---|---|
1 | 3 |
2 | 3.449490… |
3 | 3.544090… |
4 | 3.564407… |
5 | 3.568759… |
6 | 3.569692… |
7 | 3.569891… |
8 | 3.569934… |
Table 13.1: The first eight values of μn at which period-doubling bifurcations occur.

If we define μn to be the value of μ at which the period-2 n−1 cycle bifurcates to a period-2 n cycle, then we have determined analytically that μ1=3 and μ2=1+√6. We list in Table 13.1, the first eight approximate values of μn, computed numerically.
We can compute a bifurcation diagram for the logistic map. For 2<μ<4, we plot the iterates from the map, discarding initial transients. The bifurcation diagram is shown in Fig. 13.10. Notice the uncanny similarity between the bifurcation diagram for the logistic map, and the one we have previously computed for the damped, driven pendulum equation, Fig. 13.9. Also note that the computation of Fig. 13.10, being on the order of seconds, is substantially faster than that of Fig. 13.9, which took about an hour, because a one-dimensional map is much faster to compute than the Poincaré section of a pair of coupled first-order differential equations.
Computation of the Feigenbaum constant
Period doubling in the logistic map enables an accurate computation of the Feigenbaum constant δ, defined as
δ=limn→∞δn
where
δn=μn+1−μnμn+2−μn+1
Table 13.1 lists the known first eight values of μn at the bifurcation points. These values, and those at even larger values of n, are in fact very difficult to compute with high precision because of the slow convergence of the iterates at the bifurcation points. Rather, we will instead compute the values of μ at what are called superstable cycles. This now well-known method for computing δ was first described by Keith Briggs (1989).
Recall that for the general one-dimensional map
xn+1=fμ(xn),
a perturbation ϵn to a fixed point x∗ decays as
ϵn+1=f′μ(x∗)ϵn.
At a so-called superstable fixed point, however, f′μ(x∗)=0 and the perturbation decays very much faster as
ϵn+1=12f′′μ(x∗)ϵ2n.
What are the superstable fixed points of the logistic map? Now, the values xi that are in the orbit of a period −2n cycle are fixed points of the composite map
xn+1=gμ(xn),
where gu=fu∘fu∘⋯∘fu, where the composition is repeated 2n times. The orbit of a superstable cycle, then, consists of superstable fixed points of (13.37). If the period- 2n cycle has orbit x0,x1,…,x2n−1, we have fμ(x0)=x1,fμ(x1)=x2,…,fμ(x2n−1)=x0, and by the chain rule,
g′μ(xi)=f′μ(x0)f′μ(x1)⋯f′μ(x2n−1),
for all xi in the orbit of the period- 2n cycle. With
f′μ(x)=μ(1−2x)
we have g′μ(x)=0 for x=1/2. Therefore, if x0=1/2, say, is in the orbit of a period- 2n cycle, then this cycle is superstable. At the bifurcation point creating a period-2 n cycle, the cycle has marginal stability and g′μ(x0)=1. As μ increases, g′μ(x0) decreases, and eventually the period-2 n cycle loses stability when g′μ(x0)=−1. At some intermediate value of μ, then, there exists a value of x0 with g′μ(x0)=0, and here we may assign x0=1/2. Therefore, every period- 2n cycle contains a value of μ for which x0=1/2 is in the orbit of the cycle.
Accordingly, we define mn to be the value of μ at which x0=1/2 is in the orbit of the period- 2n cycle. We can modify the definition of the Feigenbaum constant (13.33) to be
δn=mn+1−mnmn+2−mn+1.
Though the values of δn computed from (13.33) and (13.40) will differ slightly, the values as n→∞ should be the same.
We can easily determine the first two values of mn. For the period-1 cycle, we have
12=m012(1−12)
or m0=2, as confirmed from Fig. 13.10. To determine m1, we make use of the period-2 cycle given by (13.21), and Fig. (13.10), which shows that the smaller root passes through x0=1/2. Therefore,
12=12m1((m1+1)−√(m1+1)(m1−3));
and solving for m1, we obtain the quadratic equation
m21−2m1−4=0
with solution m1=1+√5≈3.2361. Further values of mn will be computed numerically.
To determine mn, we need to solve the equation G(μ)=0, where
G(μ)=gμ(1/2)−12
and where as before, gμ(x) is the composition of fμ(x) repeated 2n times. The roots of (13.44) are given by m0,m1,…mn, so that the desired root is the largest one.
We shall use Newton’s method, §2.2, to solve (13.44). To implement Newton’s method, we need to compute both G(μ) and G′(μ). Define N=2n. Then using the logistic map
xn+1=μxn(1−xn),
and iterating with x0=1/2, we obtain x1,x2,…,xN. This orbit is superstable if xN=1/2. Therefore, we have
G(μ)=xN−1/2
Moreover,
G′(μ)=x′N′
where the derivative is with respect to μ. From (13.45), we have
x′n+1=xn(1−xn)+μx′n(1−xn)−μxnx′n=xn(1−xn)+μx′n(1−2xn).
Since we always choose x0=1/2, independent of μ, we have as a starting value x′0=0. Therefore, to compute both xN and x′N′ we iterate 2n times the coupled map equations
xn+1=μxn(1−xn)x′n+1=xn(1−xn)+μx′n(1−2xn)
n | mn | δn |
---|---|---|
0 | 2 | 4.7089430135405 |
1 | 1+√5 | 4.6807709980107 |
2 | 3.4985616993277 | 4.6629596111141 |
3 | 3.5546408627688 | 4.6684039259164 |
4 | 3.5666673798563 | 4.6689537409802 |
5 | 3.5692435316371 | 4.6691571813703 |
6 | 3.5697952937499 | 4.6691910014915 |
7 | 3.5699134654223 | 4.6691994819801 |
8 | 3.5699387742333 | 4.6692010884670 |
9 | 3.5699441946081 | 4.6692015881423 |
10 | 3.5699453554865 | 4.6692023902759 |
11 | 3.5699456041111 | 4.6691974782669 |
12 | 3.5699456573588 | 4.6693329633696 |
13 | 3.5699456687629 | |
14 | 3.5699456712052 |
Table 13.2: The first fourteen values of mn, and estimates of the Feigenbaum delta.
with initial values x0=1/2 and x′0=0. Newton’s method then solves for mn by iterating
μ(i+1)=μ(i)−xN−1/2x′N,
until convergence of μ(i) to mn. In double precision, we have been able to achieve a precision of about 14 digits, which we find can be obtained in fewer than 5 iterations of Newton’s method.
For Newton’s method to converge to mn, we need a good guess for μ(0). We can use the previous best estimate for the Feigenbaum delta to predict mn. From (13.40), we find
μ(0)=mn−1+mn−1−mn−2δn−2.
Although we can not compute δn−2 without knowing mn, we can nevertheless use the estimate δn−2≈δn−3. The computation, then, starts with n=2, and we can begin by taking δ−1=4.4, so that, for example,
μ(0)=3.2361+3.2361−24.4=3.5170.
Using this algorithm, we have produced Table 13.2 for mn, with corresponding calculations of δn. As the values of mn converge, the corresponding values of δn begin to lose precision. It would appear that our best estimate from the table is δ≈4.66920, compared to the known value of δ=4.669201609102990…, computed by a different algorithm capable of achieving higher precision.
Strange attractor of the chaotic pendulum
After the period-doubling cascade, the pendulum motion becomes chaotic. We choose parameter values q=4,f=1.5, and ω=2/3 in the chaotic regime, and after discarding an initial transient of 256 forcing periods, compute a Poincaré section of the phase-space trajectory in the θ−˙θ plane, sampling points every forcing period. The values of the periodic variable θ are mapped onto the interval −π<θ<π. The full Poincaré section consisting of 50,000 points is shown in the top drawing of Fig. 13.11, and a blowup of the points within the drawn rectangle (from a sample of 200,000 points over the entire attractor) is shown in the bottom drawing.

From Fig. 13.11, it appears that the Poincaré section has structure on all scales, which is reminiscent of the classical fractals discussed in §12.7. The set of points shown in Fig. 13.11 is called a strange attractor, and will be seen to have a fractional correlation dimension.
Using the algorithm for computing a correlation dimension discussed in §12.7, we draw a

log-log plot of the correlation integral C(r) versus r, shown in Fig. 13.12. A least-squares fit of a straight line to the middle region of the plot yields a correlation dimension for the Poincaré section of approximately D=1.25.