C
C <<<<<<<<<<<<<<<<<<<<<< PENDEQN - PENDEQN - PENDEQN >>>>>>>>>>>>>>>>>>>>>>
C
SUBROUTINE PENDEQN(T, Y, DYDT, NEQN, ERRMSG)
C
C Calculate the derivative of Y at T for NEQN equations and store them in
C DYDT. Note that this subroutine has been set up for use with the pendulum
C ODE problem.
C
IMPLICIT REAL*8 (A-H, O-Z)
DIMENSION Y(NEQN), DYDT(NEQN)
CHARACTER ERRMSG*80
C
COMMON /PEND/ GRAV, RLENGTH
C
C Check input.
C
IF (RLENGTH .EQ. 0.D0) THEN
ERRMSG = 'RLENGTH has a value of zero in PENDEQN.'
RETURN
ENDIF
IF (NEQN .NE. 2) THEN
ERRMSG = 'NEQN must equal 2 in PENDEQN.'
RETURN
ENDIF
C
C Note that our 2 differential equations are
C d-theta/dt = omega = DYDT(1)
C d-omega/dt = -(grav/rlength) * theta = DYDT(2)
C Set up the initial values for DYDT and for Y:
C Y(1) = theta(i)
C Y(2) = omega(i)
C
DYDT(1) = Y(2)
DYDT(2) = -(GRAV/RLENGTH) * Y(1)
C
RETURN
END
C
C <<<<<<<<<<<<<< RK4_PENDULUM - RK4_PENDULUM - RK4_PENDULUM >>>>>>>>>>>>>>>>
C
SUBROUTINE RK4_PENDULUM(Y, DYDT, NEQN, T, H, YOUT, PENDEQN,
+ ERRMSG)
C
C Solve up to 10 ordinary differential equations using 4th-order Runge-Kutta.
C
C Variable definitions:
C
C Y: Dependent variable array of size NEQN. [THETA(I), OMEGA(I)]
C DYDT: Value of time derivatives of our dependent variables.
C [d-THETA(I)/dt, d-OMEGA(I)/dt]
C NEQN: The number of equations being solve which is also the number
C of dependent variables stored in Y.
C T: The current value of the independent variable. [TIME(I)]
C H: The step size. [TSTEP]
C YOUT: The new values of Y (NEQN element array).
C ERMSG: A string containing any errors that are discovered in the
C calculations.
C
C The user supplies the subroutine PENDEQN(T, Y, DYDT, NEQN, SUBR, ERRMSG)
C which returns derivatives DYDT at T.
C
IMPLICIT REAL*8 (A-H, O-Z)
C
EXTERNAL PENDEQN
C
C Set the maximum number of functions.
C
PARAMETER (NMAX=10)
C
DIMENSION Y(NEQN), DYDT(NEQN), YOUT(NEQN)
DIMENSION YT(NMAX), DYT(NMAX), DYM(NMAX)
CHARACTER ERRMSG*80
C
C Check number of functions that were passed.
C
IF (NEQN .GT. NMAX) THEN
ERRMSG = 'Number of functions passed exceeds NMAX in '
1 // 'RK4_PENDULUM.'
RETURN
ENDIF
C
HH = H * 0.5D0
H6 = H / 6.D0
TH = T + HH
C
C First step.
C
DO 11 I = 1, NEQN
YT(I) = Y(I) + HH*DYDT(I)
11 CONTINUE
C
C Second step.
C
CALL PENDEQN(TH, YT, DYM, NEQN, ERRMSG)
IF (ERRMSG(1:1) .NE. ' ') RETURN
DO 12 I = 1, NEQN
YT(I) = Y(I) + HH*DYM(I)
12 CONTINUE
C
C Third step.
C
CALL PENDEQN(TH, YT, DYM, NEQN, ERRMSG)
IF (ERRMSG(1:1) .NE. ' ') RETURN
DO 13 I = 1, NEQN
YT(I) = Y(I) + H*DYM(I)
DYM(I) = DYT(I) + DYM(I)
13 CONTINUE
C
C Fourth step.
C
CALL PENDEQN(T+H, YT, DYT, NEQN, ERRMSG)
IF (ERRMSG(1:1) .NE. ' ') RETURN
DO 14 I = 1, NEQN
YOUT(I) = Y(I) + H6*(DYDT(I) + DYT(I) + 2.D0*DYM(I))
14 CONTINUE
C
RETURN
END