Hi all,(adsbygoogle = window.adsbygoogle || []).push({});

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:

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?

- 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.0_dp
- B = -1.0_dp
- C = 0.0_dp
- D = 0.0_dp
- y = 1.0_dp
- z = 0.0_dp
- t = 0.0_dp
- do i=1,999
- call rk2(z, y, t, n)
- n = n + 1.0_dp
- write(100,*) z, y
- 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
- h = 0.1_dp
- t = n*h
- 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

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?

**Physics Forums | Science Articles, Homework Help, Discussion**

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

Have something to add?

Draft saved
Draft deleted

**Physics Forums | Science Articles, Homework Help, Discussion**