TransWikia.com

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

Computational Science Asked on October 23, 2021

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,
$$
mddot x = -gammadot x + 2ax – 4bx^3 + F_0cos(omega t)
$$

which I’ve separated into two equations by the usual approach
$$
dot x_1 = x_2
$$

and
$$
mdot x_2 = -gamma x_2 + 2ax_1-4bx_1^3+F_0cos(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=2pi/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 << "tt xt pn";
    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

Attempt at Duffing equation strange attractor

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.
all points

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

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:
Nonzero gamma

Model 2:
Negative a

Model 3:
Nonzero F0, omega

One 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.

Answered by Maxim Umansky on October 23, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP