This semester I’m teaching our Analytical Mechanics course, and I just finished a day on Introduction to Chaos. In this class, we are using the SageMath language to do some simple numerical analysis, so I used Sage to demonstrate some of the interesting behavior you find in chaotic systems. Since I’m learning Sage myself, I thought I would post the result of that class here, to demonstrate some of the Codes and the kinds of plots you can get with simple tools like Sage.

Before getting into the science, SageMath is a free, open-source mathematics software which includes things like Maxima, Python, and the GSL. It’s great because it’s extremely powerful and can be used right in a web browser, thanks for the Sage Cell Server. So I did all of this right in front of my students, to demonstrate how easy this tool is to use.

For the scientific background, I am going to do the same example of the driven, damped pendulum found in *Classical Mechanics* by John Taylor (although the exact same system can be found in *Analytic Mechanics*, by Hand and Finch). So, I didn’t create any of this science, I’m just demonstrating how to study it using Sage.

First, some very basic background. The equation of motion for a driven, damped pendulum of length and mass being acted upon by a driving force is

here is the damping term and , the ratio of the forcing amplitude to the weight of the pendulum. In order to get this into Sage, I’m going to rewrite it as a system of first-order linear differential equations,

This is a typical trick to use numerical integrators, basically because it’s easy to integrate first-order equations, even if they are nonlinear.

It’s easiest to find chaos right near resonance, so let’s pick the parameters and . This means the -axis will display in units of the period, 1 s. We also take . The first plot will be this system when the driving force is the same as the weight. That is, , and code + result is Figure 1 shown below.

`from sage.calculus.desolvers import desolve_system_rk4`

x,y,t=var('x y t')

w=2*pi

w0=3*pi

g=3/4*pi

f=1.0

P=desolve_system_rk4([-2*g*x-w0^2*sin(y)+f*w0^2*cos(w*t),x],[x,y],[0,0,0],ivar=t,end_points=[0,15],step=0.01)

Q=[[i,k] for i,j,k in P]

intP=spline(Q)

plot(intP,0,15)

Figure 2 is a plot with the driving force slightly bigger than the weight, .

This demonstrates an *attractor*, meaning the steady-state solution eventually settles down to oscillate around . We can check this is actually still periodic by asking Sage for the value of at =30 s, =31 s, etc., by calling this line instead of the plot command above

`[intP(i) for i in range(30,40)]`

(Note that we also have to change the range of integration from to .) The output is shown in Figure 3; the period is clearly 1.0 s out to four significant figures.

Next, let’s increase the forcing to . The result is shown in Figure 7. The attractor is still present (now with a value of around ), but the behavior is much more dramatic. In fact, you might not even be convinced that the period is still 1.0 s, since the peaks look to be different values. We can repeat our experiment from above, and ask Sage to print out the value of for integer timesteps between =30 and =40. The result is shown in Figure 4. The actual period appears to be 2.0 s, since the value of does not repeat exactly after 1.0 s. This is called *Period Doubling*.

In Figure 8, I’ve displayed a plot with , and it’s immediately obvious that the oscillatory motion now has period 3.0 s. We can check this by playing the same game, shown in Figure 6.

Now we are in a position to see some unique behavior. I am going to overlay a new solution onto this one, but give the second solution a different initial value, instead of . The code I am adding is

`P2=desolve_system_rk4([-2*b*x-w0^2*sin(y)+g*w0^2*cos(w*t),x],[x,y],[0,0,-pi/2],`

ivar=t,end_points=[0,15],step=0.01)

Q2=[[i,k] for i,j,k in P2]

intP2=spline(Q2)

plot(intP,0,15)+plot(intP2,0,15,linestyle=":", color=''red'')

The result is shown in Figure 8. Here we can see the first example of the sensitivity to initial conditions. The two solutions diverge markedly once you have a slightly different initial condition, heading towards two very different attractors. Let’s plot the difference between the two oscillators,

but with only a very small difference in the initial conditions, . The code follows:

`#plot(intP,0,15)+plot(intP2,0,15,linestyle=":", color=''red'')`

plot(lambda x: abs(intP(x)-intP2(x)),0,15)

This is shown in Figure 9. It clearly decays to zero, but that’s hard to see so let’s plot it on a log scale, shown in Figure 10.

`#plot(intP,0,15)+plot(intP2,0,15,linestyle=":", color=''red'')`

plot_semilogy(lambda x: abs(intP(x)-intP2(x)),0,15)

Now, let’s see what happens if we do this same thing, but make the force parameter over the critical value of . This is displayed in Figure 11, for . We get completely the *opposite* behavior, the differences in the oscillators are driven *away* from each other due to their small initial separations. This is the essence of “Jurrasic Park Chaos” – a small change in the initial conditions (like a butterfly flapping it’s wings in Malaysa) causes a large change in the final outcome (a change in the weather pattern over California).