!     Last change:  DOS  29 Jul 2000    4:46 pm
!                                       *** copyright 2000 ***
! *** filename false.f95                *** John F. Monahan **
!                                       **********************
      PROGRAM FALSE
!  EXAMPLE OF POOR CONVERGENCE OF REGULA FALSI/FALSE POSITION
!   REMEDIED BY ILLINOIS METHOD
!  THE FUNCTION IS G  -- START WITH EXTERNAL, SAYING THAT THE
!                        ARGUMENT PASSED IS A SUBPROGRAM
      IMPLICIT NONE
      REAL, EXTERNAL :: G
      REAL YL,YR
!
  21  FORMAT(' ESTIMATE OF ROOT=',F12.6,4X,'VALUE OF F THERE=',F12.6)
!
      OPEN( UNIT=6, FILE='false.out' )
!
!  DO 30 ITERATIONS OF REGULA FALSI SEARCH ON (YL,YR)
      YL = 0.
      YR = 2.
      CALL REGULA(G,YL,YR,30,.00001,.00001)                        
      WRITE(6,21) YL,YR                                                 
!
!  DO 30 ITERATIONS OF ILLINOIS SEARCH ON (YL,YR)
      YL = 0.
      YR = 2.
      CALL ILLINI(G,YL,YR,30,.00001,.00001)
      WRITE(6,21) YL,YR
      STOP                                                              
      END PROGRAM FALSE
      REAL FUNCTION G(X)
!  FUNCTION FOR WHICH ROOT IS TO BE FOUND
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      G = X*X*X - 1. 
      RETURN                                                            
      END FUNCTION G
      SUBROUTINE REGULA(F,XL,XR,MIT,EPSX,EPSF)
!  USES REGULA FALSI, SECANT FORMULA BUT USE LAST 2 PTS STRADDLING ROOT
!  RETURNS ROOT IN XL, VALUE OF F THERE IN XR
!  EPSF: STOP IF ABS(F(X)) .LE. EPSF
!  EPSX: STOP IF CHANGE IN X IS .LE. EPSX
      IMPLICIT NONE
      REAL F
      REAL, INTENT(IN OUT) :: XL,XR
      REAL, INTENT(IN) :: EPSX,EPSF
      INTEGER, INTENT(IN) :: MIT
      REAL ESPX,ESPF,FL,S,FR,XN,FN
      INTEGER IT
!                   SET LOWER LIMIT ON EPSILONS
      ESPX = MAX(.000001,EPSX)
      ESPF = MAX(.000001,EPSF)
!                   GET DIRECTION
      FL=F(XL)
      S=1.                                                              
!  INCREASING (S=1) OR DECREASING (S=-1) ?
      IF( FL .GT. 0. ) S = -1.
      FR=F(XR)
!                   START LOOP
      DO IT = 1,MIT
!                   TEST FOR CONVERGENCE
      IF(XR-XL .LT. ESPX) EXIT
!                   FIND NEW POINT
      XN = XL - FL*(XL-XR)/(FL-FR)
!                   EVALUATE AT NEW POINT
      FN = F(XN)
      write(6,21) IT,XN,FN,XL,XR
  21  format(' REGULA IT=',I2,' XNEW=',F8.5,' F(XN)=',F9.6,' INT IS ('&
     &,F8.5,' , ',F8.5,' )')
!                   STOP IF ROOT IS CLOSE ENOUGH
      IF( ABS(FN) .LT. ESPF ) EXIT
!
!                   BRANCH
      IF( FN*S .LT. 0. ) THEN
      XL=XN
      FL=FN                                                             
                         ELSE
      XR=XN
      FR=FN                                                             
                         ENDIF
      END DO  ! LOOP ON IT
!                   MOVE ROOT TO XL, FUNCTION VAL TO XR
      XL=XN
      XR=FN                                                             
      RETURN                                                            
      END SUBROUTINE REGULA
      SUBROUTINE ILLINI(F,XL,XR,MIT,EPSX,EPSF)
!  USES 'ILLINOIS METHOD' OF ROOT-FINDING, A VARIANT OF REGULA FALSI
!  RETURNS ROOT IN XL, VALUE OF F THERE IN XR
!  EPSF: STOP IF ABS(F(X)) .LE. EPSF
!  EPSX: STOP IF CHANGE IN X IS .LE. EPSX
      IMPLICIT NONE
      REAL F
      REAL, INTENT(IN OUT) :: XL,XR
      REAL, INTENT(IN) :: EPSX,EPSF
      INTEGER, INTENT(IN) :: MIT
      REAL ESPX,ESPF,FL,S,FR,XN,FN
      INTEGER IT,LAST
!                   SET LOWER LIMIT ON EPSILONS
      ESPX = MAX(.000001,EPSX)
      ESPF = MAX(.000001,EPSF)
!                   GET DIRECTION
      FL=F(XL)
      S=1.
!  INCREASING (S=1) OR DECREASING (S=-1) ?
      IF( FL .GT. 0. ) S = -1.
      FR=F(XR)
      LAST = 2                ! FIRST TIME THROUGH
!                   START LOOP
      DO IT = 1,MIT
!                   TEST FOR CONVERGENCE
      IF(XR-XL .LT. ESPX) EXIT
!                   FIND NEW POINT
      XN = XL - FL*(XL-XR)/(FL-FR)
!                   EVALUATE AT NEW POINT
      FN = F(XN)
      write(6,21) IT,XN,FN,XL,XR
  21  format(' ILLINI IT=',I2,' XNEW=',F8.5,' F(XN)=',F9.6,' INT IS ('&
     &,F8.5,' , ',F8.5,' )')
!                   STOP IF ROOT IS CLOSE ENOUGH
      IF( ABS(FN) .LT. ESPF ) EXIT
!
!                   BRANCH
      IF( FN*S .LT. 0. ) THEN
      XL=XN
      FL=FN                                                             
!  IF MOVED THIS POINT LAST TIME, DECREASE OTHER POINT BY HALF
      IF( LAST .EQ. 1 ) FR = FR/2.
      LAST = 1                ! LAST = 1 MOVED LEFT ENDPOINT LAST TIME
                         ELSE
      XR=XN
      FR=FN                                                             
!  IF MOVED THIS POINT LAST TIME, DECREASE OTHER POINT BY HALF
      IF(LAST.EQ.0) FL=FL/2.                                            
      LAST=0                  ! LAST = 0 MOVED RIGHT ENDPOINT LAST TIME
                         ENDIF
      END DO  ! LOOP ON IT
!                   MOVE ROOT TO XL, FUNCTION VAL TO XR
      XL=XN
      XR=FN                                                             
      RETURN                                                            
      END SUBROUTINE ILLINI
! *** end of filename false.f95         *********************