A Study of the Damped, Driven Pendulum 

  following John Taylor's Classical Mechanics 

 

We are interested in the behavior of this pendulum, damped, and driven sinusoidally. 

The equation is given by him as  

                        

To make it easier to type, I have used below  x  for φ, b for β, and g for γ. 

 

The first task is to initialize the equation, i.e., to tell Maple the desired initial conditions for the differential equation. 

> `:=`(ic, (D(x))(0) = 0., x(0) = 0.); 1
 

(D(x))(0) = 0., x(0) = 0. (1)
 

Next we must give the desired values for the parameters that define the problem.   

All the ones given below are simply those used by Taylor, as defined by him on p. 464. 

> `:=`(omega, `+`(`*`(2, `*`(evalf(Pi))))); 1; `:=`(omega0, `+`(`*`(1.5, `*`(omega)))); 1; `:=`(b, `+`(`*`(`/`(1, 4), `*`(omega0)))); 1; `:=`(g, .2); 1
 

 

 

 

6.283185308
9.424777962
2.356194490
.2 (2)
 

and then the differential equation, in Maple input form, and with a name: 

> `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1
 

`+`(diff(diff(x(t), t), t), `*`(4.712388980, `*`(diff(x(t), t))), `*`(88.82643963, `*`(sin(x(t))))) = `+`(`*`(17.76528793, `*`(cos(`+`(`*`(6.283185308, `*`(t))))))) (3)
 

> `:=`(sol, dsolve({de, ic}, x(t), numeric, maxfun = 0)); 1
 

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (4)
 

This uses the Maple command  dsolve  which solves differential equations, and it chooses to uses the fourth-fifth order Runge-Kutta method, 

since this is a fairly straightforward, although nonlinear, differential equation. 

The option maxfun=0 tells it not to stop calculating until it finds the answer that has been requested. 

 

After having acquired the numerical solutions to the differential equation, we will want to plot them; therefore, we  

now bring up some extra plot functions. 

> with(plots); 1
 

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fiel...
[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fiel...
[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fiel...
[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fiel...
[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fiel...
[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fiel...
(5)
 

We may now simply create a plot of the solution for the equation, with these initial conditions, and parameter values. 

As we will see below, after the transient, which is not very long-lived, the solution is nicely (singly) periodic. 

> odeplot(sol, [t, x(t)], 0 .. 12, numpoints = 500, title = `driven, damped pendulum,\n with gamma only = 0.2`); 1
 

Plot_2d
 

The values of phi only vary between about -0.3 and +0.3, so the pendulum is not being driven enough that it can go over the top. 

The graph certainly does look periodic, after about t=1.5.   

However, there are various other ways to ask if this is so. 

We now create 3-dimensional plots of the same data, and view them so that we have a phase-plane plot, i.e., we see v versus x. 

> `:=`(p[1], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 12, numpoints = 500, color = blue)); -1; display([seq(p[i], i = 1 .. 1)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pe...
`:=`(p[1], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 12, numpoints = 500, color = blue)); -1; display([seq(p[i], i = 1 .. 1)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pe...
`:=`(p[1], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 12, numpoints = 500, color = blue)); -1; display([seq(p[i], i = 1 .. 1)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pe...
 

Plot
 

One can see the transient in the above graph, where it begins in the middle, since the initial conditions were (0,0), and is then forced outward, 

and indeed "overshoots" a bit, into the wider circle, before it "calms down" and then moves in the ellipse shown in dark color---because there  

are very many points there. 

This is perhaps more clear if we look only at those times after the transient has decayed, as we do below: 

> `:=`(pt[1], odeplot(sol, [[t, x(t), diff(x(t), t)]], 2 .. 12, numpoints = 500, color = blue)); -1; display([seq(pt[i], i = 1 .. 1)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `...
`:=`(pt[1], odeplot(sol, [[t, x(t), diff(x(t), t)]], 2 .. 12, numpoints = 500, color = blue)); -1; display([seq(pt[i], i = 1 .. 1)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `...
`:=`(pt[1], odeplot(sol, [[t, x(t), diff(x(t), t)]], 2 .. 12, numpoints = 500, color = blue)); -1; display([seq(pt[i], i = 1 .. 1)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `...
 

Plot
 

This is actually a 3-dimensional graph, which we are looking at head-on. 

We can see this for sure by taking out the orientation command, when we will see that it is actually a cylinder. 

 

> display([seq(pt[i], i = 1 .. 1)], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=0.2, without transient`, style = point, symbol = cross, symbolsize = 10); 1
display([seq(pt[i], i = 1 .. 1)], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=0.2, without transient`, style = point, symbol = cross, symbolsize = 10); 1
 

Plot
 

We now want to increase---noticeably---the value of the forcing parameter. 

This will correspond to Taylor's Fig. 12.3, on p. 466. 

> `:=`(g, .9); 1
 

.9 (6)
 

Having changed the value of the parameter, g, it is "safest" to re-define the equation for the differential equation, and its solution: 

Note that we have NOT changed the initial conditions. 

> `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, dsolve({de, ic}, x(t)...
 

 

`+`(diff(diff(x(t), t), t), `*`(4.712388980, `*`(diff(x(t), t))), `*`(88.82643963, `*`(sin(x(t))))) = `+`(`*`(79.94379567, `*`(cos(`+`(`*`(6.283185308, `*`(t)))))))
proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (7)
 

As before we now re-create his graph.  Again it has just one period, although the transient takes a little longer to disappear. 

Notice that the value of the angle is maximum slightly less than 2.5, so that the pendulum is still not able to swing over the top. 

> odeplot(sol, [t, x(t)], 0 .. 12, numpoints = 1000); 1
 

Plot_2d
 

Again we also display the state-space plots. 

> `:=`(p[2], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 12, numpoints = 500, color = black)); -1; display([seq(p[i], i = 2 .. 2)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `p...
`:=`(p[2], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 12, numpoints = 500, color = black)); -1; display([seq(p[i], i = 2 .. 2)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `p...
`:=`(p[2], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 12, numpoints = 500, color = black)); -1; display([seq(p[i], i = 2 .. 2)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `p...
 

Plot
 

> `:=`(pt[2], odeplot(sol, [[t, x(t), diff(x(t), t)]], 2 .. 12, numpoints = 500, color = black)); -1; display([seq(pt[i], i = 2 .. 2)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = ...
`:=`(pt[2], odeplot(sol, [[t, x(t), diff(x(t), t)]], 2 .. 12, numpoints = 500, color = black)); -1; display([seq(pt[i], i = 2 .. 2)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = ...
`:=`(pt[2], odeplot(sol, [[t, x(t), diff(x(t), t)]], 2 .. 12, numpoints = 500, color = black)); -1; display([seq(pt[i], i = 2 .. 2)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = ...
 

Plot
 

Continuing onward, we now change the forcing parameter so that it is stronger than gravity, i.e., g is now greater than 1. 

However, we will skip Taylor's graph of g=1.06, and go onward to his g=1.073 [Figure 12.5, p. 469] 

 

g:=1.073 

 

> `:=`(g, 1.073); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
`:=`(g, 1.073); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
 

 

 

1.073
`+`(diff(diff(x(t), t), t), `*`(4.712388980, `*`(diff(x(t), t))), `*`(88.82643963, `*`(sin(x(t))))) = `+`(`*`(95.31076972, `*`(cos(`+`(`*`(6.283185308, `*`(t)))))))
proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (8)
 

> odeplot(sol, [t, x(t)], 0 .. 40, numpoints = 1000); 1
 

Plot_2d
 

At this point in Taylor he shows an enlarged view of the bottoms of these oscillations, to show that there truly are two periods. 

We do the same thing below. 

The two periods are really very obvious.  Thinking about why we can't really see that in the graph above, I presume that they simply didn't have enough points chosen. 

> sol(25), sol(25.5), sol(26); 1
 

[t = 25., x(t) = -8.64041673555413681, diff(x(t), t) = 11.6944671377156198], [t = 25.5, x(t) = -6.06232712563896480, diff(x(t), t) = -18.4978598627688342], [t = 26., x(t) = -8.71660149782448279, diff(...
[t = 25., x(t) = -8.64041673555413681, diff(x(t), t) = 11.6944671377156198], [t = 25.5, x(t) = -6.06232712563896480, diff(x(t), t) = -18.4978598627688342], [t = 26., x(t) = -8.71660149782448279, diff(...
(9)
 

> sol(25.25); 1
 

[t = 25.25, x(t) = -4.48444212210043780, diff(x(t), t) = 8.39258812523918784] (10)
 

> odeplot(sol, [t, x(t)], 25 .. 40, numpoints = 1000, view = [25 .. 40, -9.5 .. -9.0]); 1
 

Plot_2d
 

Notice, as Taylor remarks, that this system now does rather complicated looking things before it finally settles down 

to something that is reasonably well-behaved.  We will see however, especially below, that there are now two periods. 

In addition, our pendulum is now able to swing "over the top," so that the maximum of the initial upward swings is more  

than 14, which is more than 4 pi, i.e., more than two complete swings over the top, each time coming back and forth around. 

The final periodic situation is down around -7, i.e., it is swinging back and forth steadily, but in a place that is one period  

backward from its initial position. 

 

The phase-plane plots show a lot of this erratic behavior, and then, in the one that omits the transient, 

shows clearly the two different periods. 

> `:=`(p[3], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 50, numpoints = 1000, color = green)); -1; `:=`(pt[3], odeplot(sol, [[t, x(t), diff(x(t), t)]], 25 .. 50, numpoints = 500, color = green)); -1;...
`:=`(p[3], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 50, numpoints = 1000, color = green)); -1; `:=`(pt[3], odeplot(sol, [[t, x(t), diff(x(t), t)]], 25 .. 50, numpoints = 500, color = green)); -1;...
`:=`(p[3], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 50, numpoints = 1000, color = green)); -1; `:=`(pt[3], odeplot(sol, [[t, x(t), diff(x(t), t)]], 25 .. 50, numpoints = 500, color = green)); -1;...
 

> display([seq(p[i], i = 3 .. 3)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=1.073`, style = point, symbol = cross, symbolsize = 10); 1; display(pw[...
display([seq(p[i], i = 3 .. 3)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=1.073`, style = point, symbol = cross, symbolsize = 10); 1; display(pw[...
display([seq(p[i], i = 3 .. 3)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=1.073`, style = point, symbol = cross, symbolsize = 10); 1; display(pw[...
display([seq(p[i], i = 3 .. 3)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=1.073`, style = point, symbol = cross, symbolsize = 10); 1; display(pw[...
display([seq(p[i], i = 3 .. 3)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=1.073`, style = point, symbol = cross, symbolsize = 10); 1; display(pw[...
 

 

 

Plot
Plot
Plot
 

So we see that letting Maple expand the scaling allows us to see the two periods better. 

Now we continue to follow Taylor, and increase the forcing strength a bit more, which should give us period 3. 

g:= 1.077 

 

> `:=`(g, 1.077); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
`:=`(g, 1.077); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
 

 

 

1.077
`+`(diff(diff(x(t), t), t), `*`(4.712388980, `*`(diff(x(t), t))), `*`(88.82643963, `*`(sin(x(t))))) = `+`(`*`(95.66607548, `*`(cos(`+`(`*`(6.283185308, `*`(t)))))))
proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (11)
 

> odeplot(sol, [t, x(t)], 0 .. 20, numpoints = 1000); 1
 

Plot_2d
 

This time we see that there is indeed a periodic structure; 

however, it amounts to oscillations between almost 15 and about 4, i.e, between about 5Pi and 1.3Pi, so that  it is regularly swinging 

all the way over the top, changing its mind and turning around and going back over the top, etc. 

 

In the next command I am simply saving this graph so that I can compare it later with a different set of initital conditions. 

> `:=`(ic1g1077, odeplot(sol, [t, x(t)], 0 .. 20, numpoints = 1000)); -1
 

The next thing Taylor does is to try to prove seriously that it really does have 3 periods.   

Remembering that the basic period of the driver is  T = 1, i.e., omega = 2*Pi, he evaluates the solution at integer values of t, beginning with t=30, as we do below, showing the same effect. 

> seq(rhs(sol(`+`(n, 30))[2]), n = 0 .. 12); 1
 

13.8122290832031000, 7.75847279777866384, 6.87274040053776768, 13.8122289756238016, 7.75847260376588466, 6.87274071353994920, 13.8122289316491412, 7.75847251933632798, 6.87274085290107984, 13.81222891...
13.8122290832031000, 7.75847279777866384, 6.87274040053776768, 13.8122289756238016, 7.75847260376588466, 6.87274071353994920, 13.8122289316491412, 7.75847251933632798, 6.87274085290107984, 13.81222891...
13.8122290832031000, 7.75847279777866384, 6.87274040053776768, 13.8122289756238016, 7.75847260376588466, 6.87274071353994920, 13.8122289316491412, 7.75847251933632798, 6.87274085290107984, 13.81222891...
(12)
 

> Digits; 1
 

10 (13)
 

The graph below is, perhaps, a better, or perhaps a worse, way to see the same phenomenon! 

> odeplot(sol, [t, x(t)], 30 .. 40, numpoints = 1000, view = [30 .. 40, 6 .. 15]); 1
 

Plot_2d
 

Probably a better way to see that, graphically, is with our first instance of a Poincare plot, below. 

There are two, somewhat different ways to view it.   

In the first graph below I simply plot the values of phi = x, themselves 

> pointplot([seq(rhs(sol(`+`(n, 30))[2]), n = 0 .. 39)], color = red); 1
 

Plot_2d
 

However, in this graph I plot things in the real phase-plane mode, where I plot phi horizongtally and phidot vertically! 

> pointplot([seq([rhs(sol(`+`(n, 30))[2]), rhs(sol(`+`(n, 30))[3])], n = 0 .. 39)], view = [0 .. 15, -4 .. 18], color = red); 1
 

Plot_2d
 

Now, however, before doing that, let's also look at the phase-plane plots for this situation, 

where we will see that, yes, there was a long transient behavior, but there are also two distinct periods. 

> `:=`(p[4], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = brown)); -1; `:=`(pt[4], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = brown)); -1; ...
`:=`(p[4], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = brown)); -1; `:=`(pt[4], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = brown)); -1; ...
`:=`(p[4], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = brown)); -1; `:=`(pt[4], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = brown)); -1; ...
`:=`(p[4], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = brown)); -1; `:=`(pt[4], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = brown)); -1; ...
`:=`(p[4], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = brown)); -1; `:=`(pt[4], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = brown)); -1; ...
 

 

Plot
Plot
 

If you're quite careful, you can, perhaps, make out three different ovals (differently-scaled ellipses) that constitute the pattern just above. 

 

Now, Taylor wants to look at this very same value of the forcing strength, but with a fairly-different initial condition. 

The earlier initial condition was to start at rest while hanging straight down, i.e., at the bottom; 

this new initial condition starts out, from rest, hanging straight out, along the horizontal. 

We will find that the two generate very different results. 

> `:=`(ic2, (D(x))(0) = 0., x(0) = `+`(`-`(`*`(`/`(1, 2), `*`(evalf(Pi)))))); 1
 

(D(x))(0) = 0., x(0) = -1.570796327 (14)
 

We don't change the differential equation this time, because all we're doing is changing the initial conditions. 

> `:=`(sol, dsolve({de, ic2}, x(t), numeric, maxfun = 0)); 1
 

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (15)
 

> odeplot(sol, [t, x(t)], 0 .. 20, numpoints = 1000, color = green); 1
 

Plot_2d
 

Now, I'm going to save that graph as well, and then try to look more at the structure, in detail. 

Taylor claims that with these initial conditions, this same forcing is ATTRACTED to a two-period attractor. 

Let's try to actually see that! 

> `:=`(ic2g1077, odeplot(sol, [t, x(t)], 0 .. 20, numpoints = 1000, color = green)); -1
 

Once again this is what the ode solver is giving us as solutions. 

> sol(30); 1
 

[t = 30., x(t) = -.107086698453751932, diff(x(t), t) = 17.8941059651618026] (16)
 

So we try below to pull out 40 repeats of the driving period.   

> pointplot([seq(rhs(sol(`+`(n, 30))[2]), n = 0 .. 39)], color = green); 1
 

Plot_2d
 

However, that doesn't look very good.  BUT, perhaps the problem is just with some minor sort of computer non-round-off,  

as we might guess by looking at the tickmarks in the margins. 

Therefore, we next look a little bit more at the solutions themselves. 

> sol(30), sol(40), sol(41); 1
 

[t = 30., x(t) = -.107086698453751932, diff(x(t), t) = 17.8941059651618026], [t = 40., x(t) = -.107081061641198794, diff(x(t), t) = 17.8940581672670618], [t = 41., x(t) = -.405005533642967320, diff(x(...
[t = 30., x(t) = -.107086698453751932, diff(x(t), t) = 17.8941059651618026], [t = 40., x(t) = -.107081061641198794, diff(x(t), t) = 17.8940581672670618], [t = 41., x(t) = -.405005533642967320, diff(x(...
(17)
 

And then we plot the phase-plane points for 40 periods, well after the transients have stopped.  At this point we see that we do have just TWO periods. 

> pointplot([seq([rhs(sol(`+`(n, 30))[2]), rhs(sol(`+`(n, 30))[3])], n = 0 .. 39)], view = [-.42 .. -.10, 17.8 .. 19], color = green); 1
 

Plot_2d
 

>
 

Remembering the previous situation, where it oscillated up and over all the time, this is also periodic, but the oscillations are only  

between Pi and -Pi, i.e., all the way to the top and back down and to the top from the other direction. 

 

To understand that better, we now plot the two graphs together, which was the purpose of saving them, earlier. 

Note, on the graph below, that they do actually start out in almost the same way, but then diverge sharply,  

each staying in its own periodic mode. 

In this view it is not obvious, but in fact each of them is doubly periodic, as is visible in the phase-plane views. 

> display([ic1g1077, ic2g1077], title = `comparison of 2 initial conditions`); 1
 

Plot_2d
 

Below we will first look at the phase-plane plots for the  

> `:=`(p[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = orange)); -1; `:=`(pt[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = blue)); -1; ...
`:=`(p[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = orange)); -1; `:=`(pt[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = blue)); -1; ...
`:=`(p[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = orange)); -1; `:=`(pt[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = blue)); -1; ...
`:=`(p[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = orange)); -1; `:=`(pt[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = blue)); -1; ...
`:=`(p[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = orange)); -1; `:=`(pt[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = blue)); -1; ...
`:=`(p[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 25, numpoints = 1000, color = orange)); -1; `:=`(pt[5], odeplot(sol, [[t, x(t), diff(x(t), t)]], 5 .. 25, numpoints = 500, color = blue)); -1; ...
 

 

Plot
Plot
 

> display([seq(pt[i], i = 4 .. 5)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=1.077, \n comparing 2  sets of  initial conditions`, style = point, sy...
display([seq(pt[i], i = 4 .. 5)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=1.077, \n comparing 2  sets of  initial conditions`, style = point, sy...
 

Plot
 

In the comparison of the phase-plane views, we see that this change of the initial condition---both from rest, but at different initial locations---caused 

really quite different behaviors. 

 

Now, let's go onward, to a rather larger value for the forcing strength: 

g:=1.081 

We are looking now at the lst two of the sequence of figures in Figure 12.8, on p. 472. 

> `:=`(g, 1.081); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
`:=`(g, 1.081); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
 

 

 

1.081
`+`(diff(diff(x(t), t), t), `*`(4.712388980, `*`(diff(x(t), t))), `*`(88.82643963, `*`(sin(x(t))))) = `+`(`*`(96.02138124, `*`(cos(`+`(`*`(6.283185308, `*`(t)))))))
proc (x_rkf45) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; `:=`(_EnvDSNumericSaveDigits, Digits); `:=`(Digits, 14); if... (18)
 

> odeplot(sol, [t, x(t)], 28 .. 40, numpoints = 1000, color = black); 1
 

Plot_2d
 

Below we narrow in on the bottom, where we can indeed see 4 distinct periods. 

> odeplot(sol, [t, x(t)], 28 .. 40, numpoints = 1000, color = black, view = [30 .. 42, -2.4 .. -1.3]); 1
 

Plot_2d
 

> `:=`(pt[6], odeplot(sol, [[t, x(t), diff(x(t), t)]], 28 .. 40, numpoints = 500, color = blue)); -1; display([seq(pt[i], i = 6 .. 6)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = ...
`:=`(pt[6], odeplot(sol, [[t, x(t), diff(x(t), t)]], 28 .. 40, numpoints = 500, color = blue)); -1; display([seq(pt[i], i = 6 .. 6)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = ...
`:=`(pt[6], odeplot(sol, [[t, x(t), diff(x(t), t)]], 28 .. 40, numpoints = 500, color = blue)); -1; display([seq(pt[i], i = 6 .. 6)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = ...
 

Plot
 

We see that this one has 4 distinct periods! 

 

We then go further, and increase the strength a bit more: 

g:=1.0826,  

keeping the latest initial conditions, 

therefore repeating the last of the figures in that set, on p. 472. 

> `:=`(g, 1.0826); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, d...
`:=`(g, 1.0826); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, d...
 

 

 

1.0826
`+`(diff(diff(x(t), t), t), `*`(4.712388980, `*`(diff(x(t), t))), `*`(88.82643963, `*`(sin(x(t))))) = `+`(`*`(96.16350354, `*`(cos(`+`(`*`(6.283185308, `*`(t)))))))
proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (19)
 

> odeplot(sol, [t, x(t)], 28 .. 40, numpoints = 1000, color = blue); 1
 

Plot_2d
 

> odeplot(sol, [t, x(t)], 28 .. 56, numpoints = 2000, color = blue, view = [28 .. 56, -2.5 .. -1.5]); 1
 

Plot_2d
 

> `:=`(pt[7], odeplot(sol, [[t, x(t), diff(x(t), t)]], 28 .. 80, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 7 .. 7)], orientation = [0, 90], axes = framed, labels = [t, x, v], title =...
`:=`(pt[7], odeplot(sol, [[t, x(t), diff(x(t), t)]], 28 .. 80, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 7 .. 7)], orientation = [0, 90], axes = framed, labels = [t, x, v], title =...
`:=`(pt[7], odeplot(sol, [[t, x(t), diff(x(t), t)]], 28 .. 80, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 7 .. 7)], orientation = [0, 90], axes = framed, labels = [t, x, v], title =...
 

Plot
 

However, the Poincare plot approach does indeed show all 8 of them. 

> pointplot([seq([rhs(sol(`+`(n, 30))[2]), rhs(sol(`+`(n, 30))[3])], n = 0 .. 40)], view = [-.6 .. 0, 17.5 .. 19.5], color = blue); 1
 

Plot_2d
 

We see quite a few different periods in there---Taylor claims 8, although I'm not sure I see that many. 

 

Next, we go to  

g:=1.105, 

which is larger than the value of 1.0829 where chaos sets in. 

 

> `:=`(g, 1.105); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
`:=`(g, 1.105); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
 

 

 

1.105
`+`(diff(diff(x(t), t), t), `*`(4.712388980, `*`(diff(x(t), t))), `*`(88.82643963, `*`(sin(x(t))))) = `+`(`*`(98.15321579, `*`(cos(`+`(`*`(6.283185308, `*`(t)))))))
proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (20)
 

> odeplot(sol, [t, x(t)], 8 .. 40, numpoints = 1000, color = red); 1
 

Plot_2d
 

> `:=`(pt[8], odeplot(sol, [[t, x(t), diff(x(t), t)]], 14 .. 21, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 8 .. 8)], orientation = [0, 90], axes = framed, labels = [t, x, v], title =...
`:=`(pt[8], odeplot(sol, [[t, x(t), diff(x(t), t)]], 14 .. 21, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 8 .. 8)], orientation = [0, 90], axes = framed, labels = [t, x, v], title =...
`:=`(pt[8], odeplot(sol, [[t, x(t), diff(x(t), t)]], 14 .. 21, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 8 .. 8)], orientation = [0, 90], axes = framed, labels = [t, x, v], title =...
 

Plot
 

One can see the different places the pendulum is going pretty easily here, since we have only gone through 7 periods of the forcing system. 

 

Below we will go through 192 periods, and it gets closer to filling in a decent-sized region of phase space. 

> `:=`(pt[9], odeplot(sol, [[t, x(t), diff(x(t), t)]], 8 .. 200, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 9 .. 9)], orientation = [0, 90], axes = framed, labels = [t, x, v], title =...
`:=`(pt[9], odeplot(sol, [[t, x(t), diff(x(t), t)]], 8 .. 200, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 9 .. 9)], orientation = [0, 90], axes = framed, labels = [t, x, v], title =...
`:=`(pt[9], odeplot(sol, [[t, x(t), diff(x(t), t)]], 8 .. 200, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 9 .. 9)], orientation = [0, 90], axes = framed, labels = [t, x, v], title =...
 

Plot
 

>
 

> pointplot([seq([rhs(sol(`+`(n, 20))[2]), rhs(sol(`+`(n, 20))[3])], n = 0 .. 40)], color = blue, view = [-3.4 .. 2.7, -22 .. 22]); 1
 

Plot_2d
 

Let's try doing quite a few more 

> pointplot([seq([rhs(sol(`+`(n, 20))[2]), rhs(sol(`+`(n, 20))[3])], n = 0 .. 400)], color = blue, view = [-3.4 .. 2.7, -22 .. 22]); 1
 

Plot_2d
 

Yet one more time, increasing the forcing strength,  

this time to  

g:=1.130. 

 

> `:=`(g, 1.13); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, dso...
`:=`(g, 1.13); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, dso...
 

 

 

1.13
`+`(diff(diff(x(t), t), t), `*`(4.712388980, `*`(diff(x(t), t))), `*`(88.82643963, `*`(sin(x(t))))) = `+`(`*`(100.3738768, `*`(cos(`+`(`*`(6.283185308, `*`(t)))))))
proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (21)
 

> odeplot(sol, [t, x(t)], 8 .. 40, numpoints = 1000, color = red); 1
 

Plot_2d
 

> odeplot(sol, [t, x(t)], 8 .. 40, numpoints = 4000, color = red, view = [20 .. 40, -3 .. 2]); 1
 

Plot_2d
 

> `:=`(pt[10], odeplot(sol, [[t, x(t), diff(x(t), t)]], 8 .. 40, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 10 .. 10)], orientation = [0, 90], axes = framed, labels = [t, x, v], title...
`:=`(pt[10], odeplot(sol, [[t, x(t), diff(x(t), t)]], 8 .. 40, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 10 .. 10)], orientation = [0, 90], axes = framed, labels = [t, x, v], title...
`:=`(pt[10], odeplot(sol, [[t, x(t), diff(x(t), t)]], 8 .. 40, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 10 .. 10)], orientation = [0, 90], axes = framed, labels = [t, x, v], title...
 

Plot
 

Here we may again put out a phase-plane Poincare plot, and see that there are only 3 periods. 

> pointplot([seq([rhs(sol(`+`(n, 20))[2]), rhs(sol(`+`(n, 20))[3])], n = 0 .. 40)], color = blue, view = [-2.5 .. 8, -20 .. 20]); 1
 

Plot_2d
 

As noted in Taylor, on p. 481, this is not chaotic, even though the forcing is larger than the onset of chaos. 

Instead there seem to be 3 distinct periods. 

As he notes, there are regions of periodicity between regions of chaos. 

 

Therefore, we now go much larger in our forcing strength,  

to g:=1.503. 

 

> `:=`(g, 1.503); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
`:=`(g, 1.503); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
 

 

 

1.503
`+`(diff(diff(x(t), t), t), `*`(4.712388980, `*`(diff(x(t), t))), `*`(88.82643963, `*`(sin(x(t))))) = `+`(`*`(133.5061388, `*`(cos(`+`(`*`(6.283185308, `*`(t)))))))
proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error (22)
 

> odeplot(sol, [t, x(t)], 0 .. 40, numpoints = 1000, color = blue); 1
 

Plot_2d
 

> `:=`(no1icchaotic, odeplot(sol, [t, x(t)], 0 .. 40, numpoints = 1000, color = blue)); -1
 

Notice that this one begins near zero, actually at -Pi/2, and then starts revolving rapidly, and erratically, turning round and round, all the way  

down about -65, i.e., about 20 turns in the same direction, and then bounces around a bit, and then starts turning back around and around  

the other direction, and is erratic. 

 

The erratic nature is seen much better in the phase-plane plot below. 

> `:=`(pt[11], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 40, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 11 .. 11)], orientation = [0, 90], axes = framed, labels = [t, x, v], title...
`:=`(pt[11], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 40, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 11 .. 11)], orientation = [0, 90], axes = framed, labels = [t, x, v], title...
`:=`(pt[11], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 40, numpoints = 1000, color = blue)); -1; display([seq(pt[i], i = 11 .. 11)], orientation = [0, 90], axes = framed, labels = [t, x, v], title...
 

Plot
 

Below we plot the angular velocity, rather than the angle itself. 

> odeplot(sol, [t, diff(x(t), t)], 0 .. 40, numpoints = 1000, color = blue); 1
 

Plot_2d
 

>
 

APoincare phase-plane plot shows almost as many disparate dots as there were period repeats. 

> pointplot([seq([rhs(sol(`+`(n, 20))[2]), rhs(sol(`+`(n, 20))[3])], n = 0 .. 40)], color = blue, view = [-70 .. 0, -25 .. 25]); 1
 

Plot_2d
 

And we can do that in more detail. 

> pointplot([seq([rhs(sol(`+`(n, 20))[2]), rhs(sol(`+`(n, 20))[3])], n = 0 .. 400)], color = blue, view = [-70 .. 0, -25 .. 25]); 1
 

Plot_2d
 

As yet a different way of looking at the problem, since the pendulum is rotating over the top so much, let us, instead, look at how  

its angular velocity is changing over time, which we do above., 

 

However, at this point we want to begin to study something slightly different, namely how do the solutions differ from one another 

when their initial conditions are very close. 

Therefore, I have saved the plot of x(t) above, and now want to create a separate plot for a different initial condition, but  

one which differs only by one thousandth of a radian, as given below. 

 

> `:=`(ic3, (D(x))(0) = 0., x(0) = `+`(`-`(`*`(`/`(1, 2), `*`(evalf(Pi)))), 0.1e-2)); 1
 

(D(x))(0) = 0., x(0) = -1.569796327 (23)
 

> `:=`(sol, dsolve({de, ic3}, x(t), numeric, maxfun = 0)); 1
 

proc (x_rkf45) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; `:=`(_EnvDSNumericSaveDigits, Digits); `:=`(Digits, 14); if... (24)
 

> `:=`(no2icchaotic, odeplot(sol, [t, x(t)], 0 .. 40, numpoints = 1000, color = green)); -1
 

> display([no1icchaotic, no2icchaotic]); 1
 

Plot_2d
 

We see that they start out so close that we cannot tell the difference between them, for the first 10 cycles, and then they diverge  

markedly, and, presumably, forever after. 

On that score, I note that this graph has been run out to 40 cycles, whereas the one Taylor presents is only 25; the divergence  

continues! 

We can also see the difference in the phase plane plots, below. 

> `:=`(pt[12], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 40, numpoints = 1000, color = green)); -1; display([seq(pt[i], i = 12 .. 12)], orientation = [0, 90], axes = framed, labels = [t, x, v], titl...
`:=`(pt[12], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 40, numpoints = 1000, color = green)); -1; display([seq(pt[i], i = 12 .. 12)], orientation = [0, 90], axes = framed, labels = [t, x, v], titl...
`:=`(pt[12], odeplot(sol, [[t, x(t), diff(x(t), t)]], 0 .. 40, numpoints = 1000, color = green)); -1; display([seq(pt[i], i = 12 .. 12)], orientation = [0, 90], axes = framed, labels = [t, x, v], titl...
 

Plot
 

And now we will display the two phase plane plots together, in the two colors: 

> display([seq(pt[i], i = 11 .. 12)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=1.503,\n comparison`, style = point, symbol = cross, symbolsize = 10...
display([seq(pt[i], i = 11 .. 12)], orientation = [0, 90], axes = framed, labels = [t, x, v], title = `pendulum motion,\n for gamma=1.503,\n comparison`, style = point, symbol = cross, symbolsize = 10...
 

Plot
 

It is perhaps worth retreating a bit at this point, since those two graphs diverged so very much when their initial conditions differed only by .001 radian. 

Let us go back and find such a situation when gamma was much smaller. 

 

> `:=`(g, 1.077); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
`:=`(g, 1.077); 1; `:=`(de, `+`(diff(x(t), t, t), `*`(2, `*`(b, `*`(diff(x(t), t)))), `*`(`^`(omega0, 2), `*`(sin(x(t))))) = `*`(g, `*`(`^`(omega0, 2), `*`(cos(`*`(omega, `*`(t))))))); 1; `:=`(sol, ds...
 

 

 

1.077
`+`(diff(diff(x(t), t), t), `*`(4.712388980, `*`(diff(x(t), t))), `*`(88.82643963, `*`(sin(x(t))))) = `+`(`*`(95.66607548, `*`(cos(`+`(`*`(6.283185308, `*`(t)))))))
proc (x_rkf45) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; `:=`(_EnvDSNumericSaveDigits, Digits); `:=`(Digits, 14); if... (25)
 

> `:=`(ic3g1077, odeplot(sol, [t, x(t)], 0 .. 20, numpoints = 1000, color = brown)); -1
 

> display([ic2g1077, ic3g1077]); 1
 

Plot_2d
 

If you look very, very carefully, in the graph above you can flecks of green here and there among the brown lines. 

The green and the brown lines differ in their initial conditions by .001 radian, and you can see they never differ very much. 

It's only evident at all at the tips of the lines; that is probably an artifact of the computer, which of course only samples the  

integration every where, and needs more sampling at the places where it is turning sharply. 

>
 

>