Explanation Request Poincaré Map

Question

How do I compute this Poincare Map example (Duffing oscillator): https://en.wikipedia.org/wiki/Duffing_equation#/media/File:Forced_Duffing_equation_Poincar%C3%A9_section.png

This question purpose is to check wether I understand the concept of the Poincaré Map for non-autonomous sytems, because I'll need this for a college project. I start to write the equation as a system of ODE:

\begin{cases} x'=y \\ y'=\gamma cos(\omega t)-\delta y-\alpha x - \beta x^3 \end{cases}

I solve this numerically with ode45 from MATLAB, and I plot every point $(x(\frac{2\pi}{\omega}k), y(\frac{2\pi}{\omega}k))$. I'm using the same parameters used in the example $(\alpha,\beta,\delta,\gamma,\omega)=(1,5,0.02,8,0.5)$ but I'm not obtaining anything similar to the image. I'd like to know what is wrong in here.

Answer 1

You need to do a scatter plot, or a general plot without lines, only with markers. And you need to discard an initial segment where the solution converges towards the attractor.

It appears as if the solution converges somewhat to a limit cycle, so use different initial points to fill the attractor faster.

% set up
a=1; b=5; c=8; d=0.02; w=0.5;
f=@(t,x) [x(2); c*cos(w*t)-d*x(2)-a*x(1)-b*x(1).^3];
opts = odeset('AbsTol',1e-4,'RelTol',1e-5);
% compose the plot, 
% iterate over solutions from random initial points
hold on;
for batch=1:50
  batch
  % initial segment
  t = [0:1:30]*2*pi/w;
  [T,X] = ode45(f,t,4*rand(2,1)-2,opts);
  % get plot points
  t = [0:1:150]*2*pi/w;
  [T,X] = ode45(f,t,X(end,:),opts);
  plot(X(:,1),X(:,2),'.');
end%for
hold off;