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
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)
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
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP