DECLARE FUNCTION fct# (Num%, x#) DECLARE FUNCTION ipegasus% (Num%, x1#, x2#, f2#, iter%) '******************************************************************************* '* Test program for pegasus (6 functions defined in fonction.pas) * '* --------------------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* File tpegasus.lst contains: * '* * '* --------------------------------------------------------------------------- * '* Pegasus method for searching real roots of real valued, nonlinear functions * '* --------------------------------------------------------------------------- * '* * '* f(x) = 5 * x - exp (x) * '* Starting value x0 = 0 * '* Return code = 1 * '* Root = .2591711018190738 * '* Function value = 8.53267109746092D-17 * '* Iterations = 7 * '* * '* f(x) = (((((x-6)*x+15)*x-20)*x+15)*x-6)*x+1 * '* Starting value x0 = .5 * '* Return code = -1 * '* Root = 1.5 * '* Function value = .015625 * '* Iterations = 0 * '* * '* f(x) = sin (x) * '* Starting value x0 = 3 * '* Return code = 1 * '* Root = 3.141592653589793 * '* Function value = 1.22514845490862D-16 * '* Iterations = 6 * '* * '* f(x) = 1 + sin (x) * '* Starting value x0 = -2 * '* Return code = -1 * '* Root = -1.5 * '* Function value = 2.505013395945569D-03 * '* Iterations = 0 * '* * '* f(x) = exp(x) -(1.0 + x + x*x*0.5) * '* Starting value x0 = 2 * '* Return code = -1 * '* Root = 3 * '* Function value = 11.58553692318767 * '* Iterations = 0 * '* * '* f(x) = (x-1#)*(x-1#)*( sin(PI*x) - log(abs(2#*x/(x+1#))) * '* Starting value x0 = 2 * '* Return code = -1 * '* Root = 3 * '* Function value = -1.621860432432656 * '* Iterations = 0 * '* --------------------------------------------------------------------------- * '* * '* Ref.: "Numerical algorithms with C, By Gisela Engeln-Muellges and Frank * '* Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '******************************************************************************* 'Program Test_Pegasus DEFDBL A-H, O-Z DEFINT I-N DIM rc AS INTEGER DIM texte AS STRING OPEN "tpegasus.lst" FOR OUTPUT AS #2 'output file PRINT #2, " ---------------------------------------------------------------------------" PRINT #2, " Pegasus method for searching real roots of real valued, nonlinear functions" PRINT #2, " ---------------------------------------------------------------------------" FOR i = 1 TO 6 IF i = 1 THEN 'select suitable example x1 = 0# x2 = 1# NumFunc = 1 texte = "f(x) = 5 * x - exp (x)" ELSEIF i = 2 THEN x1 = .5# x2 = 1.5# NumFunc = 2 texte = "f(x) = (((((x-6)*x+15)*x-20)*x+15)*x-6)*x+1" ELSEIF i = 3 THEN x1 = 3# x2 = 3.5# NumFunc = 3 texte = "f(x) = sin (x)" ELSEIF i = 4 THEN x1 = -2# x2 = -1.5# NumFunc = 4 texte = "f(x) = 1 + sin (x)" ELSEIF i = 5 THEN x1 = 2# x2 = 3# NumFunc = 5 texte = "f(x) = exp(x) -(1.0 + x + x*x*0.5)" ELSEIF i = 6 THEN x1 = 2# x2 = 3# NumFunc = 6 texte = "f(x) = (x-1#)*(x-1#)*( sin(PI*x) - log(abs(2#*x/(x+1#)))" END IF PRINT #2, PRINT #2, " "; texte PRINT #2, " Starting value x0 = "; x1 rc = ipegasus(NumFunc, x1, x2, f, iter) PRINT #2, " Return code = "; rc PRINT #2, " Root = "; x2 PRINT #2, " Function value = "; f PRINT #2, " Iterations = "; iter NEXT i PRINT #2, " ---------------------------------------------------------------------------" CLOSE #2 CLS PRINT PRINT " Results in tpegasus.lst." PRINT END FUNCTION fct (Num, x) PI = 4# * ATN(1#) IF Num = 1 THEN fct = 5 * x - EXP(x) ELSEIF Num = 2 THEN fct = (((((x - 6) * x + 15) * x - 20) * x + 15) * x - 6) * x + 1 ELSEIF Num = 3 THEN fct = SIN(x) ELSEIF Num = 4 THEN fct = 1# + SIN(x) ELSEIF Num = 5 THEN fct = EXP(x) - (1# + x + x * x / 2#) ELSEIF Num = 6 THEN fct = ((x - 1#) ^ 2) * (SIN(PI * x) - LOG(2# * x / (x + 1#))) END IF END FUNCTION FUNCTION ipegasus (Num, x1, x2, f2, iter) '*====================================================================* '* * '* pegasus computes one zero of the continuous function fct, * '* provided that the two starting values x1 and x2 satisfy: * '* fct(x1) * fct(x2) <= 0 . * '* * '*====================================================================* '* * '* Applications: * '* ============ * '* Determine one root of the continuous function fct, if an * '* inclusion interval [x1, x2] is known for the root. * '* fct is defined in unit Fonction.pas. * '* * '*====================================================================* '* * '* Input parameters: * '* ================ * '* Num Integer; * '* Number of example (1 to 6) * '* x1,x2 Double; * '* Starting values with fct(x1) * fct(x2) <= 0. * '* * '* Output parameters: * '* ================= * '* x2 Double; * '* Computed approximation for the root of fct * '* * '* f2 Double; * '* Functional value at the approximate root, this * '* must be nearly zero. * '* iter integer; * '* Number of iterations that were performed. * '* * '* Return values: * '* ============= * '* = -1 No inclusion: fct(x2) * fct(x1) > 0 * '* = 0 Root has been found with ABS(f2) < FCTERR * '* = 1 Break-off with * '* ABS(xnew-xold) < ABSERR + xnew * RELERR, * '* check functional value * '* = 2 Iteration limit reached * '* * '*====================================================================* '* * '* Constants used: ABSERR, RELERR, MACH_EPS, EPSROOT, ITERMAX * '* ============== * '* * '*====================================================================* ITERMAX = 300 ' Maximal number of iteraions ABSERR = 0# ' Admissable absolute error XMACHEPS = .000000000000001#' Small real number RELERR = 4# * XMACHEPS ' Admissable relative error FCTERR = XMACHEPS * XMACHEPS' Maximal function errror DIM rc AS INTEGER rc = 2 iter = 0 ' Initialize iteration counter f1 = fct(Num, x1) ' Function values at x1, x2 f2 = fct(Num, x2) IF f1 * f2 > 0# THEN ipegasus = -1 ' No inclusion -> Error EXIT FUNCTION END IF IF f1 * f2 = 0# THEN ' One starting value is a root IF f1 = 0# THEN x2 = x1 f2 = 0# END IF ipegasus = 0 EXIT FUNCTION END IF WHILE iter <= ITERMAX ' Pegasus iteration iter = iter + 1 s12 = (f2 - f1) / (x2 - x1) ' Secant slope x3 = x2 - f2 / s12 ' new approximation f3 = fct(Num, x3) IF f2 * f3 <= 0# THEN ' new inclusion interval x1 = x2 f1 = f2 ELSE f1 = f1 * f2 / (f2 + f3) END IF x2 = x3 f2 = f3 IF ABS(f2) < FCTERR THEN ' Root found rc = 0 GOTO 10 END IF ' Break-off with small step size IF ABS(x2 - x1) <= ABS(x2) * RELERR + ABSERR THEN rc = 1 GOTO 10 END IF WEND 10 IF ABS(f1) < ABS(f2) THEN ' Choose approximate root with ' least magnitude function value x2 = x1 f2 = f1 END IF ipegasus = rc END FUNCTION