shoot.for

C
C*****************************  ABSTRACT  *******************************
C
C     THIS PROGRAM SOLVES THE BVP POSED IN EXAMPLE 5.6 USING THE SHOOTING
C  METHOD.  THE BLACK BOX FUNCTION CREATED BY THE INTEGRATION OF THE BVP
C  TRANSFORMED INTO AN IVP IS SOLVED USING THE SECANT METHOD.
C
C
C****************************  NOMENCLATURE  ****************************
C
C  DXMAX- THE MAXIMUM ALLOWED CHANGE IN X BETWEEN ITERATIONS
C  ERLIM- CONVERGENCE CRITERIA
C  F- THE NAME OF THE SUBROUTINE THAT CALCULATES THE VALUE OF THE BLACK
C     FUNCTION
C  ICMAX- THE MAXIMUM NUMBER OF ITERATIONS ALLOWED
C  IPRINT- PRINT FLAG =1 PRINT INTERMEDIATE VALUES; =0 NO PRINT
C  X1-	THE INITIAL VALUE THE ROOT OF THE NONLINEAR EQUATION
C
C************************************************************************
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Y(2),CAY(2),Y1(2),DYDX(2)
      COMMON /ONE/RCAY,D
      EXTERNAL F
C
C  SET SECANT PARAMETERS
C
      ERLIM=1.D-6
      ICMAX=50
      DXMAX=2.
      IPRINT=1
C
C  SET INITIAL GUESS FOR DY/DX AT X=0
C
      X1=1.
C
C  CALL SECANT METHOD
C
      CALL SECANT(F,ERLIM,ICMAX,DXMAX,IPRINT,X1)
      STOP
      END
C
C*****************************  ABSTRACT  *****************************
C
C     THIS SUBROUTINE CALLS SUBROUTINE RKUTTA TO INTEGRATE THE BVP
C  WHICH HAS BEEN TRANSFORMED INTO AN IVP.
C
C****************************  NOMENCLATURE  **************************
C
C  DXPNT- THE PRINT INTERVAL
C  FUNCT- THE NAME OF THE SUBROUTINE THAT CALCULATES F(X,Y)
C  IPNT- PRINT FLAG (=0 NO PRINT; =1 PRINT)
C  N- THE NUMBER OF DEPENDENT VARIABLES
C  P- THE MAXIMUM PERCENTAGE CHANGE IN A DEPENDENT VARIABLE
C  X- THE INDEPENDENT VARIABLE
C  XMAX-  THE FINAL VALUE OF X
C  Y(I)-  THE ITH DEPENDENT VARIABLE
C
C**********************************************************************
C
      SUBROUTINE F(XV,FV)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Y(2),Y1(2),CAY(2),DYDX(2)
      EXTERNAL FUNCT
C
C  SPECIFY INPUT DATA FOR RUNGE-KUTTA INTEGRATOR
C
      N=2
      XMAX=1.
      DXPNT=.1
      IPNT=1
      P=10.
C  SET INITIAL CONDITIONS
      X=0.0
      Y(1)=0.0
      Y(2)=XV
C
C  CALL RUNGE-KUTTA INTEGRATOR
C
      CALL RKUTTA(FUNCT,N,P,X,Y,XMAX,DXPNT,IPNT,CAY,Y1,DYDX)
C
C  SET FUNCTION VALUE AS VALUE OF DYDX AT X=1
C
      FV=Y(1)-1.
      RETURN
      END
C
C*************************  ABSTRACT  ********************************
C
C     THIS SUBROUTINE CALCULATES THE DERIVATIVE OF Y(I) WITH RESPECT TO
C  X, DYDX(I).  THE USER SUPPLIES THE EQUATIONS FOR DYDX(I).
C
C*********************************************************************
C
C
      SUBROUTINE FUNCT(X,Y,DYDX)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Y(2),DYDX(2)
      COMMON /ONE/RCAY,D
      DYDX(2)=-Y(1)*Y(2)-5.
   10 DYDX(1)=Y(2)
      RETURN
      END
C
C*****************************  ABSTRACT  *****************************
C
C     THIS SUBROUTINE CALLS SUBROUTINE RK1 AT EACH PRINT INTERVAL
C  AND PRINTS OUT THE INTERMEDIATE RESULTS IF IPNT=1.
C
C****************************  NOMENCLATURE  **************************
C
C  DXPNT- THE PRINT INTERVAL
C  DXP- THE INTERNALLY CALCULATED PRINT INTERVAL
C  F- THE NAME OF THE SUBROUTINE THAT CALCULATES F(X,Y)
C  IPNT- PRINT FLAG (=0 NO PRINT; =1 PRINT)
C  N- THE NUMBER OF DEPENDENT VARIABLES
C  NS- THE NUMBER OF PRINT INTERVAL STEPS
C  PRO- THE PRINT INTERVAL NEAR XMAX
C  X- THE INDEPENDENT VARIABLE
C  XF- THE VALUE OF X AT THE END OF THE PRINT INTERVAL
C  XO- THE VALUE OF X AT THE BEGINNING OF THE PRINT INTERVAL
C  XMAX-  THE FINAL VALUE OF X
C  Y(I)-  THE ITH DEPENDENT VARIABLE
C
C**********************************************************************
C
      SUBROUTINE RKUTTA(F,N,P,X,Y,XMAX,DXPNT,IPNT,CAY,Y1,DYDX)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Y(1),DYDX(1),CAY(1),Y1(1)
      EXTERNAL F
C  PRINT INITIAL CONDITIONS IF IPNT = 1
      IF(IPNT.EQ.1)WRITE(6,10)X,(Y(J),J=1,N)
C  CALCULATE THE NUMBER OF STEPS PER DXPNT
      DXP=DXPNT
      NS=IFIX(XMAX/DXPNT)+1
      XO=X
      PRO=XMAX-DXPNT*DFLOAT(NS-1)
C
C  CALL SUBROUTINE RK TIL X=XMAX
C
      DO 1000 I=1,NS
      IF(I.EQ.NS)DXP=PRO
      IF((I.EQ.NS).AND.(PRO.LT.1.E-6))GO TO 1000
      XF=XO+DXP
C  CALL SUBROUTINE RK
      CALL RK(F,N,P,XO,XF,Y,CAY,Y1,DYDX)
C  PRINT INTERMEDIATE RESULTS IF IPNT=1
      IF(IPNT.EQ.1)WRITE(6,10)XF,(Y(J),J=1,N)
   10 FORMAT( 10X,3H X=,D12.5,3X,3H Y=,5(D14.7,3X))
C  INCREMENT X
 1000 XO=XO+DXP
      RETURN
      END
C
C*************************  ABSTRACT  ********************************
C
C     THIS APPLIES A FOURTH ORDER RUNGE-KUTTA METHOD FOR THE INTEGRATION
C  OF A SYSTEM OF ODE'S FROM X=XO TO X=XF.
C
C****************************  NOMENCLATURE  **************************
C
C  CAY(I)- A VECTOR USED TO STORE THE SLOPES OF DEPENDENT VARIABLE I
C	   AT THE BEGINNING, MIDDLE AND END OF A STEP DX
C  DX- THE STEP SIZE FOR X
X  DXM- THE MINIMUM STEP SIZE DX
C  DYDX(I)- THE DERIVATIVE OF Y(I) WITH RESPECT TO X
C  F- THE NAME OF THE SUBROUTINE THAT CALCULATES F(X,Y)
C  IPNT- PRINT FLAG (=0 NO PRINT; =1 PRINT)
C  N- THE NUMBER OF DEPENDENT VARIABLES
C  NS- THE NUMBER OF PRINT INTERVAL STEPS
C  PRO- THE PRINT INTERVAL NEAR XMAX
C  X- THE INDEPENDENT VARIABLE
C  XF- THE VALUE OF X AT THE END OF THE PRINT INTERVAL
C  XO- THE VALUE OF X AT THE BEGINNING OF THE PRINT INTERVAL
C  Y(I)- THE ITH DEPENDENT VARIABLE
C  Y1(I)- AN INTERNALLY USED VALUE OF THE ITH DEPENDENT VARIABLE
C
C**********************************************************************
C
      SUBROUTINE RK(F,N,P,XO,XF,Y,CAY,Y1,DYDX)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Y(1),DYDX(1),CAY(1),Y1(1)
      EXTERNAL F
      X=XO
C
C  BEGIN ITERATIVE CALCULATION LOOP
C
    1 CALL F(N,X,Y,DYDX)
C
C  STEP SIZE CALCULATION
C
      DDD=DYDX(1)
      IF(DABS(DYDX(1)).LT.1.E-10)DDD=1.E-10
      DXM=DABS(P*Y(1)/DDD/100.)
C  DETERMINE THE MINIMUM STEP SIZE DX BASED UPON A P% CHANGE
      DO 2 I=2,N
      IF(DYDX(I).LT.1.E-10)GO TO 2
      TEST=DABS(P*Y(I)/DYDX(I)/100.)
      IF(TEST.LT.DXM)DXM=TEST
    2 CONTINUE
C  MAKE SURE THAT X DOES NOT EXCEED XF
      IF(DXM.LT.1.E-4)DXM=1.E-4
      DX=DXM
      XT=X+DX
      IF(XT.GT.XF)DX=XF-X
C  SET K1
      DO 33 I=1,N
      CAY(I)=DYDX(I)
   33 Y1(I)=Y(I)+DYDX(I)*DX/2.
      X1=X+DX/2.
      CALL F(N,X1,Y1,DYDX)
C  SET K2
      DO 34 I=1,N
      CAY(I)=CAY(I)+2.*DYDX(I)
   34 Y1(I)=Y(I)+DYDX(I)*DX/2.
      CALL F(N,X1,Y1,DYDX)
C  SET K3
      DO 35 I=1,N
      CAY(I)=CAY(I)+2.*DYDX(I)
   35 Y1(I)=Y(I)+DYDX(I)*DX
      X1=X+DX
      CALL F(N,X1,Y1,DYDX)
C  SET K4 AND CALCULATE NEW VALUES OF Y(I)
      DO 36 I=1,N
      CAY(I)=CAY(I)+DYDX(I)
   36 Y(I)=Y(I)+DX*CAY(I)/6.
C  INCREMENT X
  100 X=X+DX
      IF(X.LT.XF)GO TO 1
      RETURN
      END
C
C*****************************  ABSTRACT  *******************************
C
C     THIS SUBROUTINE APPLIES THE SECANT METHOD FOR THE SOLUTION OF A
C  SINGLE NONLINEAR EQUATION.
C
C****************************  NOMENCLATURE  ****************************
C
C  DX- THE CALCULATED CHANGE IN THE VALUE OF X
C  DXMAX- THE MAXIMUM ALLOWED CHANGE IN X BETWEEN ITERATIONS
C	  (USER SPECIFIED)
C  ERLIM- CONVERGENCE CRITERIA
C  FO- THE NEXT TO THE LATEST FUNCTION VALUE
C  F1- THE LATEST FUNCTION VALUE
C  FV- THE FUNCTION VALUE
C  IC- THE ITERATION COUNTER
C  ICMAX- THE MAXIMUM NUMBER OF ITERATIONS ALLOWED
C  IPRINT- PRINT FLAG =1 PRINT INTERMEDIATE VALUES; =0 NO PRINT
C  X1-	THE LATEST VALUE THE ROOT OF THE NONLINEAR EQUATION
C  XO-	THE NEXT TO THE LATEST VALUE THE ROOT OF THE NONLINEAR EQUATION
C
C************************************************************************
C
      SUBROUTINE SECANT(F,ERLIM,ICMAX,DXMAX,IPRINT,X1)
      IMPLICIT REAL*8(A-H,O-Z)
      EXTERNAL F
      IC=0
C
C  BEGIN ITERATIVE LOOP
C
   10 IC=IC+1
C  CALL SUBRUTINE F
      CALL F(X1,F1)
C  INITIALIZE SECANT METHOD
      IF(IC.EQ.1)XO=X1*.95
      IF(IC.EQ.1)CALL F(XO,FO)
C  APPLY SECANT APPROXIATION
      DX=-F1*(X1-XO)/(F1-FO)
      IF(ABS(DX).GT.DXMAX)DX=DXMAX*ABS(DX)/DX
      XC=X1+DX
C  TRANSFER X1,F1 TO XO,FO
      XO=X1
      FO=F1
      X1=XC
C  PRINT OUT INTERMEDIATE RESULTS
      IF(IPRINT.EQ.1)WRITE(6,11)IC,XC,F1
   11 FORMAT( 5X,' IC=',I3,5X,' X=',F7.5,5X,' F(X)=',D10.3)
C  CHECK FOR EXCESSIVE ITERATIONS
      IF(IC.GT.ICMAX)WRITE(6,22)
   22 FORMAT( ' SECANT METHOD DID NOT CONVERGE')
      IF(IC.GT.ICMAX)STOP
C  CHECK FOR CONVERGENCE
      IF(ABS(DX).GT.ERLIM)GO TO 10
C
C  CONVERGED SOLUTION OBTAINED; RETURNING TO CALLING PROGRAM
C
      CALL F(XC,FV)
      RETURN
      END