FORTRAN 90 How to plot the Poincaré section?

[youngfreedman]

Hi all,

I'm writing a fortran 90 program to simulate the Duffing oscillator. I have it working and have produced a cos wave which I expected for D (the driving force) being set to 0. Here's the code:

1. program rungekutta
2. implicit none
3. integer, parameter :: dp = selected_real_kind(15,300)
4. integer :: i, n
5. real(kind=dp) z, y, t, A, C, D, B, omega, h
6. open(unit=100, file="rungekutta.dat",status='replace')
7. n = 0
8. !constants
9. omega = 1.0_dp
10. A = 0.0_dp
11. B = -1.0_dp
12. C = 0.0_dp
13. D = 0.0_dp
14. y = 1.0_dp
15. z = 0.0_dp
16. t = 0.0_dp
17. do i=1,999
18. call rk2(z, y, t, n)
19. n = n + 1.0_dp
20. write(100,*) z, y
21. end do
22. 
    
23. contains
24. subroutine rk2(z, y, t, n) !subroutine to implement runge-kutta algorithm
25. implicit none
26. integer, parameter :: dp = selected_real_kind(15,300)
27. integer, intent(inout) :: n
28. real(kind=dp) :: k1y, k1z, k2y, k2z, y, z, t
29. h = 0.1_dp
30. t = n*h
31. k1y = dydt(y,z,t)*h
32. k1z = dzdt(y,z,t)*h
33. k2z = dzdt(y + (0.5_dp*k1y), z + (0.5_dp*k1z), t + (0.5_dp*h))*h
34. k2y = dydt(y, z +(0.5_dp*k1z), t)*h
35. y = y + k2y
36. z = z + k2z
37. end subroutine
38. 
    
39. !2nd order ODE split into 2 for coupled Runge-Kutta, useful to define 2 functions
40. function dzdt(y,z,t)
41. real(kind=dp) :: y, z, t, dzdt
42. dzdt = -A*y**3.0_dp + B*y - C*z + D*sin(omega*t)
43. end function
44. 
    
45. function dydt(y,z,t)
46. real(kind=dp) :: z, dydt, y, t
47. dydt = z
48. end function
49. end program

So next I need to set the value for the driving force D to a small value and set it at the unstable equilibrium. What does that mean? How do I set it at the "unstable equilibrium"? Is that just y=0?

My second question is about the limit cycle. What does this look like? What should I expect? The question says to adjust the parameters until I observe one, but I don't know what one looks like. I've tried googling it but can't find anything useful.

Finally, and most importantly (as it asks me to do this often), it asks me to plot the Poincaré section. To do this, I have to take snap shots of the phase-space diagram at certain time intervals such that t*omega = 2*n*pi, n = 0,1,2,3... How do I do this, and what should it look like?

 

[The Bill]

Finally, and most importantly (as it asks me to do this often), it asks me to plot the Poincaré section. To do this, I have to take snap shots of the phase-space diagram at certain time intervals such that t*omega = 2*n*pi, n = 0,1,2,3... How do I do this, and what should it look like?

A Poincaré section of a solution to a differential equation is a cross section of the solution in the state space. The times that a solution crosses the subspace chosen for the Poincaré section are not equally spaced in general.

The Poincaré sections I've seen for the Duffing equation are of the dependent variable and its first derivative when the second derivative crosses a certain value, often zero. These crossings probably won't be exactly equally spaced in time when the oscillations are displaying interesting behaviors.