# Forced damped pendulum and route to chaos
Using Newton's second law we can derived following equation of motion.
\begin{equation}
\ddot{\theta} + \gamma\dot{\theta} + \dfrac{g}{l}\sin{\theta} = \dfrac{F}{ml}\cos{\omega t}
\tag{1}
\end{equation}
Or can be rewrite as normalize unit
\begin{equation}
\ddot{\theta} + \dfrac{1}{q}\dot{\theta} + \sin{\theta} = f\cos{\omega_D t}
\tag{2}
\end{equation}
Thus, we can seperate this second order equation into three first order. Using these system of ODE we can applying RK4 method to solve it.
\begin{equation}
\begin{split}
&\dot{\theta} = \omega\\
&\dot{\omega} = -\dfrac{\omega}{q}-\sin{\theta}+f\cos{\phi}\\
&\dot{\phi} = \omega_D\\
\end{split}
\tag{3}
\end{equation}
Now, the parameter space has three component that we can mess around, since the nolinearity of this system by slightly varing the value of parameter we can get a totally different and sometimes chaotic result.
For the following result I use below a set of paramter and varing f value.
\begin{equation*}
q=2, \quad\quad\omega_D=\dfrac{2}{3}
\end{equation*}

## Result

### f<1
You can see if f is under 1, the natural frequency of this system we get quite regular behavior, nothing chaotic at all. The Poincare
section are mark as the blue cross ($t_n=5,100$) and phase space are the red line ($5<t<100$), again verified the observation of chaotic system.



### f>1
However, when f is close to 1 or greater than 1, we get very different outcome, some may settle to a periodic motion others become quite messy. Moreover, the Poincare section of this chaotic state look very periodic as if some pattern are undereath the mysterious behavior.




### Bifurcation diagram

If we then vary the value of f from 0 to 2 and plot the Poincare's section of each f we can get above beautiful graph (bifurcation diagram). What this graph tell us is that the periodic behavior last till the frequency of external force close to natural frequency in this case is 1. As f increase, one branch become two branch, two branch become four branch and the route leads to chaotic. The gray part of the graph indicate chaotic behavior, we also see that there are certain range of value of f in which the system has a stable preiodic behavior.

## Code

program sysrk4
!========================================================
!   purpose: solve system of ode
!
!   methods: rk4
!
!   input: f value
!
!   output: in the order of theta, omega, phi
!
!       date        programer       description of change
!       ====        =========       =====================
!     11/17/20       morrish        original code
!========================================================
    implicit none
    double precision, external :: f1, f2, f3
    double precision :: w(4, 3), h, f, pi=4*atan(1.d0)
    double precision, dimension(:), allocatable :: theta, omega, phi
    integer :: i, n
    character(len=10) :: filename

    ! read input
    write(*, *) 'input f value of the system'
    read(*, *) f
    write(filename, '(1I3.3)') nint(f*100.d0)

    ! time step
    h = 1.d-4
    n = nint(1000/h)
    allocate(theta(n))
    allocate(omega(n))
    allocate(phi(n))
    ! inital value
    theta(1) = .1d0
    omega(1) = 0.d0
    phi(1) = 0.d0

    do i = 1, n-1
        w(1, 1) = h * f1(omega(i))
        w(1, 2) = h * f2(omega(i), theta(i), phi(i), f)
        w(1, 3) = h * f3()

        w(2, 1) = h * f1(omega(i)+0.5d0*w(1, 1))
        w(2, 2) = h * f2(omega(i)+0.5d0*w(1, 1), &
                         theta(i)+0.5d0*w(1, 2), &
                         phi(i)+0.5d0*w(1, 3), f)
        w(2, 3) = h * f3()

        w(3, 1) = h * f1(omega(i)+0.5d0*w(2, 1))
        w(3, 2) = h * f2(omega(i)+0.5d0*w(2, 1), &
                         theta(i)+0.5d0*w(2, 2), &
                         phi(i)+0.5d0*w(2, 3), f)
        w(3, 3) = h * f3()

        w(4, 1) = h * f1(omega(i)+w(3, 1))
        w(4, 2) = h * f2(omega(i)+w(3, 1), theta(i)+w(3, 2), &
                         phi(i)+w(3, 3), f)
        w(4, 3) = h * f3()

        theta(i+1) = theta(i) + (w(1, 1) + 2.d0 * w(2, 1) + &
                     2.d0 * w(3, 1) + w(4, 1)) / 6.d0
        if (theta(i+1) < -pi) then
            theta(i+1) = pi + (theta(i+1) + pi)
        elseif (theta(i+1) > pi) then
            theta(i+1) = - pi + (theta(i+1) - pi)
        end if
        omega(i+1) = omega(i) + (w(1, 2) + 2.d0 * w(2, 2) + &
                     2.d0 * w(3, 2) + w(4, 2)) / 6.d0
        phi(i+1) = phi(i) + (w(1, 3) + 2.d0 * w(2, 3) + &
                   2.d0 * w(3, 3) + w(4, 3)) / 6.d0

    end do

    open(66, file='data/'//trim(filename)//'.dat', status='unknown')
    do i = 1, n
        if ((phi(i) >= 10.d0 / 3.d0).and.(phi(i) <= 200.d0 / 3.d0)) then
            write(66, *) theta(i), omega(i), phi(i)
        end if
    end do
end program sysrk4

function f1(omega)
    implicit none
    double precision :: f1, omega
    f1 = omega
end function

function f2(omega, theta, phi, f)
    implicit none
    double precision :: f2, omega, theta, phi, f
    double precision :: q=2.d0
    f2 = - omega / q - dsin(theta) + f * dcos(phi)
end function

function f3()
    implicit none
    double precision :: f3
    f3 = 2.d0 / 3.d0
end function