我正在尝试用Runge-Kutta四阶的数字求解双摆。
要求解的微分方程在以下页面上：
http://www.myphysicslab.com/dbl_pendulum.html
甚至还有一个动画显示角度 1 和 2。我想通过 Runge-Kutta 方法求解页面下方出现的方程。
将这些结果打印到一个文件中，然后使用 gnuplot 绘制它们。但是，我的程序无法运行，找不到原因，希望您能帮助我。
我在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
以下是该方法的一个很好的描述：http://mathworld.wolfram.com/Runge-KuttaMethod.html
运行程序时出现的错误是：
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
双摆的力学问题和通过龙格-库塔的四阶数值方法对运动方程初始值进行数值积分。选择的语言是FORTRAN和gfortran编译器。这些图是在 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
图形是在 gnuplot 的帮助下通过.dat文件生成的。图1是xy平面中相应质量m1和m2的位置图。图 2 由配置空间中的广义坐标（角度 1 和角度 2）图组成。