!
! =====> Program - P125.F90
!
PROGRAM Rooter
! methods.
! Function types!
REAL Bisect, NwtRph, Secant
REAL Tol, X1, X2
PRINT *, 'This is Program >> P125 = Root finding demo'
!
! Tell program where data for READ * is coming from
OPEN(UNIT=5, FILE='P125.DAT') ! UNIT=5 is the default input
!
PRINT *, 'Program to demonstrate numerical root finding'
PRINT *, 'F = x**3 - Exp(-x) for example'
PRINT *
PRINT *, 'Enter tolerance for roots: '
READ *, Tol
Print *, Tol
PRINT *
PRINT *, 'Enter low limit, high limit: '
READ *, X1, X2
Print *, X1, X2
PRINT *
PRINT *, 'Bisection: ', Bisect(X1, X2, Tol)
PRINT *, 'Enter estimated root: '
READ *, X1
Print *, X1
PRINT *
PRINT *, 'Newton: ', NwtRph(X1, Tol)
PRINT *, 'Enter 1st, 2nd root estimates: '
READ *, X1, X2
Print *, X1, X2
PRINT *
PRINT *, 'Secant: ', Secant(X1, X2, Tol)
STOP
END
! SUPPORT ROUTINES
REAL FUNCTION F(x)
! The function of which a root (zero) is to be found.
F = x**3 - Exp(-x)
RETURN
END
REAL FUNCTION FPrime(x)
! The derivative of F.
FPrime = 3*x*x + Exp(-x)
RETURN
END
SUBROUTINE Swap(x1, x2)
! Exchanges two real values x1 and x2.
REAL x1, x2, Temp
Temp = x1
x1 = x2
x2 = Temp
RETURN
END
! ALGORITHMS
REAL FUNCTION Bisect(x_a, x_b, Tol)
! Finds one real root of the function F(x)
! on the interval [x_a, x_b].
REAL x_a, x_b, Tol, &
x_Lo, x_Hi, x_New, &
F_Lo, F_Hi, F_New
x_Lo = x_a
x_Hi = x_b
F_Lo = F(x_Lo)
F_Hi = F(x_Hi)
! Start of main Bisection loop
10 x_New = (x_Lo + x_Hi)/2.0
F_New = F(x_New)
! The point (X_New, F_New) replaces the
! interval boundary with the same sign of F,
! to ensure that the new interval still
! contains a sign change.
IF (F_Lo*F_New < 0.0) THEN
x_Hi = x_New
F_Hi = F_New
ELSE
x_Lo = x_New
F_Lo = F_New
END IF
! Repeat loop if have not fulfilled
! any convergence criterion
IF (Abs(F_New) > Tol .AND. &
(x_Hi - x_Lo) > Tol) GO TO 10
Bisect = X_New
RETURN
END
REAL FUNCTION NwtRph(x0, Tol)
! Finds one real root of the function F(x) using
! an initial guess at the root, x0.
REAL x0, Tol, &
x_New, x_Old, &
F_x, Fp_x
x_New = x0
F_x = F(x0)
! Main loop---iterate until one convergence
! criterion met.
10 x_Old = x_New
F_x = F(x_Old)
Fp_x = FPrime(x_Old)
x_New = x_Old - F_x/Fp_x
! Repeat loop if have not converged.
IF (Abs(F_x) > Tol .AND. &
Abs(x_New - x_Old) > Tol) GO TO 10
NwtRph = x_New
RETURN
END
REAL FUNCTION Secant(x_g1, x_g2, Tol)
! Finds one real root of the function F(x) using
! two initial guesses at the root, x_g1 and x_g2.
REAL x_g1, x_g2, Tol, &
x0, x1, x_New, &
F0, F1
x0 = x_g1
x1 = x_g2
F0 = F(x0)
F1 = F(x1)
! Ensure that x1 is the better guess to start
IF (Abs(F0) < Abs(F1)) THEN
CALL Swap(x0, x1)
CALL Swap(F0, F1)
END IF
! Start of main Secant loop
! Estimate new x as the point on the x-axis
! intersected by the line joining (x0, F0)
! and (x1, F1)
10 x_New = x1 - F1*(x1 - x0)/(F1 -F0)
! Update both estimates, assuming new x is
! better than the two previous guesses.
x0 = x1
F0 = F1
x1 = x_New
F1 = F(x1)
! Repeat loop if have not converged.
IF (Abs(F1) > Tol .AND. &
Abs(x1 - x0) > Tol) GO TO 10
Secant = x1
RETURN
END
DATA:
0.01
-2.0 2.0
.1
.1 .2
OUTPUT:
Program entered
This is Program >> P125 = Root finding demo
Program to demonstrate numerical root finding
F = x**3 - Exp(-x) for example
Enter tolerance for roots:
9.9999998E-03
Enter low limit, high limit:
-2.0000000 2.0000000
Bisection: 0.7734375
Enter estimated root:
0.1000000
Newton: 0.7728938
Enter 1st, 2nd root estimates:
0.1000000 0.2000000
Secant: 0.7726500
Fortran-90 STOP
Page builder: Charles Boivin
Last modified: 11/07/95