How to plot the Poincaré Section? (Duffing Oscillator)

Question

I've written a program that successfully shows a simple limit cycle for the Duffing equation. However, I now need to plot the Poincaré section for this case.

I need to do this by taking snapshots of the Phase-Space diagram at regular time intervals, such that t*omega = 2*pi*n. As I have omega set to 1 for this case, this is just when t = 2*pi*n. I've attempted this, but am not getting the Poincaré section I expect.

Here's my code:

program rungekutta
implicit none
integer, parameter :: dp = selected_real_kind(15,300)
integer :: i, n
real(kind=dp) z, y, t, A, C, D, B, omega, h
open(unit=100, file="rungekutta.dat",status='replace')
n = 0
!constants
omega = 1.0_dp
A = 0.25_dp
B = 1.0_dp
C = 0.1_dp
D = 1.0_dp
y = 0.0_dp
z = 0.0_dp
t = 0.0_dp
do i=1,1000
call rk2(z, y, t, n)
n = n + 1.0_dp
write(100,*) y, z
end do

contains
subroutine rk2(z, y, t, n) !subroutine to implement runge-kutta algorithm
implicit none
integer, parameter :: dp = selected_real_kind(15,300)
integer, intent(inout) :: n
real(kind=dp) :: k1y, k1z, k2y, k2z, y, z, t, pi
pi = 4.0*ATAN(1.0)
h = 0.1_dp
t = n*2*pi
k1y = dydt(y,z,t)*h
k1z = dzdt(y,z,t)*h
k2z = dzdt(y + (0.5_dp*k1y), z + (0.5_dp*k1z), t + (0.5_dp*h))*h
k2y = dydt(y, z +(0.5_dp*k1z), t)*h
y = y + k2y
z = z + k2z
end subroutine

!2nd order ODE split into 2 for coupled Runge-Kutta, useful to define 2 
functions
function dzdt(y,z,t)
real(kind=dp) :: y, z, t, dzdt
dzdt = -A*y**3.0_dp + B*y - C*z + D*sin(omega*t)
end function

function dydt(y,z,t)
real(kind=dp) ::  z, dydt, y, t
dydt = z
end function
end program

I have also attached an image of what my Poincaré section looks like:

.

This is y on the x axis vs dydt.

And an image of what I'd expect:

Answer 1

In your rk2 routine you perform one step of step length 0.1. Thus the plot is the full trajectory of the solution at that resolution. However the intend seems to be to integrate over a full period length. This would require a loop in that routine.

In other words, what you want is the plot of (y(n*T), z(n*T)) where T is one of the periods of the system, per your code T=2*p. What you actually compute is (y(n*h), z(n*h)) where h=0.1 is the step size of a single step of RK2.

Also the arguments of k2y need to be corrected as per the comment of user5713492

With a corrected integrator you should get something like the following picture:

where the red squares are the points at t=n*2*pi. The indicated step size by the dots on the solution curve is the same h=0.1, the integration is over t=0..300.

def RK2(f,u,times,subdiv = 1):
     uout = np.zeros((len(times),)+u.shape)
     uout[0] = u;
     for k in range(len(times)-1):
         t = times[k]
         h = (times[k+1]-times[k])/subdiv
         for j in range(subdiv):
            k1 = f(u,t)*h
            k2 = f(u+0.5*k1, t+0.5*h)*h
            u, t = u+k2, t+h
         uout[k+1]=u
     return uout

def plotphase(A,B,C,D):
     def derivs(u,t): y,z = u; return np.array([ z, -A*y**3 + B*y - C*z + D*np.sin(t) ])
     N=60
     u0 = np.array([0.0, 0.0])
     t  = np.arange(0,300,2*np.pi/N); 
     u  = RK2(derivs, u0, t, subdiv = 10)
     plt.plot(u[:-2*N,0],u[:-2*N,1],'.--y', u[-2*N:,0],u[-2*N:,1], '.-b', lw=0.5, ms=2);
     plt.plot(u[::N,0],u[::N,1],'rs', ms=4); plt.grid(); plt.show()

plotphase(0.25, 1.0, 0.1, 1.0)