# 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?

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) • Can you show the phase diagram in$x, \dot{x}$space including all points on the trajectory? The cases to look at would those with only one term on the right-hand side nonzero, then it will be easier to see if this makes sense. Jul 19, 2020 at 3:06 • I've added a plot for all zero except$\gamma$and$F_0$as well as the original plot with all the points. – tmph Jul 19, 2020 at 5:41 ## 1 Answer 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. • Thanks for the tip. I've edited my post to include each of these models (you may want to check your first case, as I am pretty sure your solution is wrong). I don't think the problem is with solving the differential equation itself; I think the problem is how I am getting the strange attractor terms (I don't know what they really are, to be honest and I think I am misunderstanding the textbook). – tmph Jul 20, 2020 at 0:17 • Ok, I corrected the first case. Are your numerical solutions plotted with analytic solutions? And the match is pretty much exact? Then the only one left is the cubic term. Jul 20, 2020 at 0:35 • I have added a solution for the cubic term, so it should be also checked, just in case (although in the code it does not look there is anything special about the cubic term). Jul 20, 2020 at 1:11 • Yes, all the other terms were plotted with analytic solutions (sorry if that's not clear) and seem to agree. The cubic terms gave me reasonable solutions when I played with the parameters (e.g.$b=0.5$), and it gave me qualitatively what I would expect based on e.g. Wolfram Alpha's output. However, using your parameters I get a solution that seems to blow up very quickly and give NaNs. – tmph Jul 20, 2020 at 1:32 • Actually, it seems to me that I am only getting solutions if I assume$b>0$. I don't know why that would be the case; I get it just blows up really quickly. Nevertheless, I am only interested in$b>0\$ anyway.
– tmph
Jul 20, 2020 at 1:37