Duffing attractor


Author: GC

Date: July 19, 2016

Next SI model


We consider the Duffing equation with a non autonomous term: \[ \ddot{x} + k \dot{x} + x^3 = C \cos(t). \] The following program computes its orbit passing through $(3, 4)$ at $t=0$, with $k=0.05$ and $C=7.5$. One possibility to create an animation is to generate a set of figures, and to use the following bash command:

convert -delay 40 -loop 0 *.png duffing.gif



Duffing attractor




Source code

#!/usr/bin/env python3

"""
Duffing attractor
"""


# scientific libraries
from matplotlib import pyplot as plt
import numpy as np
from scipy.integrate import odeint

# custom arrow
def myarrow(plt,begin,end,dim,col):
    return plt.arrow(begin[0],begin[1],end[0],end[1],head_width=dim[0],head_length=dim[1],color=col)

# parameters of the system
k = 0.05
C = 7.5

# initial condition
X0 = [3, 4]

# equations 
def duffing(X,t):
    x, y = X    
    dx = y
    dy = -k*y - x**3 + C*np.cos(t)
    return [dx, dy]

# figure
fig, ax = plt.subplots()
for T in range(2, 100):
    time = np.arange(0,T,0.01)
    orbit = odeint(duffing, X0, time)
    x, y = orbit.T
    ax.plot(x, y, 'b')
    myarrow(plt,[x[-2],y[-2]], [(x[-1]-x[-2])/20,(y[-1]-y[-2])/20], [0.15, 0.15], 'b')
    ax.set_ylim((-8, 8))
    ax.set_xlim((-4, 4))
    ax.set_xlabel(r"$x$",fontsize=24)
    ax.set_ylabel(r"$y$",fontsize=24)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    if T<10:fig.savefig('~/img0'+str(T)+'.png')
    else:fig.savefig('~/img'+str(T)+'.png')
    ax.clear()


Next SI model