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)