# Homework Help: FORTRAN 90 How to plot the Poincaré section?

Tags:
1. Nov 27, 2017

### 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. 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. Nov 27, 2017

### The Bill

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.