How do I plot a Poincare plot for a non-autonomous ODE's?

I have system of damped driven double pendulum which is given by a set of 5 ODE's. I want to plot the Poincare map of the solution at different parameters of the system. I have read some references where the system is Hamiltonian and for a fixed Energy Poincare map is obtained. As for the technique l Iearned that to plot the Poincare map one needs to collect the points ( and ) when the lower pendulum crosses the vertically down position with positive velocity. Now the problem is how to write a code for this? I can integrate the dynamical set of equations to see the time evolution of data and may a get a periodic solution (ofcourse have to discared the transient) . Now to collect the data for ( and ) when is zero and is positive may not be possible as in my time evolution data there may not be points with is zero.

It would be very helpful if one can suggest what would be best practice to plot the Poincare map for this system.

Answer 1

[Bjorn Gustavsson]

This sound like a perfect example of when one should use the events capabilities of the ODE-solvers. If you look at the code of the ballode demo you'll see how to set this up. The difference in your case is that you simply track and collect the events. Something like this would do it for an events-function:

function [value,isterminal,direction] = events(t,y)
% Locate the time when theta2 passes through zero and dtheta2 is positive direction
% and stop integration.
value = y(2);     % detect zero-crossing = 0, if you have theta2 in the second component of y
isterminal = 0;   % stop the integration
direction = 1;   % positive direction, or -1 if I got this one wrong to get dtheta2 > 0 

Then set the options-struct to something like this:

refine = 4;
options = odeset('Events',@events,'Refine',refine);

Then it should be possible to run the ode-function something like this:

tstart = 0;
tfinal = 300;        % whatever time-span you're interested in
y0 = [0; 0,10,20,0]; % whatever your initial conditions are
[t,y,te,ye,ie] = ode23(@f,[tstart tfinal],y0,options); % or one of the other odeNN-functions

This would give you the solution for all parameters at all the theta_2 crossings in variable ye at time te.

HTH

[Bjorn Gustavsson]

Did this help solve your problem?

[Rebeka]

Thanks a lot. Yes,this solves my problem.