Numerical evaluation of the period of a limit cycle

Question

How can I calculate all the periods of the limit cycle of the Ueda-Duffing equation with forcing:

$\ddot{x} + k \dot{x} + x^3 = B \cos(t) $

for each set of parameters $(k, B)$ ?

Edit:

The equation exhibits sub-harmonic resonance for some sets of parameters (and chaotic behaviour too). E.g. for $k=0.08, B=0.2$ there are 5 coexisting attractors of period $2n\pi $ with $n=1, 2, 3$

Answer 1

I have done some exploratory python script that uses at its core a boundary value solver to find the closed loops from a systematic sweep of the relevant part of the phase space:

k=0.08; B=0.2
def odesys(u,t): return [u[1], B*np.cos(t)-k*u[1]-u[0]**3 ]
def bc(ya, yb): return yb-ya 

def norm(a): return max(abs(a))

rows, cols = 2, 3
points = [np.array([10.0,10.0])]; # have some non-relevant point 
for p in range(rows*cols): # loop over periods T=(p+1)*2*pi
    plt.subplot(rows,cols,p+1);
    n=p+1;
    t_init = np.linspace(0, 2*n*np.pi, n*10+1)
    t_sol  = np.linspace(0, 2*n*np.pi, n*500+1)
    for x in np.linspace(-0.5,0.5,5+1):
      for v in np.linspace(-0.5,0.5,4+1):
        u0 = np.array([ x, v ]);
        u_init = odeint(odesys, u0, t_init); 
        res = solve_bvp(lambda t,u:odesys(u,t), bc, t_init, u_init.T, tol=5e-3, max_nodes=n*500+1)
        print res.message,"\n",n,": ",res.sol(0)
        if res.success and min( norm(pp - res.sol(0)) for pp in points ) > 5e-3:
            res = solve_bvp(lambda t,u:odesys(u,t), bc, res.x, res.y, tol=1e-5, max_nodes=n*10000+1)
            print res.message, "\n refine to",res.sol(0) 
            if res.success and min( norm(pp - res.sol(0)) for pp in points ) > 1e-4:
                for j in range(n): points.append(res.sol(2*np.pi*j))
                u_sol = res.sol(t_sol);
                plt.plot(u_sol[0], u_sol[1], label="$u_0=$(%.6f,%.6f)"%(u_sol[0,0], u_sol[1,0])); 
    plt.grid(); plt.legend(), plt.title("T=%d$\cdot2\pi$"%n);
plt.show()

For the parameters posted in the question, $k=0.08$, $B=0.2$, screening up to period 12 (times $2\pi$) gave only periodic orbits for up to period 3.

Increasing the forcing and reducing the friction coefficient to $k=0.02$, $B=1.5$ periodic orbits were found for periods $1,3,6$ with as expected a more chaotic behavior