In [1]:
using Plots

Swing into Chaos¶

First example: a real pendulum, i.e. without small angle approx., with air resistance, and with an external driving force.

$$ \begin{align} \frac{d^2}{dt^2} \theta = -\frac{g}{L} \, \sin(\theta) - b \, \dot{\theta} + F_{\rm ext} \, \cos(\omega_d t) \end{align} $$

In [2]:
# real pendulum

g = 1.0
L = 1.0
ww = sqrt(g /L)

# external force and driving frequency
wd = 2.0/3.0

b = 0.5

function ff(t, y; Fext=0.0)
    ret = zeros(2)
    ret[1] = y[2]
    # small angle approximation (linear)
    # ret[2] = -g / L * y[1] - b * y[2] + Fext * cos(wd * t)
    # realistic case (non-linear)
    ret[2] = -g / L * sin(y[1]) - b * y[2] + Fext * cos(wd * t)
    return ret
end

Nmax = 1250
tf = 40.0
dt = tf/(Nmax-1.0)

tvec = range(0.0, tf, step = dt)
yvec = [zeros(2) for tt in tvec]

data = zeros(Nmax, 3)

# initial condition
yvec[1] = [pi-0.01, 0.0]
data[1, :] = [tvec[1], yvec[1]...]

for i1 in 1:Nmax-1

    tt = tvec[i1]
    yy = yvec[i1]

    # RK4 on the fly
    f1 = ff(tt, yy)
    f2 = ff(tt+dt/2, yy+dt*f1/2)
    f3 = ff(tt+dt/2, yy+dt*f2/2)
    f4 = ff(tt+dt, yy+dt*f3)

    yvec[i1+1] = yy + dt * (f1+2*f2+2*f3+f4)/6

    data[i1+1, 1:3]= [tt, yvec[i1+1]...]
    
end
    
p1 = plot(
    tvec, data[:, 2],
    label="num. sol.", alpha = 0.4, color=:red,
    xlabel="t", ylabel="theta(t)"
)

# ref. result
plot!(tvec, t -> data[1,2] * cos(ww * t), 
    label = "natural frequency", color="gray", line=(dash=:dash))

# p2 = plot(data[:, 2], data[:,3])
Out[2]:
No description has been provided for this image
In [3]:
function pendulum(; b=0.5, Fext=0.0, wd=2.0/3.0, tf = 70.0, yinit=[-pi+0.01, 0.0])
    
    g = 1.0
    L = 1.0
    ww = sqrt(g /L)
    
    function ff(t, y)
        ret = zeros(2)
        ret[1] = y[2]
        # small angle approximation (linear)
        # ret[2] = -g / L * y[1] - b * y[2] + Fext * cos(wd * t)
        # realistic case (non-linear)
        ret[2] = -g / L * sin(y[1]) - b * y[2] + Fext * cos(wd * t)
        return ret
    end

    dt = 0.01

    Nmax = round(Int, tf/dt)
    
    tvec = range(0.0, tf, step = dt)
    yvec = [zeros(2) for tt in tvec]
    
    data = zeros(Nmax, 3)
    
    # initial condition
    yvec[1] = yinit
    data[1, :] = [tvec[1], yvec[1]...]
    
    for i1 in 1:Nmax-1
    
        tt = tvec[i1]
        yy = yvec[i1]
    
        # RK4 on the fly
        f1 = ff(tt, yy)
        f2 = ff(tt+dt/2, yy+dt*f1/2)
        f3 = ff(tt+dt/2, yy+dt*f2/2)
        f4 = ff(tt+dt, yy+dt*f3)
    
        yvec[i1+1] = yy + dt * (f1+2*f2+2*f3+f4)/6
    
        data[i1+1, 1:3]= [tt, yvec[i1+1]...]
        
    end

    return data
    
end
Out[3]:
pendulum (generic function with 1 method)
In [4]:
data1 = pendulum(Fext=0.0)
data2 = pendulum(Fext=0.5)
data3 = pendulum(Fext=1.18)
    
p1 = plot(
    data1[:,1], data1[:,2],
    label="Fext=0.0", alpha = 0.6, color=:black,
    xlabel="t", ylabel="theta(t)"
)

plot!(data2[:,1], data2[:,2], label="Fext=0.5")
plot!(data3[:,1], data3[:,2], label="Fext=1.18")

# ref. result
plot!(data1[:,1], t -> pi*cos(t), label = "natural frequency", color="gray", line=(dash=:dash))

# p2 = plot(data3[:, 2], data3[:,3], label="phase diagram")
Out[4]:
No description has been provided for this image
In [5]:
data1 = pendulum(Fext=0.0, b=0.0)
data2 = pendulum(Fext=0.0, b=0.0, yinit=[2.25, 0.0])
data3 = pendulum(Fext=0.0)
    
p2 = plot(data1[:, 2], data1[:,3],
        label="Fext=0, b=0", alpha=0.6, color=:black, xlabel="theta", ylabel="omega")

plot!(data2[:, 2], data2[:,3], label="")
plot!(data3[:, 2], data3[:,3], label="w air resistance")
Out[5]:
No description has been provided for this image

angle wrapper¶

It is useful to restrict the angle to the range (-pi, pi].

In [6]:
# put the angle from -pi to pi
wrap = theta -> mod(theta + pi, 2*pi) - pi
Out[6]:
#17 (generic function with 1 method)
In [7]:
wrap(pi)
Out[7]:
-3.141592653589793
In [8]:
data1 = pendulum(Fext=0.0)
data2 = pendulum(Fext=0.5)
data3 = pendulum(Fext=1.15, tf=10*(3*pi))

p2 = plot(data1[:, 2], data1[:,3],
        label="Fext=0.0", alpha=0.6, color=:black, xlabel="theta", ylabel="omega")

plot!(data2[:, 2], data2[:,3], label="Fext=0.5")
plot!(wrap.(data3[:, 2]), data3[:,3], label="Fext=1.15", alpha=0.65)
scatter!((data1[1,2], data1[1,3]), label="start", color=:black)
Out[8]:
No description has been provided for this image

Tools to analyze Chaos¶

  • phase diagram: (y, ydot)-diagram
  • Poincare map Lyapunov exponent
In [9]:
function pendulum_strob(; Fext=0.5, Nsample=80)
    
    g = 1.0; L = 1.0; ww = sqrt(g/L)
    b = 0.5; wd = (2/3)*ww

    Td = 2*pi / wd

    function evolve(tf, ti, yinit)
        
        function ff(t, y)
            return [y[2], -g/L*sin(y[1]) - b*y[2] + Fext*cos(wd*t)]
        end
        function rk4(t, y, dt)
            f1 = ff(t, y)
            f2 = ff(t + dt/2, y + dt*f1/2)
            f3 = ff(t + dt/2, y + dt*f2/2)
            f4 = ff(t + dt,   y + dt*f3)
            return y + dt*(f1 + 2*f2 + 2*f3 + f4)/6
        end

        # check that it is small enough
        dt = 0.005
        
        tv = range(ti, tf, step=dt)
        yy = yinit
    
        for tt in tv
            yy = rk4(tt, yy, dt)
        end
    
        return yy
    end

    # sample stroboscopically
    tt = 0.0
    yy = [0.0, 0.0]
    
    # remove transient
    tt = 10.0*Td
    yy = evolve(tt, 0.0, yy)
    
    data = zeros(Nsample, 3)
    
    data[1, :] = [tt, yy...]
    
    for nn in 1:Nsample-1
    
        yy = evolve(tt+Td, tt, yy)
        tt += Td
    
        data[nn+1, :] = [tt, yy...]
        
    end

    return data
end

data1 = pendulum_strob(Fext=0.5)
data2 = pendulum_strob(Fext=1.0)
# we do 1500 periods
data3 = pendulum_strob(Fext=1.15, Nsample=1500)

# Create a custom colormap along the rainbow
mycolor = cgrad([:red, :orange, :yellow, :green, :blue, :indigo, :purple])
     
scatter(data1[:, 2], data1[:, 3], marker=:circle, alpha=0.8,
        xlabel="theta", ylabel="theta dot", color=:grey,
        title="Stroboscopic Poincaré Section", label="Fext=0.5",
        legend=true)
scatter!(data2[:,2], data2[:,3], alpha=0.2, label="Fext=1.0")
scatter!(wrap.(data3[:,2]), data3[:,3], alpha=0.1, label="Fext=1.15",
        marker_z = data3[:, 1], colormap = mycolor)
Out[9]:
No description has been provided for this image

further notes¶

  • Poincare Section brings order to chaos!
  • if you dial $F_{\rm ext}$ further up it becomes periodic again!
  • lines on phase diagram don't cross, but only in the volume: (t, theta, theta dot), not applicable in the projected plane.

Van der Pol equation¶

$$ \begin{align} \frac{d^2}{dt^2} y = -y - \mu \, (y^2-1) \frac{d y}{dt} \end{align} $$

This system is famous for its limited cycle.

In [10]:
function Pol(; mu=2.5, tf=50.0)

    function ff(t, y)
        return [y[2], -y[1]-mu*(y[1]^2-1)*y[2]]
    end

    dt = 0.1

    Nmax = round(Int, tf/dt) + 1

    tt = 0.0
    yvec = [pi/4.0, 0.0]
    
    data = zeros(Nmax, 3)
    
    # initial condition
    data[1, :] = [tt, yvec...]'
    
    for i1 in 2:Nmax

        yy = yvec
    
        # RK4 on the fly
        f1 = ff(tt, yy)
        f2 = ff(tt+dt/2, yy+dt*f1/2)
        f3 = ff(tt+dt/2, yy+dt*f2/2)
        f4 = ff(tt+dt, yy+dt*f3)
    
        yvec = yy + dt * (f1+2*f2+2*f3+f4)/6
    
        data[i1, :]= [tt+dt, yvec...]'

        tt += dt
        
    end

    return data
    
end

data1 = Pol(mu=0.0)
data2 = Pol(mu=0.15)
data3 = Pol(mu=2.5)

p2 = plot(data1[:, 2], data1[:,3],
        label="mu=0.0", alpha=0.6, color=:black, xlabel="theta", ylabel="omega")

plot!(data2[:, 2], data2[:, 3], label="mu=0.15", color=:green)

mycolor = cgrad([:red, :orange, :yellow, :green, :blue, :indigo, :purple])
scatter!(data3[:, 2], data3[:,3], label="mu=2.5", alpha=0.4, marker_z=data3[:,1], colormap=mycolor)
Out[10]:
No description has been provided for this image