Chaos and SageMath

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 L and mass m being acted upon by a driving force F(t)=F_0\cos\omega t is

\ddot{\phi}+2\gamma\dot{\phi}-\omega_0^2\sin(\phi)=f\omega_0^2\cos(\omega t)

\gamma here is the damping term and f=F_0/(mg), 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,

\dot y=x

\dot x=f\omega_0^2\cos(\omega t)+\omega_0^2\sin(y)-2\gamma x

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 \omega_0=3\pi and \omega=2\pi. This means the t-axis will display in units of the period, 1 s. We also take \gamma=3/4\pi. The first plot will be this system when the driving force is the same as the weight. That is, f=1, and code + result is Figure 1 shown below.

from sage.calculus.desolvers import desolve_system_rk4
x,y,t=var('x y t')
Q=[[i,k] for i,j,k in P]

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

Figure 1, f=1

Figure 2, f=1.06







This demonstrates an attractor, meaning the steady-state solution eventually settles down to oscillate around \phi\approx 2\pi. We can check this is actually still periodic by asking Sage for the value of \phi at t=30 s, t=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 [0,15] to [30,40].) The output is shown in Figure 3; the period is clearly 1.0 s out to four significant figures.


Figure 3: The value of \phi at integer steps between t=30 and t=40 for \gamma=1.06.


Figure 4: The value of \phi at integer steps between t=30 and t=40 for \gamma=1.073


Figure 6: The value of \phi at integer steps between t=30 and t=40 for \gamma=1.077.









Next, let’s increase the forcing to \gamma=1.073. The result is shown in Figure 7. The attractor is still present (now with a value of around \phi=-2\pi), 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 \phi for integer timesteps between t=30 and t=40. The result is shown in Figure 4. The actual period appears to be 2.0 s, since the value of \phi does not repeat exactly after 1.0 s. This is called Period Doubling.


Figure 7: f=1.073


Figure 8: f=1.077









In Figure 8, I’ve displayed a plot with f=1.077, 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, \phi(0)=-\pi/2 instead of \phi(0)=0. The code I am adding is

Q2=[[i,k] for i,j,k in P2]
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, \Delta\phi(0)=1\times 10^{-4}. The code follows:

#plot(intP,0,15)+plot(intP2,0,15,linestyle=":", color=''red'')
plot(lambda x: abs(intP(x)-intP2(x)),0,15)



Figure 8: Two oscillators with f=1.073 but with different initial values, \Delta\phi(0)=-\pi/2.


Figure 9: A plot of the difference in the oscillators over time, \Delta\phi(t), for f=1.077 and a very small initial separation, \Delta \phi(0)=1\times 10^-4.









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 f=1.0829. This is displayed in Figure 11, for f=1.105. 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).


Figure 10: A log plot of two oscillators with f=1.073 but with very small differences in their initial conditions, \Delta\phi(0)=1\times 10^{-4}.


Figure 11: Finally, a log plot of two overcritical oscillators (f=1.105) and very small differences in their initial conditions, \Delta\phi(0)=1\times 10^{-4}


I have a troll!

As a practicing physicist I’ve met/interacted with a few people who can be considered Trolls – or maybe “quacks” is a more appropriate term since rather than just arguing for the sake of arguing they believe they have discovered some important feature about the world which they alone understand. I’ve been e-mailed by quacks, seen some quack seminars (always the best), and have now had a public debate with one – “public” as in the venue is the comment section on Youtube.

Of course, who cares, everyone on the internet is crazy. Well, this is my first experience so I’m recording it. I’m making some teaching videos for a partially flipped class we are teaching at Merrimack College. Last week, I posted a video about time dilation for my class to watch this week:

(click the youtube link on the bottom right to see the troll I am referring to),

Pretty quickly, I had someone named “Pentcho Valev” asking why the speed of light was constant. I was split between thinking “wow someone doesn’t understand but really wants to know more!” and thinking “uh oh”. In retrospect, I should have known what was happening as soon as a read these two lines:

“To put it simply, the frequency shifts because the speed of light shifts.”

“An alternative explanation of the frequency shift (the only salvation for relativity) involves the assumption that the motion of the observer has somehow changed the wavelength of the incoming light. […] This assumption is so obviously absurd that relativists never state it explicitly. Yet without it relativity collapses.”

Doing my due diligence as a physicist and teacher, I attempted to reason with him. But, he’s a troll and it didn’t work. Meh, no big deal.

BUT, it turns out Pentcho Valev is an entire internet quack phenomenon! There is even an entire (albeit out-of-date) website, outlining his “scholarly activities”:

So there are lots of ways one can decide informally “they have made it” – that is, not an award or publication or something. Maybe you get recognized at a conference by someone who knows your work, or the subject of something you published is a topic of debate *without* you having to inject it into the discussion manually. Well, I’m trying to teach early-career STEM majors the basics of mechanics – how to solve problems, how to use conceptual and analytic reasoning, and how to avoid common pitfalls and misunderstandings. And I’ve had a famous troll pay attention.

I’m going to count this as “I’ve made it”.

Physics in Music Instructional Videos

I’ve posted a set of short demonstration videos on the physics of music to YouTube over the past few weeks. The first two were by request from an instructor I work with at The Expert TA, and then I got inspired during one of my practice sessions.


The first is on Overtones and Beats, which I demonstrate on my double bass. It’s classic topic; I discuss how to produce the overtones series on a stringed instrument, and I also talk a little bit about beats.


The second is on Vibrato, and is just a simple illustration of what vibrato is and what it sounds like on a stringed instrument (my double bass again). I slow things down as best I can so you can really hear what is going on.


The third is on Sympathetic Resonance on a classical guitar. While I was playing on day I started noticing a variety of resonance occurring, so I thought it would be interesting to make a short video discussing the phenomena. I’ve been thinking in sympathetic resonance from a compositional standpoint for a while, and it’s a topic which might interest students as well. Resonance on a guitar is not quite as dramatic as the Tacoma Narrows bridge, but I’m thinking it might be something they can relate a little better too…


The quality of the videos is good enough but not fantastic – I’ve just done them now because I have time now. If I have more time (and technology) in the future I might update/improve them.