!
! Program - P125.F90
!
PROGRAM Rooter
! methods.
! Function types!
IMPLICIT NONE
REAL :: Bisect, NwtRph, Secant
REAL :: Tol, X1, X2
!
INTERFACE
SUBROUTINE Swap(x1, x2)
IMPLICIT NONE
REAL ,INTENT(IN OUT) :: x1, x2
END SUBROUTINE Swap
END INTERFACE
!
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 PROGRAM Rooter
! SUPPORT ROUTINES
REAL FUNCTION F(x)
! The function of which a root (zero) is to be found.
F = x**3 - Exp(-x)
RETURN
END FUNCTION F
REAL FUNCTION FPrime(x)
! The derivative of F.
FPrime = 3*x*x + Exp(-x)
RETURN
END FUNCTION FPrime
SUBROUTINE Swap(x1, x2)
! Exchanges two real values x1 and x2.
IMPLICIT NONE
REAL ,INTENT(IN OUT) :: x1, x2
REAL :: Temp
Temp = x1
x1 = x2
x2 = Temp
RETURN
END SUBROUTINE Swap
! 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].
IMPLICIT NONE
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 FUNCTION Bisect
REAL FUNCTION NwtRph(x0, Tol)
! Finds one real root of the function F(x) using
! an initial guess at the root, x0.
IMPLICIT NONE
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 FUNCTION NwtRph
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.
IMPLICIT NONE
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 FUNCTION Secant
#include <stdio.h>
#include <math.h>
/* Program to demonstrate numerical root finding methods. Starting
from some initial guesses at a root (i.e. value of x such that
F(x) = 0), these algorithms iteratively improve the estimate of
the root until some convergence criterion is fulfilled.
Three techniques are presented. The Bisection method is easy to
understand and safe, but converges slowly. The Newton-Raphson
method converges quickly, but requires evaluation of the derivative
of F. The Secant method converges nearly as quickly as Newton-
Raphson, but does not require any derivative evaluations.
*/
#define MinTol 1.0E-10
#define TrueRoot 0.772883
#define MaxStep 20
double f(double x){
/* The function for which a root is to be found */
return(x*x*x - exp(-x));
}
double f_prime(double x){
/* The derivative of f */
return(3*x*x + exp(-x));
}
void swap(double *a, double *b){
/* Exchange two real variables */
double temp;
temp = *a;
*a = *b;
*b = temp;
}
double bisect(double x0, double x1, double tolerance){
/* Finds one real root of the function F(x) on the interval
[X_Lo, X_Hi], which must contain a sign change in F (i.e.
this will not work for roots of even multiplicity).
Proceeds by halving the interval [X_Lo, X_Hi]. Depending
on the sign of the function at the midpoint of the interval,
either the low or high limit is discarded (in order to preserve
the sign change in the interval).
The algorithm terminates when either the size of the interval
shrinks below Tolerance, or the magnitude of F at the midpoint
of the interval is less than Tolerance. */
double root;
int i;
for(i=1; i<=MaxStep; i++){
root=(x0 + x1) /2;
if(f(x0)*f(root) > 0.0)
x0=root;
else
x1=root;
if( (x1-root) < tolerance || fabs(f(root)) < tolerance) break;
}
return(root);
}
double newton_raphson(double x, double tolerance){
/* Finds one real root of the function F(x) using the Newton-
Raphson method and an initial guess at the root, x0.
While this method converges rapidly to the root, it does
require the evaluation of the derivative of F, F'(x).
The algorithm terminates when either the difference between
new and old estimates of the root falls below Tolerance, or
the magnitude of the function at the latest estimate is less
than Tolerance.
As a fail-safe, the algorithm is abandoned after MaxStep
(a global constant) iterations. */
double root;
int i;
for(i=1; i<=MaxStep; i++){
root = x - (f(x) / f_prime(x));
if( (root-x) < tolerance || fabs(f(root)) < tolerance) break;
x=root;
}
return(root);
}
double secant(double x0, double x1, double tolerance){
/* Finds one real root of the function F(x) using the Secant
method and two initial guesses at the root, x0 and x1.
This method has similar convergence properties to Newton-
Raphson, but does not require the evaluation of F'.
The algorithm terminates when either the difference between
new and old estimates of the root falls below Tolerance, or
the magnitude of the function at the latest estimate is less
than Tolerance.
As a fail-safe, the algorithm is abandoned after MaxStep
(a global constant) iterations.
If you understand the Newton-Raphson method, can you see how
the f(x0) * (x2-x1)/(f(x2) - f(x1)) term in this method replaces
the f(x)/f_prime(x) ? This secant method is actually approximating
the slope, to avoid evaluating the derivative (which may at times
be impossible to do.) */
double root;
int i;
for(i=1; i<=MaxStep; i++){
root = x1 - f(x0)*(x1-x0)/(f(x1)-f(x0));
if( fabs(x1-x0) < tolerance || fabs(f(root)) < tolerance) break;
x0 = x1;
x1= root;
}
return(root);
}
main(){
double x0, x1, tolerance, root;
printf("Root-finding Program\n\n\n");
printf(" This program finds the root of the function x^3 - e^(-x). It\n");
printf("will prompt you for two guesses, as well as the tolerance, and\n");
printf("then it will use the bisection method, the secant method, and \n");
printf("the Newton-Raphson method to find the root. The actual root \n");
printf("is approximately x = 0.772883 .\n\n\n");
printf("Enter first point value x=");
scanf("%lf",&x0);
printf("Enter last point value x=");
scanf("%lf",&x1);
printf("Enter Tolerance :");
scanf("%lf",&tolerance);
if(x0 > x1) swap(&x0, &x1);
if( fabs(f(x0)) <= tolerance || fabs(f(x1)) <= tolerance)
printf("One of the endpoints is a root.\n");
else if( f(x0)*f(x1) > 0 )
printf("There is no single root in that interval.\n");
else if(tolerance < MinTol)
printf("That is too small a tolerance.\n");
else if( (x1 - x0) < tolerance )
printf("The specified interval is smaller than the tolerance.\n");
else{
/* Approximate the root with bisection method */
root=bisect(x0, x1, tolerance);
/* If the root returned does not yield zero within tolerance,
it means the function could not find the root with the allowed
number of iterations. */
if(f(root) > tolerance){
printf("The bisection method could not find a root within\n");
printf("tolerance with %d iterations.\n",MaxStep); }
else{
printf("\nBisection method yields;\n");
printf("Root=%8.5lf True Root=%8.5lf Error=%8.5lf",
root, TrueRoot, (root-TrueRoot) ); }
/* Approximate the root with Newton-Raphson method */
root=newton_raphson(x0, tolerance);
if(f(root) > tolerance){
printf("\nThe Newton-Raphson method could not find a root within\n");
printf("tolerance with %d iterations.\n",MaxStep); }
else{
printf("\nNewton-Raphson method yields;\n");
printf("Root=%8.5lf True Root=%8.5lf Error=%8.5lf",
root, TrueRoot, (root-TrueRoot) );}
/* Approximate the root with the secant method */
root=secant(x0, x1, tolerance);
if(f(root) > tolerance){
printf("\nThe secant method could not find a root within\n");
printf("tolerance with %d iterations.\n",MaxStep); }
else{
printf("\nSecant method yields;\n");
printf("Root=%8.5lf True Root=%8.5lf Error=%8.5lf\n",
root, TrueRoot, (root-TrueRoot) ); }
} /* END the long if-else block */
} /* END main() */
PROGRAM RootDemo;
{ Program to demonstrate numerical root finding methods. Starting
from some initial guesses at a root (i.e. value of x such that
F(x) = 0), these algorithms iteratively improve the estimate of
the root until some convergence criterion is fulfilled.
Three techniques are presented. The Bisection method is easy to
understand and safe, but converges slowly. The Newton-Raphson
method converges quickly, but requires evaluation of the derivative
of F. The Secant method converges nearly as quickly as Newton-
Raphson, but does not require any derivative evaluations.
This program was compiled with Turbo Pascal V5.5, but should be
compilable with Turbo Pascal versions 3.0 and above. Deviations from
standard Pascal are also minimal, particularly in the coding of the
algorithms themselves. }
CONST
N_Alg = 3; { Number of algorithms implemented }
MinTol = 1.0E-10;
TrueRoot = 0.772883; { Found by experiment!! }
MaxStep = 20;
TYPE
{ Record containing information specific to an algorithm }
AlgRec = RECORD
Title: String[80]; { Identifying text string }
X0, X1: Real; { Initial guess(es) at root }
END;
VAR
Tolerance, { Convergence criterion
(allowed error in Root) }
Root: Real; { Algorithm's estimate of root }
Alg,
Step: Integer; { Number of steps to get Root }
Ch: Char;
AlgData: ARRAY [1..N_Alg] OF AlgRec;
FUNCTION F(x: Real): Real;
{ The function of which a root (zero) is to be found. }
BEGIN
F:= x*x*x - exp(-x);
END; { FUNCTION F }
FUNCTION F_Prime(x: Real): Real;
{ The derivative of F. }
BEGIN
F_Prime:= 3*x*x + exp(-x);
END; { FUNCTION F_Prime }
PROCEDURE Swap(VAR x1, x2: Real);
{ Exchanges two real values x1 and x2. }
VAR
Temp: Real;
BEGIN
Temp:= x1;
x1:= x2;
x2:= Temp;
END; { PROCEDURE Swap }
FUNCTION Bisect(
X_Lo, X_Hi,
Tolerance: Real;
VAR Step: Integer): Real;
{ Finds one real root of the function F(x) on the interval
[X_Lo, X_Hi], which must contain a sign change in F (i.e.
this will not work for roots of even multiplicity).
Proceeds by halving the interval [X_Lo, X_Hi]. Depending
on the sign of the function at the midpoint of the interval,
either the low or high limit is discarded (in order to preserve
the sign change in the interval).
The algorithm terminates when either the size of the interval
shrinks below Tolerance, or the magnitude of F at the midpoint
of the interval is less than Tolerance. }
VAR
X_New,
F_Lo, F_Hi,
F_New: Real;
BEGIN
{ Ensure proper ordering of initial interval. }
IF X_Lo > X_Hi THEN
Swap(X_Lo, X_Hi);
Tolerance:= Abs(Tolerance);
IF Tolerance < MinTol THEN
Tolerance:= MinTol;
Step:= 0;
F_Lo:= F(X_Lo);
F_Hi:= F(X_Hi);
{ First check whether one of the initial endpoints is a root. }
IF Abs(F_Lo) < Tolerance THEN
BEGIN
Bisect:= X_Lo;
Exit;
END
ELSE IF Abs(F_Hi) < Tolerance THEN
BEGIN
Bisect:= X_Hi;
Exit;
END
{ Complain if interval does not contain a change of sign of F. }
ELSE IF NOT (F_Lo*F_Hi < 0) THEN
BEGIN
WriteLn('***WARNING (Bisect):');
WriteLn(' No sign change of F on interval [',
X_Lo:7:3, ',', X_Hi:7:3, '].');
Step:= -1;
Exit;
END;
REPEAT
X_New:= (X_Lo + X_Hi)/2;
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 THEN
BEGIN
X_Hi:= X_New;
F_Hi:= F_New;
END
ELSE
BEGIN
X_Lo:= X_New;
F_Lo:= F_New;
END;
Step:= Succ(Step);
UNTIL (Abs(F_New) < Tolerance) OR
((X_Hi - X_Lo) < Tolerance);
Bisect:= X_New;
END; { FUNCTION Bisect }
FUNCTION Newton_Raphson(
x0, Tolerance: Real;
VAR Step: Integer): Real;
{ Finds one real root of the function F(x) using the Newton-
Raphson method and an initial guess at the root, x0.
While this method converges rapidly to the root, it does
require the evaluation of the derivative of F, F'(x).
The algorithm terminates when either the difference between
new and old estimates of the root falls below Tolerance, or
the magnitude of the function at the latest estimate is less
than Tolerance.
As a fail-safe, the algorithm is abandoned after MaxStep
(a global constant) iterations. }
VAR
X_New,
X_Old,
F_x,
Fp_x: Real;
BEGIN
Tolerance:= Abs(Tolerance);
IF Tolerance < MinTol THEN
Tolerance:= MinTol;
Step:= 0;
X_New:= x0;
X_Old:= X_New + 2*Tolerance;
F_x:= F(x0);
WHILE (Step < MaxStep) AND
(Abs(F_x) > Tolerance) AND
(Abs(X_New - X_Old) > Tolerance) DO
BEGIN
X_Old:= X_New;
F_x:= F(X_Old);
Fp_x:= F_Prime(X_Old);
IF Abs(Fp_x) < MinTol THEN
BEGIN
WriteLn('***WARNING (Newton_Raphson):');
WriteLn(' |F''(', X_Old:7:3, ')| --> 0.');
Step:= -1;
Exit;
END;
X_New:= X_Old - F_x/Fp_x;
Step:= Succ(Step);
END; { WHILE }
IF Step >= MaxStep THEN
BEGIN
WriteLn('***WARNING (Newton_Raphson):');
WriteLn(' Maximum of ', MaxStep:1, ' iterations exceeded.');
Step:= -1;
END
ELSE
Newton_Raphson:= X_New;
END; { Newton_Raphson }
FUNCTION Secant(
x0, x1, Tolerance: Real;
VAR Step: Integer): Real;
{ Finds one real root of the function F(x) using the Secant
method and two initial guesses at the root, x0 and x1.
This method has similar convergence properties to Newton-
Raphson, but does not require the evaluation of F'.
The algorithm terminates when either the difference between
new and old estimates of the root falls below Tolerance, or
the magnitude of the function at the latest estimate is less
than Tolerance.
As a fail-safe, the algorithm is abandoned after MaxStep
(a global constant) iterations. }
VAR
x_New,
F0, F1: Real;
BEGIN
Tolerance:= Abs(Tolerance);
IF Tolerance < MinTol THEN
Tolerance:= MinTol;
{ Complain if initial guesses are the same (within Tolerance). }
IF Abs(x0 - x1) < Tolerance THEN
BEGIN
WriteLn('***WARNING (Secant):');
WriteLn(' Initial guesses are not distinct.');
Step:= -1;
Exit;
END;
F0:= F(x0);
F1:= F(x1);
{ Ensure that x0 is the better guess to start. }
IF Abs(F0) < Abs(F1) THEN
BEGIN
Swap(x0, x1);
Swap(F0, F1);
END;
{ Check if x1 is a root }
IF Abs(F1) < Tolerance THEN
BEGIN
Secant:= x1;
Exit;
END;
Step:= 0;
WHILE (Step < MaxStep) AND
(Abs(F1) > Tolerance) AND
(Abs(x1 - x0) > Tolerance) DO
BEGIN
F0:= F1 - F0;
IF Abs(F0) < MinTol THEN
BEGIN
WriteLn('***WARNING (Secant):');
WriteLn(' |F(', x1:7:3, ') - F(', x0:7:3,')| --> 0.');
Step:= -1;
Exit;
END;
{ Estimate new x as the point on the x-axis intersected by
the line joining (x0, F0) and (x1, F1) }
x_New:= x1 - F1*(x1 - x0)/F0;
{ Update both estimates, assuming new x is better than the
two previous guesses. }
x0:= x1;
F0:= F1;
x1:= x_New;
F1:= F(x1);
Step:= Succ(Step);
END; { WHILE }
IF Step >= MaxStep THEN
BEGIN
WriteLn('***WARNING (Secant):');
WriteLn(' Maximum of ', MaxStep:1, ' iterations exceeded.');
Step:= -1;
END
ELSE
Secant:= x1;
END; { Secant }
BEGIN { Main }
AlgData[1].Title:= 'Bisection';
AlgData[2].Title:= 'Newton-Raphson';
AlgData[3].Title:= 'Secant';
REPEAT
REPEAT
Write('Enter tolerance for roots: ');
ReadLn(Tolerance);
UNTIL Tolerance > 0;
FOR Alg:= 1 TO N_Alg DO
WITH AlgData[Alg] DO
BEGIN
WriteLn;
WriteLn(Title+' Method');
CASE Alg OF
1: BEGIN
Write('Enter low limit, high limit: ');
ReadLn(X0, X1);
END;
2: BEGIN
Write('Enter estimated root: ');
ReadLn(X0);
END;
3: BEGIN
Write('Enter 1st, 2nd root estimates: ');
ReadLn(X0, X1);
END;
END; { CASE Alg }
END; { WITH AlgData[Alg] }
WriteLn; WriteLn;
WriteLn('Algorithm':15, 'Steps':8, 'Root':12, 'F(Root)':12,
'True':12, 'Error':12);
WriteLn;
FOR Alg:= 1 TO N_Alg DO
WITH AlgData[Alg] DO
BEGIN
CASE Alg OF
1: Root:= Bisect(X0, X1, Tolerance, Step);
2: Root:= Newton_Raphson(X0, Tolerance, Step);
3: Root:= Secant(X0, X1, Tolerance, Step);
END; { CASE Alg }
IF Step >= 0 THEN
WriteLn(Title:15, Step:8, Root:12:8, F(Root):12:8,
TrueRoot:12:8, (Root - TrueRoot):12:8);
END; { WITH AlgData[Alg] }
WriteLn; WriteLn;
Write('Another run? (y/[N]): ');
ReadLn(Ch);
UNTIL NOT ((Ch = 'Y') OR (Ch = 'y'));
END.
Last modified: 08/07/97