rk78.f

Go to the documentation of this file.

      SUBROUTINE rk78(DER,T,TOUT,N,X,H,RELERR,ABSERR,SP)  
C  VERSION OF 4/1/85
C  PURPOSE
C    RUNGE-KUTTA 7(8) METHOD AS GIVEN BY THE FEHLBERG COEFFICIENTS
C    NASA TR R-287
C  INPUT
C    DER    = NAME OF SUBROUTINE CONTAINING THE DERIVATIVES
C    T      = INITIAL INTEGRATION TIME
C    TOUT   = FINAL INTEGRATION TIME
C    N      = NUMBER OF EQUATIONS TO BE INTEGRATED
C    X      = THE STATE TO BE INTEGRATED, DIMENSION 6 HERE 
C    H      = INITIAL GUESS OF STEP-SIZE
C    RELERR = RELATIVE TOLERANCE
C    ABSERR = ABSOLUTE TOLERANCE
C    SP     = APPROXIMATE TIME THAT USEROP WILL BE CALLED
C  OUTPUT
C    T      = AT TOUT, T=TOUT
C    X      = STATE AT TOUT
C  CALL SUBROUTINES
C    DER, USEROP
C  REFERENCES
C    JPL EM 312/87-153, 20 APRIL 1987
C    NASA TR R-287, OCTOBER 1968
C  ANALYSIS
C    J. H. KWOK - JPL
C  PROGRAMMER
C    J. H. KWOK - JPL
C  PROGRAM MODIFICATIONS
C
C  DECLARED DER AS EXTERNAL TO ELIMINATE COMPILER WARNING.
C  F. S. PATT, JANUARY 13, 1995.
C
C    NONE
C  COMMENTS
C    USER MUST CALL RK78CN BEFORE CALLING RK78.
C    THE DIMENSIONS OF XDUM AND F CAN BE CHANGED FROM 6 TO ACCOMMODATE
C    ANY SYSTEM OF DIFFERENTIAL EQUATIONS.  USEROP IS A USER OPTION
C    ROUTINE WHERE THE HEALTH OF THE INTEGRATION CAN BE MONITORED BY
C    SETTING SP=ZERO.  IT CAN ALSO BE USED TO CHECK CERTAIN CONDITIONS
C    OF THE STATE, AND POSSIBLY RESTART WITH A NEW STATE OR STEP-SIZE
C    BY SETTING ISTART=1 IN USEROP.
C
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      common/felcon/ch(13),alph(13),beta(13,12),ordrcp,norder,ntimes
      dimension xdum(6),f(6,13),x(*)
      EXTERNAL der
      DATA zero,one/0.d0,1.d0/
C
C *** FACT IS A SCALING FACTOR TO REDUCE THE ESTIMATED STEPSIZE TO AVOID
C *** EXCEESIVE NUMBER OF REJECTION.
C *** SMALL IS A SMALL NUMBER COMPARED TO THE SIZE OF T
C
      DATA fact,small/0.7d0,1.d-8/
C
C *** INITIALIZATION
C
      nstp=0
      nrej=0
      sdt=dsign(one,tout-t)   
      dt=dabs(h)*sdt
      spp=dabs(sp)*sdt
      tpr=t
      istart=0
      iflag=0
C
C *** CALL USEROP AT INITIAL T
C
C      CALL USEROP(T,X,N,DT,NREJ,NSTP,ISTART,IFLAG)
C
C *** START INTEGRATION AND SAVE INITIAL INFO
C
  100 CONTINUE
      istart=0
      nstp=nstp+1
      tdum=t
      tpr=tpr+spp
      DO 110 i=1,n   
  110 xdum(i)=x(i)
C
C *** CHECK FOR T+DT PASSING TOUT
C
      IF (dabs(dt).GT.dabs(tout-t)) dt=tout-t  
C
C *** CHECK FOR REACHING TOUT
C
      IF (dabs(t-tout).LT.small) GO TO 900
C
C *** START FUNCTION EVALUATIONS
C
      CALL der(t,x,f(1,1))
      DO 300 k=2,ntimes  
      kk=k-1
      DO 310 i=1,n   
      temp=zero
      DO 320 j=1,kk  
  320 temp=temp+beta(k,j)*f(i,j)
  310 x(i)=xdum(i)+dt*temp  
      t=tdum+alph(k)*dt 
      CALL der(t,x,f(1,k))  
  300 CONTINUE
C
C *** COMPUTE NEW STATE
C
      DO 400 i=1,n  
      temp=zero 
      DO 410 l=1,ntimes 
  410 temp=temp+ch(l)*f(i,l)
  400 x(i)=xdum(i)+dt*temp
C
C *** START LOCAL TRUNCATION ERROR COMPUTATION
C
      err=relerr
      DO 500 i=1,n  
      ter=dabs((f(i,1)+f(i,11)-f(i,12)-f(i,13))*ch(12)*dt)
      tol=dabs(x(i))*relerr+abserr
      const=ter/tol
      IF (const.GT.err) err=const
  500 CONTINUE  
C
C *** ESTIMATE NEW STEP-SIZE
C
      dt=fact*dt*(one/err)**ordrcp
C
C *** IF IFLAG.EQ.1, ACCEPT THIS STEP UNCONDITIONALLY
C
      IF (iflag.EQ.1) THEN
        iflag=0
        GO TO 510
      ENDIF
C
C *** REJECT LAST STEP IF ER.LT.ONE OR ISTART.EQ.1
C
  220 CONTINUE
      IF (err.GT.one.OR.istart.EQ.1) THEN
        t=tdum
        DO 610 i=1,n
  610   x(i)=xdum(i)
        nrej=nrej+1
        nstp=nstp-1
        GO TO 100
      ENDIF
C
C *** CHECK FOR TIME TO CALL USEROP
C
  510 CONTINUE
      IF (sdt.GT.zero.AND.t.GE.tpr) GO TO 200
      IF (sdt.LT.zero.AND.t.LE.tpr) GO TO 200
      GO TO 100
  200 CONTINUE
C      CALL USEROP(T,X,N,DT,NREJ,NSTP,ISTART,IFLAG)
      IF (istart.EQ.1) GO TO 220
      GO TO 100
  900 CONTINUE
C
C *** CALL USEROP THE LAST TIME
C
C      CALL USEROP(T,X,N,DT,NREJ,NSTP,ISTART,IFLAG)
      RETURN
      END