Edward Lorenz’ Chaotic Butterfly

The butterfly effect is one of the most widely known principles of chaos theory. It has become a meme, propagating through popular culture in movies, books, TV shows and even casual conversation.

Can a butterfly flapping its wings in Florida send a hurricane to New York?

The origin of the butterfly effect is — not surprisingly — the image of a butterfly-like set of trajectories that was generated, in one of the first computer simulations of chaos theory, by Edward Lorenz.

Lorenz’ Butterfly

Excerpted from Galileo Unbound (Oxford, 2018) pg. 215

When Edward Lorenz (1917 – 2008) was a child, he memorized all perfect squares up to ten thousand.  This obvious interest in mathematics led him to a master’s degree in the subject at Harvard in 1940 under the supervision of Georg Birkhoff.  Lorenz’s master’s thesis was on an aspect of Riemannian geometry, but his foray into nonlinear dynamics was triggered by the intervention of World War II.  Only a few months before receiving his doctorate in mathematics from Harvard, the Japanese bombed Pearl Harbor. 

Lorenz left the PhD program at Harvard to join the United States Army Air Force to train as a weather forecaster in early 1942, and he took courses on forecasting and meteorology at MIT.  After receiving a second master’s degree, this time in meteorology, Lorenz was posted to Hawaii, then to Saipan and finally to Guam.  His area of expertise was in high-level winds, which were important for high-altitude bombing missions during the final months of the war in the Pacific.  After the Japanese surrender, Lorenz returned to MIT, where he continued his studies in meteorology, receiving his doctorate degree in 1948 with a thesis on the application of fluid dynamical equations to predict the motion of storms. 

One of Lorenz’ colleagues at MIT was Norbert Wiener (1894 – 1964), with whom he sometimes played chess during lunch at the faculty club.  Wiener had published his landmark book Cybernetics: Control and Communication in the Animal and Machine in 1949 which arose out of the apparently mundane problem of gunnery control during the Second World War.  As an abstract mathematician, Wiener attempted to apply his cybernetic theory to the complexities of weather, but he developed a theorem concerning nonlinear fluid dynamics which appeared to show that linear interpolation, of sufficient resolution, would suffice for weather forecasting, possibly even long-range forecasting.  Many on the meteorology faculty embraced this theorem because it fell in line with common practices of the day in which tomorrow’s weather was predicted using linear regression on measurements taken today.  However, Lorenz was skeptical, having acquired a detailed understanding of atmospheric energy cascades as larger vortices induced smaller vortices all the way down to the molecular level, dissipating as heat, and then all the way back up again as heat drove large-scale convection.  This was clearly not a system that would yield to linearization.  Therefore, Lorenz determined to solve nonlinear fluid dynamics models to test this conjecture.

Even with a computer in hand, the atmospheric equations needed to be simplified to make the calculations tractable.  Lorenz was more a scientist than an engineer, and more of a meteorologist than a forecaster.  He did not hesitate to make simplifying assumptions if they retained the correct phenomenological behavior, even if they no longer allowed for accurate weather predictions. 

He had simplified the number of atmospheric equations down to twelve.  Progress was good, and by 1961, he had completed a large initial numerical study.  He focused on nonperiodic solutions, which he suspected would deviate significantly from the predictions made by linear regression, and this hunch was vindicated by his numerical output.  One day, as he was testing his results, he decided to save time by starting the computations midway by using mid-point results from a previous run as initial conditions.  He typed in the three-digit numbers from a paper printout and went down the hall for a cup of coffee.  When he returned, he looked at the printout of the twelve variables and was disappointed to find that they were not related to the previous full-time run.  He immediately suspected a faulty vacuum tube, as often happened.  But as he looked closer at the numbers, he realized that, at first, they tracked very well with the original run, but then began to diverge more and more rapidly until they lost all connection with the first-run numbers.  His initial conditions were correct to a part in a thousand, but this small error was magnified exponentially as the solution progressed.

At this point, Lorenz recalled that he “became rather excited”.  He was looking at a complete breakdown of predictability in atmospheric science.  If radically different behavior arose from the smallest errors, then no measurements would ever be accurate enough to be useful for long-range forecasting.  At a more fundamental level, this was a break with a long-standing tradition in science and engineering that clung to the belief that small differences produced small effects.  What Lorenz had discovered, instead, was that the deterministic solution to his 12 equations was exponentially sensitive to initial conditions (known today as SIC). 

The Lorenz Equations

Over the following months, he was able to show that SIC was a result of the nonperiodic solutions.  The more Lorenz became familiar with the behavior of his equations, the more he felt that the 12-dimensional trajectories had a repeatable shape.  He tried to visualize this shape, to get a sense of its character, but it is difficult to visualize things in twelve dimensions, and progress was slow.  Then Lorenz found that when the solution was nonperiodic (the necessary condition for SIC), four of the variables settled down to zero, leaving all the dynamics to the remaining three variables. 

Lorenz narrowed the equations of atmospheric instability down to three variables: the stream function, the change in temperature and the deviation in linear temperature. The only parameter in the stream function is something known as the Prandtl Number. This is a dimensionless number which is the ratio of the kinetic viscosity of the fluid to its thermal diffusion coefficient and is a physical property of the fluid. The only parameter in the change in temperature is the Rayleigh Number which is a dimensionless parameter proportional to the difference in temperature between the top and the bottom of the fluid layer. The final parameter, in the equation for the deviation in linear temperature, is the ratio of the height of the fluid layer to the width of the convection rolls. The final simplified model is given by the flow equations

The Butterfly

Lorenz finally had a 3-variable dynamical system that displayed chaos.  Moreover, it had a three-dimensional state space that could be visualized directly.  He ran his simulations, exploring the shape of the trajectories in three-dimensional state space for a wide range of initial conditions, and the trajectories did indeed always settle down to restricted regions of state space.  They relaxed in all cases to a sort of surface that was elegantly warped, with wing-like patterns like a butterfly, as the state point of the system followed its dynamics through time.  The attractor of the Lorenz equations was strange.  Later, in 1971, David Ruelle (1935 – ), a Belgian-French mathematical physicist named this a “strange attractor”, and this name has become a standard part of the language of the theory of chaos.

The first graphical representation of the butterfly attractor is shown in Fig. 1 drawn by Lorenz for his 1963 publication.

Fig. 1 Excerpts of the title, abstract and sections of Lorenz’ 1963 paper. His three-dimensional flow equations produce trajectories that relax onto a three-dimensional “strange attractor“.

Using our modern plotting ability, the 3D character of the butterfly is shown in Fig. 2

Fig. 2 Edward Lorenz’ chaotic butterfly

A projection onto the x-y plane is shown in Fig. 3. In the full 3D state space the trajectories never overlap, but in the projection onto a 2D plane the trajectories are moving above and below each other.

Fig. 3 Projection of the butterfly onto the x-y plane centered on the origin.

The reason it is called a strange attractor is because all initial conditions relax onto the strange attractor, yet every trajectory on the strange attractor separates exponentially from neighboring trajectories, displaying the classic SIC property of chaos. So here is an elegant collection of trajectories that are certainly not just random noise, yet detailed prediction is still impossible. Deterministic chaos has significant structure, and generates beautiful patterns, without actual “randomness”.

Python Program

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 16 07:38:57 2018

@author: nolte
Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019)

Lorenz model of atmospheric turbulence
"""
import numpy as np
import matplotlib as mpl

import matplotlib.colors as colors
import matplotlib.cm as cmx

from scipy import integrate
from matplotlib import cm
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import cnames
from matplotlib import animation

plt.close('all')

jet = cm = plt.get_cmap('jet') 
values = range(10)
cNorm  = colors.Normalize(vmin=0, vmax=values[-1])
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)

def solve_lorenz(N=12, angle=0.0, max_time=8.0, sigma=10.0, beta=8./3, rho=28.0):

    fig = plt.figure()
    ax = fig.add_axes([0, 0, 1, 1], projection='3d')
    ax.axis('off')

    # prepare the axes limits
    ax.set_xlim((-25, 25))
    ax.set_ylim((-35, 35))
    ax.set_zlim((5, 55))

    def lorenz_deriv(x_y_z, t0, sigma=sigma, beta=beta, rho=rho):
        """Compute the time-derivative of a Lorenz system."""
        x, y, z = x_y_z
        return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

    # Choose random starting points, uniformly distributed from -15 to 15
    np.random.seed(1)
    x0 = -10 + 20 * np.random.random((N, 3))

    # Solve for the trajectories
    t = np.linspace(0, max_time, int(500*max_time))
    x_t = np.asarray([integrate.odeint(lorenz_deriv, x0i, t)
                      for x0i in x0])

    # choose a different color for each trajectory
    # colors = plt.cm.viridis(np.linspace(0, 1, N))
    # colors = plt.cm.rainbow(np.linspace(0, 1, N))
    # colors = plt.cm.spectral(np.linspace(0, 1, N))
    colors = plt.cm.prism(np.linspace(0, 1, N))

    for i in range(N):
        x, y, z = x_t[i,:,:].T
        lines = ax.plot(x, y, z, '-', c=colors[i])
        plt.setp(lines, linewidth=1)

    ax.view_init(30, angle)
    plt.show()

    return t, x_t


t, x_t = solve_lorenz(angle=0, N=12)

plt.figure(2)
lines = plt.plot(t,x_t[1,:,0],t,x_t[1,:,1],t,x_t[1,:,2])
plt.setp(lines, linewidth=1)
lines = plt.plot(t,x_t[2,:,0],t,x_t[2,:,1],t,x_t[2,:,2])
plt.setp(lines, linewidth=1)
lines = plt.plot(t,x_t[10,:,0],t,x_t[10,:,1],t,x_t[10,:,2])
plt.setp(lines, linewidth=1)

To explore the parameter space of the Lorenz attractor, the key parameters to change are sigma (the Prandtl number), r (the Rayleigh number) and b on line 31 of the Python code.

References

[1] E. N. Lorenz, The essence of chaos (The Jessie and John Danz lectures; Jessie and John Danz lectures.). Seattle :: University of Washington Press (in English), 1993.

[2] E. N. Lorenz, “Deterministic Nonperiodic Flow,” Journal of the Atmospheric Sciences, vol. 20, no. 2, pp. 130-141, 1963 (1963)

Georg Duffing’s Equation

Although coal and steam launched the industrial revolution, gasoline and controlled explosions have sustained it for over a century.  After early precursors, the internal combustion engine that we recognize today came to life in 1876 from the German engineers Otto and Daimler with later variations by Benz and Diesel.  In the early 20th century, the gasoline engine was replacing coal and oil in virtually all mobile conveyances and had become a major industry attracting the top mechanical engineering talent.  One of those talents was the German engineer Georg Duffing (1861 – 1944) whose unlikely side interest in the quantum mechanics revolution brought him to Berlin to hear lectures by Max Planck, where he launched his own revolution in nonlinear oscillators.

The publication of this highly academic book by a nonacademic would establish Duffing as the originator of one of the most iconic oscillators in modern dynamics.

An Academic Non-Academic

Georg Duffing was born in 1861 in the German town of Waldshut on the border with Switzerland north of Zurich.  Within a year the family moved to Mannheim near Heidelberg where Georg received a good education in mathematics as well as music.  His mathematical interests attracted him to engineering, and he built a reputation that led to an invitation to work at Westinghouse in the United States in 1910.  When he returned to Germany he set himself up as a consultant and inventor with the freedom to move where he wished.  In early 1913 he wished to move to Berlin where Max Planck was lecturing on the new quantum mechanics at the University.  He was always searching for new knowledge, and sitting in on Planck’s lectures must have made him feel like he was witnessing the beginnings of a new era.            

At that time Duffing was interested in problems related to brakes, gears and engines.  In particular, he had become fascinated by vibrations that often were the limiting factors in engine performance.  He stripped the problem of engine vibration down to its simplest form, and he began a careful and systematic study of nonlinear oscillations.  While in Berlin, he had became acquainted with Prof. Meyer at the University who had a mechanical engineering laboratory.  Meyer let Duffing perform his experiments in the lab on the weekends, sometime accompanied by his eldest daughter.  By 1917 he had compiled a systematic investigation of various nonlinear effects in oscillators and had written a manuscript that collected all of this theoretical and experimental work.  He extended this into a small book that he published with Vieweg & Sohn in 1918 to be purchased for a price of 5 Deutsch Marks [1].   The publication of this highly academic book by a nonacademic would establish Duffing as the originator of one of the most iconic oscillators in modern dynamics.

Fig. 1 Cover of Duffing’s 1918 publication on nonlinear oscillators.

Duffing’s Nonlinear Oscillator

The mathematical and technical focus of Duffing’s book was low-order nonlinear corrections to the linear harmonic oscillator.  In one case, he considered a spring that either became stiffer or softer as it stretched.  This happens when a cubic term is added to the usual linear Hooke’s law.  In another case, he considered a spring that was stiffer in one direction than another, making the stiffness asymmetric.  This happens when a quadratic term is added.  These terms are shown in Fig. 2 from Duffing’s book.  The top equation is a free oscillation, and the bottom equation has a harmonic forcing function.  These were the central equations that Duffing explored, plus the addition of damping that he considered in a later chapter as shown in Fig. 3. The book lays out systematically, chapter by chapter, approximate and series solutions to the nonlinear equations, and in special cases described analytically exact solutions (such as for the nonlinear pendulum).

Fig. 2 Duffing’s equations without damping for free oscillation and driven oscillation with quadratic (producing an asymmetric potential) and cubic (producing stiffening or softening) corrections to the spring force.
Fig. 3 Inclusion of damping in the case with cubic corrections to the spring force.

Duffing was a practical engineer as well as a mathematical one, and he built experimental systems to test his solutions.  An engineering drawing of his experimental test apparatus is shown in Fig. 4. The small test pendulum is at S in the figure. The large pendulum at B is the drive pendulum, chosen to be much heavier than the test pendulum so that it can deliver a steady harmonic force through spring F1 to the test system. The cubic nonlinearity of the test system was controlled through the choice of the length of the test pendulum, and the quadratic nonlinearity (the asymmetry) was controlled by allowing the equilibrium angle to be shifted from vertical. The relative strength of the quadratic and cubic terms was adjusted by changing the position of the mass at G. Duffing derived expressions for all the coefficients of the equations in Fig. 1 in terms of experimentally-controlled variables. Using this apparatus, Duffing verified to good accuracy his solutions for various special cases.

Fig. 4 Duffing’s experimental system he used to explore and verify his equations and solutions.

           Duffing’s book is a masterpiece of careful systematic investigation, beginning in general terms, and then breaking the problem down into its special cases, finding solutions for each one with accurate experimental verifications. These attributes established the importance of this little booklet in the history of science and technology, but because it was written in German, most of the early citations were by German scientists.  The first use of Duffing’s name associated to the nonlinear oscillator problem occurred in 1928 [2], as was the first reference to him in a work in English in a book by Timoshenko [3].  The first use of the phrase “Duffing Equation” specifically to describe an oscillator with a linear and cubic restoring force was in 1942 in a series of lectures presented at Brown University [4], and this nomenclature had become established by the end of that decade [5].  Although Duffing had spent considerable attention in his book to the quadratic term for an asymmetric oscillator, the term “Duffing Equation” now refers to the stiffening and softening problem rather than to the asymmetric problem.

Fig. 5 The Duffing equation is generally expressed as a harmonic oscillator (first three terms plus the harmonic drive) modified by a cubic nonlinearity and driven harmonically.

Duffing Rediscovered

Nonlinear oscillations remained mainly in the realm of engineering for nearly half a century, until a broad spectrum of physical scientists began to discover deep secrets hiding behind the simple equations.  In 1963 Edward Lorenz (1917 – 2008) of MIT published a paper that showed how simple nonlinearities in three equations describing the atmosphere could produce a deterministic behavior that appeared to be completely chaotic.  News of this paper spread as researchers in many seemingly unrelated fields began to see similar signatures in chemical reactions, turbulence, electric circuits and mechanical oscillators.  By 1972 when Lorenz was invited to give a talk on the “Butterfly Effect” the science of chaos was emerging as new frontier in physics, and in 1975 it was given its name “chaos theory” by James Yorke (1941 – ).  By 1976 it had become one of the hottest new areas of science. 

        Through the period of the emergence of chaos theory, the Duffing oscillator was known to be one of the archetypical nonlinear oscillators.  A particularly attractive aspect of the general Duffing equations is the possibility of studying a “double-well” potential.  This happens when the “alpha” in the equation in Fig. 5 is negative and the “beta” is positive.  The double-well potential has a long history in physics, both classical and modern, because it represents a “two-state” system that exhibits bistability, bifurcations, and hysteresis.  For a fixed “beta” the potential energy as a function of “alpha” is shown in Fig. 6.  The bifurcation cascades of the double-well Duffing equation was investigated by Phillip Holmes (1945 – ) in 1976 [6], and the properties of the strange attractor were demonstrated in 1978 [7] by Yoshisuke Ueda (1936 – ).  Holmes, and others, continued to do detailed work on the chaotic properties of the Duffing oscillator, helping to make it one of the most iconic systems of chaos theory.

Fig. 6 Potential energy of the Duffing Oscillator. The position variable is x, and changing alpha is along the other axis. For positive beta and alpha the potential is a quartic. For positive beta and negative alpha the potential is a double well.

Python Code for the Duffing Oscillator: Duffing.py

This Python code uses the simple ODE solver on the driven-damped Duffing double-well oscillator to display the configuration-space trajectories and the Poincaré map of the strange attractor.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Duffing.py
Created on Wed May 21 06:03:32 2018
@author: nolte
D. D. Nolte, Introduction to Modern Dynamics: Chaos, Networks, Space and Time, 2nd ed. (Oxford,2019)
"""
import numpy as np
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from scipy import integrate
from matplotlib import pyplot as plt
from matplotlib import cm
import time
import os

plt.close('all')

# model_case 1 = Pendulum
# model_case 2 = Double Well
print(' ')
print('Duffing.py')

alpha = -1       # -1
beta = 1         # 1
delta = 0.3       # 0.3
gam = 0.15    # 0.15
w = 1
def flow_deriv(x_y_z,tspan):
    x, y, z = x_y_z
    a = y
    b = delta*np.cos(w*tspan) - alpha*x - beta*x**3 - gam*y
    c = w
    return[a,b,c]
                
T = 2*np.pi/w

px1 = np.random.rand(1)
xp1 = np.random.rand(1)
w1 = 0

x_y_z = [xp1, px1, w1]

# Settle-down Solve for the trajectories
t = np.linspace(0, 2000, 40000)
x_t = integrate.odeint(flow_deriv, x_y_z, t)
x0 = x_t[39999,0:3]

tspan = np.linspace(1,20000,400000)
x_t = integrate.odeint(flow_deriv, x0, tspan)
siztmp = np.shape(x_t)
siz = siztmp[0]

y1 = x_t[:,0]
y2 = x_t[:,1]
y3 = x_t[:,2]
    
plt.figure(2)
lines = plt.plot(y1[1:2000],y2[1:2000],'ko',ms=1)
plt.setp(lines, linewidth=0.5)
plt.show()

for cloop in range(0,3):

#phase = np.random.rand(1)*np.pi;
    phase = np.pi*cloop/3

    repnum = 5000
    px = np.zeros(shape=(2*repnum,))
    xvar = np.zeros(shape=(2*repnum,))
    cnt = -1
    testwt = np.mod(tspan-phase,T)-0.5*T;
    last = testwt[1]
    for loop in range(2,siz):
        if (last < 0)and(testwt[loop] > 0):
            cnt = cnt+1
            del1 = -testwt[loop-1]/(testwt[loop] - testwt[loop-1])
            px[cnt] = (y2[loop]-y2[loop-1])*del1 + y2[loop-1]
            xvar[cnt] = (y1[loop]-y1[loop-1])*del1 + y1[loop-1]
            last = testwt[loop]
        else:
            last = testwt[loop]
 
    plt.figure(3)
    if cloop == 0:
        lines = plt.plot(xvar,px,'bo',ms=1)
    elif cloop == 1:
        lines = plt.plot(xvar,px,'go',ms=1)
    else:
        lines = plt.plot(xvar,px,'ro',ms=1)
        
    plt.show()

plt.savefig('Duffing')
Fig. 7 Strange attractor of the double-well Duffing equation for three selected phases.

[1] G. Duffing, Erzwungene Schwingungen bei veranderlicher Eigenfrequenz und ihre technische Bedeutung, Vieweg & Sohn, Braunschweig, 1918.

[2] Lachmann, K. “Duffing’s vibration problem.” Mathematische Annalen 99: 479-492. (1928)

[3] S. Timoshenko, Vibration Problems in Engineering, D. Van Nostrand Company, Inc.,New York, 1928.

[4] K.O. Friedrichs, P. Le Corbeiller, N. Levinson, J.J. Stoker, Lectures on Non-Linear Mechanics delivered at Brown University, New York, 1942.

[5] Kovacic, I. and M. J. Brennan, Eds. The Duffing Equation: Nonlinear Oscillators and their Behavior. Chichester, United Kingdom, Wiley. (2011)

[6] Holmes, P. J. and D. A. Rand. “Bifurcations of Duffings Equation – Application of Catastrophe Theory.” Journal of Sound and Vibration 44(2): 237-253. (1976)

[7] Ueda, Y. “Randomly Transitional Phenomena in the System Governed by Duffings Equation.” Journal of Statistical Physics 20(2): 181-196. (1979)