title: Solution of double pendulum by Runge-Kutta fourth order FORTRAN
body: <p>I'm trying to numerically solve the double pendulum by Runge-Kutta fourth order.</p>
  
  <p>The differential equations to be solved are on the following page:</p>
  
  <p><a href="http://www.myphysicslab.com/dbl_pendulum.html" rel="nofollow">http://www.myphysicslab.com/dbl_pendulum.html</a></p>
  
  <p>There is even an animation showing the angles 1 and 2. I want to solve the equations that appear further down the page by Runge-Kutta method.</p>
  
  <p>Print these results to a file and then plot them using gnuplot. But, my program does not run and can not find the reason, I hope you can help me.</p>
  
  <p>My code in Fortran:</p>
  
  <pre><code>PROGRAM MAIN
  
  !In this program the double pendulum is solved by Rugen-Kutta method of order 4
  
  !variables required
  INTEGER, PARAMETER:: PREC = SELECTED_REAL_KIND(9,99)
  REAL(PREC), DIMENSION(:), ALLOCATABLE :: theta1, theta2, omega1, omega2, t
  REAL(PREC) :: LEN1, LEN2, M1, M2, G, PI, DT
  !The following variables are for the Runge-Kutta
  REAL(PREC) :: K11, K12, K13, K14, K21, K22, K23, K24
  REAL(PREC) :: K31, K32, K33, K34, K41, K42, K43, K44
  INTEGER :: I, J, N
  
  !Definition of variables needed
  PI= 4.*ATAN(1.)
  I=0
  G=9.8
  DT=0.0003 !change over time
  T(0)= 0.0 !initial time
  
  !number of repetitions
  n=500000
  ALLOCATE (theta1(0:n),omega1(0:n),theta2(0:n),omega2(0:n),t(0:n))
  
  !Initials values are given by the user
  PRINT*, 'First pendulum'
  CALL INICIALIZA(THETA1, OMEGA1, N, LEN1, M1)
  PRINT*, ' '
  PRINT*, 'Second pendulum'
  CALL INICIALIZA(THETA2, OMEGA2, N, LEN2, M2)
  
  !Runge-kutta method
  DO
  
  !Runge-Kutta fourth order
      !
      K11 = DT*OMEGA1(I)
      K21 = DT*OMEGA2(I)
      K31 = DT*DOM1(M1, M2, THETA1(I), THETA2(I) &amp; 
      , OMEGA1(I), OMEGA2(I), LEN1, LEN2, G)
      K41 = DT*DOM2(M1, M2, THETA1(I), THETA2(I) &amp; 
      , OMEGA1(I), OMEGA2(I), LEN1, LEN2, G)
      !
      K12 = DT*(OMEGA1(I) + K11/2.0)
      K22 = DT*(OMEGA2(I) + K21/2.0)
      K32 = DT*DOM1(M1, M2, (THETA1(I) + K11/2.0), (THETA2(I) + K21/2.0) &amp; 
      , (OMEGA1(I) + K31/2.0), (OMEGA2(I) + K41/2.0), LEN1, LEN2, G)
      K42 = DT*DOM2(M1, M2, (THETA1(I) + K11/2.0), (THETA2(I) + K21/2.0) &amp; 
      , (OMEGA1(I) + K31/2.0), (OMEGA2(I) + K41/2.0), LEN1, LEN2, G)
      !
      K13 = DT*(OMEGA1(I) + K12/2.0)
      K23 = DT*(OMEGA1(I) + K22/2.0)
      K33 = DT*DOM1(M1, M2, (THETA1(I) + K12/2.0), (THETA2(I) + K22/2.0) &amp; 
      , (OMEGA1(I) + K32/2.0), (OMEGA2(I) + K42/2.0), LEN1, LEN2, G)
      K43 = DT*DOM2(M1, M2, (THETA1(I) + K12/2.0), (THETA2(I) + K22/2.0) &amp; 
      , (OMEGA1(I) + K32/2.0), (OMEGA2(I) + K42/2.0), LEN1, LEN2, G)
      !
      K14 = DT*(OMEGA1(I) + K13)
      K24 = DT*(OMEGA1(I) + K23)
      K34 = DT*DOM1(M1, M2, (THETA1(I) + K13), (THETA2(I) + K23) &amp; 
      , (OMEGA1(I) + K33), (OMEGA2(I) + K43), LEN1, LEN2, G)
      K44 = DT*DOM2(M1, M2, (THETA1(I) + K13), (THETA2(I) + K23) &amp; 
      , (OMEGA1(I) + K33), (OMEGA2(I) + K43), LEN1, LEN2, G)
      !
      THETA1(I+1) = THETA1(I)+((K11+(2.0*(K12+K13))+K14)/6.0)
      THETA2(I+1) = THETA2(I)+((K21+(2.0*(K22+K23))+K24)/6.0)
      OMEGA1(I+1) = OMEGA1(I)+((K31+(2.0*(K32+K33))+K34)/6.0)
      OMEGA2(I+1) = OMEGA2(I)+((K41+(2.0*(K42+K43))+K44)/6.0)
      !
      if (theta1(i+1) &gt; PI ) theta1(i+1)=theta1(i+1)-2.*PI
      if (theta1(i+1) &lt; -PI) theta1(i+1)=theta1(i+1)+2.*PI
      if (theta2(i+1) &gt; PI ) theta2(i+1)=theta2(i+1)-2.*PI
      if (theta2(i+1) &lt; -PI) theta2(i+1)=theta2(i+1)+2.*PI    
  
      t(i+1) = t(i) + dt
  
      IF (i &gt;= n-1) EXIT
  
      i=i+1
  
  ENDDO
  !The results are saved in a file
  OPEN (UNIT=10,FILE='dou2.dat',STATUS='UNKNOWN')
  
  do j=0,n
      wRITE(10,*) theta1(j),theta2(j),omega1(j),omega2(j),t(j)
  end do
  
  CLOSE(10)
  
  END PROGRAM MAIN
  
  !Subroutine to initial values
  SUBROUTINE inicializa(theta, omega, n ,length, m)
  INTEGER, PARAMETER:: PREC = SELECTED_REAL_KIND(9,99)
  INTEGER, INTENT (IN):: n
  REAL(PREC), DIMENSION(0:n):: theta, omega
  REAL(PREC):: length, m
  print*,'starting angle'
  read*, theta(0)
  print*,'Initial angular velocity'
  read*, omega(0)
  print*,'Lenght of pendulum'
  read*, length
  print*, 'Mass pendulum'
  read*, m
  END SUBROUTINE inicializa
  
  !Functions that determine the derivative of omega (the angular velocity)
  
  REAL FUNCTION DOM1(N1, N2, X1, X2, Y1, Y2, L1, L2, A)
      INTEGER, PARAMETER:: PREC = SELECTED_REAL_KIND(9,99)
      REAL(PREC) :: N1, N2, X1, X2, Y1, Y2, L1, L2, A
      DOM1 = ((-A*((2.0*N1)+N2)*SIN(X1))-(N2*A*SIN(X1-(2.0*X2))) &amp; 
      -(2.0*SIN(X1-X2)*N2*((Y2*Y2*L2)+(Y1*Y1*L1*COS(X1-X2))))) &amp;
      /(L1*((2.0*N1)+N2-(N2*COS((2.0*X1)-(2.0*X2)))))
  END FUNCTION DOM1
  
  REAL FUNCTION DOM2(N1, N2, X1, X2, Y1, Y2, L1, L2, A)
      INTEGER, PARAMETER:: PREC = SELECTED_REAL_KIND(9,99)
      REAL(PREC) :: N1, N2, X1, X2, Y1, Y2, L1, L2, A
      DOM2 = (2.0*SIN(X1-X2)*((Y1*Y1*L1*(N1+N2))+(A*(N1+N2)*COS(X1)) &amp;
      +(Y2*Y2*L2*N2*COS(X1-X2))))/(L2*((2.0*N1)+N2 &amp;
      -N2*COS((2.0*X1)-(2.0*X2))))
  END FUNCTION DOM2
  </code></pre>
  
  <p>Here is a good description of the method: <a href="http://mathworld.wolfram.com/Runge-KuttaMethod.html" rel="nofollow">http://mathworld.wolfram.com/Runge-KuttaMethod.html</a></p>
  
  <p>The error that I get when running the program is:</p>
  
  <pre><code>Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
  Backtrace for this error:
     #0  0xB763D163
     #1  0xB763D800
     #2  0xB774A3FF
     #3  0x8048F17 in MAIN__ at doupend2.f90:?
   Violación de segmento
  </code></pre>