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

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 = x;
dxdt = (1/m)*(-gam*x+2*a*x-4*b*x*x*x+F0*cos(omega*t));
}

void write_duffing(const state_type &x, const double t, ofstream& outfile) {
outfile << t << "\t" << x << "\t" << x << 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?