1. PF Contest - Win "Conquering the Physics GRE" book! Click Here to Enter
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

FORTRAN 90 How to plot the Poincaré section?

  1. Nov 27, 2017 #1
    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. contains
    23. subroutine rk2(z, y, t, n) !subroutine to implement runge-kutta algorithm
    24. implicit none
    25. integer, parameter :: dp = selected_real_kind(15,300)
    26. integer, intent(inout) :: n
    27. real(kind=dp) :: k1y, k1z, k2y, k2z, y, z, t
    28. h = 0.1_dp
    29. t = n*h
    30. k1y = dydt(y,z,t)*h
    31. k1z = dzdt(y,z,t)*h
    32. k2z = dzdt(y + (0.5_dp*k1y), z + (0.5_dp*k1z), t + (0.5_dp*h))*h
    33. k2y = dydt(y, z +(0.5_dp*k1z), t)*h
    34. y = y + k2y
    35. z = z + k2z
    36. end subroutine

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

    42. function dydt(y,z,t)
    43. real(kind=dp) :: z, dydt, y, t
    44. dydt = z
    45. end function
    46. 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?
  2. jcsd
  3. Nov 27, 2017 #2
    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.
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted