Calculating the Strange Attractor of the Duffing Oscillator in C++

Question

I am simultaneously trying to learn computational physics methods, chaos, and C++. I think this is the right site for the question, and I apologise if not.

I started working through Thijssen's Computational Physics textbook, and the first question (exercise 1.1b) is to solve the Duffing equation, $$ m\ddot x = -\gamma\dot x + 2ax - 4bx^3 + F_0\cos(\omega t) $$ which I've separated into two equations by the usual approach $$ \dot x_1 = x_2 $$ and $$ m\dot x_2 = -\gamma x_2 + 2ax_1-4bx_1^3+F_0\cos(\omega t). $$

I am trying to get the plot for the strange attracter (which from google looks like it might also be called the Poincaré map?), where as I understand it you just output $x$ and $\dot x$ at every $T=2\pi/\omega$, and plot $x$ vs $\dot x$. Currently my approach is to solve the equation with boost's odeint, and output every $T$ to a file "duffing.txt".

Here is my code (apologies for the (ab)use of lambda functions)

#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

#include <iostream> 
#include <fstream>

typedef boost::array<double,2> state_type;

void duffing(const state_type &x, state_type &dxdt, double t, double F0, double omega,
            double gam, double m, double a, double b) { 
    dxdt[0] = x[1];
    dxdt[1] = (1/m)*(-gam*x[1]+2*a*x[0]-4*b*x[0]*x[0]*x[0]+F0*cos(omega*t));
}

void write_duffing(const state_type &x, const double t, ofstream& outfile) { 
    outfile << t << "\t" << x[0] << "\t" << x[1] << endl;
}

int main(int argc, char **argv) { 
    state_type x = {0.5, 0.}; // initial conditions {x0,dxdt0}

    // parameters
    const double m = 1.;
    const double a = 0.25;
    const double b = 0.5;
    const double F0 = 2.0;
    const double omega = 2.4;
    const double gam = 0.1; 
    const double T = 2*M_PI/omega;

    string filename = "duffing.txt";
    double t0 = 0.0;
    double t1 = 10000*T;
    double dt = T/200.;
    auto f = [F0, omega, gam, m, a, b](const state_type &x, state_type &dxdt, double t) {
        duffing(x, dxdt, t, F0, omega, gam, m, a, b); };
    ofstream outfile;
    outfile.open(filename);
    outfile << "t\t x\t p\n";
    double last_t = 0;
    auto obs = [&outfile, T, &last_t](state_type &x, const double t){
        if (abs(t-last_t)>=T){
            write_duffing(x,t,outfile);
            last_t = t; 
        }
        };
    auto rkd = runge_kutta_dopri5<state_type>{};
    auto stepper = make_dense_output(1.0e-9, 1.0e-9, rkd);
    integrate_const(stepper,f, x, t0, t1, dt, obs);

    outfile.close();
    
    return 0;
}

Plotting with Gnuplot, however, this is the output of plot "duffing.txt" using 2:3 with linespoints

which is basically just an oval and doesn't seem chaotic at all. I have played with the parameters without much luck (the ones in the code are from the textbook, which includes a clearly chaotic plot, which I'm not sure is okay to rehost here).

It doesn't seem like the mistake is the integration routine since if I replace my equation with the Lorenz equations I get back the solution shown in the odeint examples. Am I going about printing it at the wrong time, or some other conceptual mistake?

Edit: as requested in the comments, here is the plot with all the points.

Here is also the plot for all terms on the RHS=0 except omega=2.4 and F0=2.0.

Unless I need to review my undergrad calculus, I think this is what is expected. Why am I not seeing a strange attractor for the more complicated case?

Edit 2: Here are the results for the "reduced models" as suggested by Maxim Umansky. The results seem to match! It doesn't seem to be a problem with how I set up my integrator, just something about how I am extracting the strange attractor... (each case has $x=0.5$ and $\dot x=0.5)

Model 1:

Model 2:

Model 3:

Answer 1

For debugging the code, there is a set of analytic solutions here for several reduced models corresponding to subsets of terms on the right-hand side. These analytic solutions have to be reproduced by the code. Verification testing of this kind is a standard practice for debugging simulation models.

Reduced model 1:

$ m \ddot{x} = - \gamma \dot{x} $

Solution: $ x = x_0 + v_0 \tau [1 - \exp(-t/\tau)] $

where $\tau = m/\gamma$

Reduced model 2:

$ m \ddot{x} = 2 a {x} $

Assume $a<0$, then

Solution: $ x = x_{0} \cos(\Omega t) + (v_{0}/\Omega) \sin(\Omega t), $

where $\Omega= (-2 a /m)^{1/2}$

Reduced model 3:

$ m \ddot{x} = F_0 \cos(\omega t) $

Solution: $ x = x_0 + v_0 t + \frac{F_0}{\omega^2} (1 - \cos(\omega t)), $

Reduced model 4:

$ \ddot{x} = - \beta x^3, $

where $\beta = - 4 b/m$.

This is a nonlinear problem, so finding a general solution is difficult; but we can easily find a particular solution.

Solution: $ x = \alpha / t, $

where $\alpha^2 = -2 m/\beta$, and initial conditions at $t=1$ are $x_{t=1}=\alpha$, $v_{t=1} = -\alpha$. We are interested in real-valued $\alpha$ so $\beta$ is negative (so $b$ is positive), and $\alpha$ can take one of the real-valued square root values. For example, for $m=1$, $\beta=-2$ (i.e., $b=1/2$), $\alpha=1$, and the solution is $x=\alpha/t$, for initial conditions at t=1: $x_1=1$, $v_1=-1$.

Most likely the bugs in the code will be found in the process of verifying these analytic solutions; or at least the search for bugs will be greatly simplified after these solutions are successfully reproduced.