! ! 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