Segmentation fault in double pendulum code by Runge-Kutta

Question

I'm trying to numerically solve the double pendulum by Runge-Kutta fourth order.

The differential equations to be solved are on the following page:

http://www.myphysicslab.com/dbl_pendulum.html

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.

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.

My code in Fortran:

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) & 
    , OMEGA1(I), OMEGA2(I), LEN1, LEN2, G)
    K41 = DT*DOM2(M1, M2, THETA1(I), THETA2(I) & 
    , 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) & 
    , (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) & 
    , (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) & 
    , (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) & 
    , (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) & 
    , (OMEGA1(I) + K33), (OMEGA2(I) + K43), LEN1, LEN2, G)
    K44 = DT*DOM2(M1, M2, (THETA1(I) + K13), (THETA2(I) + K23) & 
    , (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) > PI ) theta1(i+1)=theta1(i+1)-2.*PI
    if (theta1(i+1) < -PI) theta1(i+1)=theta1(i+1)+2.*PI
    if (theta2(i+1) > PI ) theta2(i+1)=theta2(i+1)-2.*PI
    if (theta2(i+1) < -PI) theta2(i+1)=theta2(i+1)+2.*PI    

    t(i+1) = t(i) + dt

    IF (i >= 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))) & 
    -(2.0*SIN(X1-X2)*N2*((Y2*Y2*L2)+(Y1*Y1*L1*COS(X1-X2))))) &
    /(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)) &
    +(Y2*Y2*L2*N2*COS(X1-X2))))/(L2*((2.0*N1)+N2 &
    -N2*COS((2.0*X1)-(2.0*X2))))
END FUNCTION DOM2

Here is a good description of the method: http://mathworld.wolfram.com/Runge-KuttaMethod.html

The error that I get when running the program is:

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

Answer 1

The mechanical problem of the double pendulum and numerical integration of the initial value of motion equations via the 4th order numerical method of Runge-Kutta. The chosen language was FORTRAN together with the gfortran compiler. The graphs were generated in gnuplot.

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
!         FEDERAL UNIVERSITY OF GOIÁS (UFG)         !
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
! Author            : VY Canis Majoris              !
! Theme             : Double Pendulum               !
! Language          : Fortran 90                    !
! Compiler          : Gfortran                      !
! Operating System  : Windows 10                    !
! Creation date     : 02/05/2023                    !
! Modification Date : 02/13/2023                    !
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!

program Pendulo_Duplo
    implicit none

    integer, parameter :: n = 100000
    real               :: h = 1.0E-3
    real               :: t0,o1_0,o2_0,w1_0,w2_0,m1,m2,l1,l2,g, eps
    real               :: t(n),o1(n),o2(n),w1(n),w2(n)! saída
    real               :: theta1(n),theta2(n),wtheta1(n),wtheta2(n)! saída
    

    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
    !  Enter the double pendulum attributes and initial conditions !
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
      t0 = 00.000
    o1_0 = -45.000
    o2_0 = 45.000
    w1_0 = 00.000
    w2_0 = 00.000
      m1 = 01.000
      m2 = 01.000
      l1 = 01.000
      l2 = 01.000
       g = 09.810


    call RK4(t0,o1_0,o2_0,w1_0,w2_0,h,n,m1,m2,l1,l2,g,t,o1,o2,w1,w2)

    stop
end program Pendulo_Duplo

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
! Differential equations of the double pendulum !
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!

! FIRST ODE
real function f1(t, o1, o2, w1, w2)
    implicit none
    real, intent(in) :: t, o1, o2, w1, w2
    f1 = w1
end function f1

! SECOND ODE
real function f2(t, o1, o2, w1, w2)
    implicit none
    real, intent(in) :: t, o1, o2, w1, w2
    f2 = w2
end function f2

! THIRD ODE
real function f3(t, o1, o2, w1, w2, m1, m2, l1, l2, g)
    implicit none
    real, intent(in) :: t, o1, o2, w1, w2
    real, intent(in) :: m1, m2, l1, l2, g

    f3 = (-G * (2 * M1 + M2) * sin(o1)&
                - M2 * G * sin(o1 - 2 * o2)&
                - 2 * sin(o1 - o2) * M2 * ((w2**2) * L2&
                + (w1**2) * L1 * cos(o1 - o2))) / (L1 * (2 * M1 + M2&
                - M2 * cos(2 * o1 - 2 * o2)))
end function f3

! FOURTH ODE
real function f4(t, o1, o2, w1, w2, m1, m2, l1, l2, g)
    implicit none
    real, intent(in) :: t, o1, o2, w1, w2
    real, intent(in) :: m1, m2, l1, l2, g

    f4 = (2 * sin(o1 - o2) * ((w1**2) * L1 * (M1&
                + M2) + G * (M1 + M2) * cos(o1)&
                + (w2**2) * L2 * M2 * cos(o1 - o2))) / (L2 * (2 * M1&
                + M2 - M2 * cos(2 * o1 - 2 * o2)))
end function f4

! CHANGE FUNCTION FROM DEGREES TO RADIANS
real function radianos(graus)
    implicit none
    real :: graus
    real(Kind = 8) :: pi
    pi = 4.0*atan(1.0d0)
    radianos = (pi/180)*graus
end function radianos

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
!  Sub Routine for the 4 Order Runge-Kutta method  !
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
subroutine RK4(t0,o1_0,o2_0,w1_0,w2_0,h,n,m1,m2,l1,l2,g,t,o1,o2,w1,w2)
    implicit none
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
    !    INPUT VARIABLES (INPUTS)   !
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
    integer, intent(in) :: n
    real,    intent(in) :: h
    real,    intent(in) :: t0, o1_0, o2_0, w1_0, w2_0
    real,    intent(in) :: m1, m2, l1, l2, g
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
    !   RETURN VARIABLES (OUTPUTS)   !
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
    real, intent(out)   :: t(n), o1(n), o2(n), w1(n), w2(n)
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
    ! AUXILIARY VARIABLES - INTERNAL STATEMENTS !
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
    integer        :: i
    real           :: k11, k21, k31, k41
    real           :: k12, k22, k32, k42
    real           :: k13, k23, k33, k43
    real           :: k14, k24, k34, k44
    real, external ::  f1,  f2,  f3,  f4, radianos
    real           :: x1(n), x2(n), y1(n), y2(n)

    !%%%%%%%%%%%%%%%%%%%%!
    ! Initial conditions !
    !%%%%%%%%%%%%%%%%%%%%!
     t(1) = t0
    o1(1) = radianos(o1_0)
    o2(1) = radianos(o2_0)
    w1(1) = w1_0
    w2(1) = w2_0

    do i = 1,n-1
        k11 = h*f1(t(i),o1(i),o2(i),w1(i),w2(i))
        k12 = h*f2(t(i),o1(i),o2(i),w1(i),w2(i))
        k13 = h*f3(t(i),o1(i),o2(i),w1(i),w2(i),m1,m2,l1,l2,g)
        k14 = h*f4(t(i),o1(i),o2(i),w1(i),w2(i),m1,m2,l1,l2,g)

        k21 = h*f1(t(i)+h/2,o1(i)+k11/2,o2(i)+k12/2,w1(i)+k13/2,w2(i)+k14/2)
        k22 = h*f2(t(i)+h/2,o1(i)+k11/2,o2(i)+k12/2,w1(i)+k13/2,w2(i)+k14/2)
        k23 = h*f3(t(i)+h/2,o1(i)+k11/2,o2(i)+k12/2,w1(i)+k13/2,w2(i)+k14/2,m1,m2,l1,l2,g)
        k24 = h*f4(t(i)+h/2,o1(i)+k11/2,o2(i)+k12/2,w1(i)+k13/2,w2(i)+k14/2,m1,m2,l1,l2,g)

        k31 = h*f1(t(i)+h/2,o1(i)+k21/2,o2(i)+k22/2,w1(i)+k23/2,w2(i)+k24/2)
        k32 = h*f2(t(i)+h/2,o1(i)+k21/2,o2(i)+k22/2,w1(i)+k23/2,w2(i)+k24/2)
        k33 = h*f3(t(i)+h/2,o1(i)+k21/2,o2(i)+k22/2,w1(i)+k23/2,w2(i)+k24/2,m1,m2,l1,l2,g)
        k34 = h*f4(t(i)+h/2,o1(i)+k21/2,o2(i)+k22/2,w1(i)+k23/2,w2(i)+k24/2,m1,m2,l1,l2,g)

        k41 = h*f1(t(i)+h,o1(i)+k31,o2(i)+k32,w1(i)+k33,w2(i)+k34)
        k42 = h*f2(t(i)+h,o1(i)+k31,o2(i)+k32,w1(i)+k33,w2(i)+k34)
        k43 = h*f3(t(i)+h,o1(i)+k31,o2(i)+k32,w1(i)+k33,w2(i)+k34,m1,m2,l1,l2,g)
        k44 = h*f4(t(i)+h,o1(i)+k31,o2(i)+k32,w1(i)+k33,w2(i)+k34,m1,m2,l1,l2,g)
   
        o1(i+1) = o1(i) + ((k11+2*(k21+k31)+k41)/6)
        o2(i+1) = o2(i) + ((k12+2*(k22+k32)+k42)/6)
        w1(i+1) = w1(i) + ((k13+2*(k23+k33)+k43)/6)
        w2(i+1) = w2(i) + ((k14+2*(k24+k34)+k44)/6)
         t(i+1) =  t(i) + h

        !write(*,*) t(i),o1(i),o2(i),w1(i),w2(i)
    end do

    ! Generation of .dat files !
    open(0,file = "pd0.dat")
        open(1,file = "pd1.dat")
            do i = 1,n-1
                
                ! Position of mass m1 in the xy plane !
                x1(i) =  l1*sin(o1(i))
                y1(i) = -l1*cos(o1(i))

                ! Position of mass m2 in the xy plane !
                x2(i) =  l1*sin(o1(i))+l2*sin(o2(i))
                y2(i) = -l1*cos(o1(i))-l2*cos(o2(i))

                !write(*,*) x1(i),y1(i),x2(i),y2(i)

                !           1    2     3     4      5      6
                write(0,*) 0.0, 0.0, x1(i), y1(i), x2(i), y2(i)
            end do

            do i = 1,n-1
                !           1      2      3     4      5
                write(1,*) t(i), o1(i), o2(i), w1(i), w2(i)
            end do
        close(1)
    close(0)

end subroutine RK4

The graphics were generated through .dat files with the help of gnuplot. Figure1 is a plot of the positions of the respective masses m1 and m2 in the xy plane. Figure2 consists of the plot of the generalized coordinates (angle1 and angle 2) in the configuration space.