'***************************************************************** '* Program to demonstrate Bisection & Quartile subroutines * '* ------------------------------------------------------------- * '* References: * '* * '* Bisection: BASIC Scientific Subroutines, Vol. II * '* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * '* * '* Quartile: After an algorithm provided By Namir Shammas. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '* ------------------------------------------------------------- * '* Example: Find a real root of f(x)=(x+1)^5 * '* * '* Sample run: * '* * '* What is the initial range (X0,X1): * '* X0 = -5 * '* X1 = 0 * '* Convergence criterion: 1e-6 * '* * '* * '* Bisection: The calculated zero is X = -.9999996423721313 * '* * '* The associated Y value is Y = 5.850012206076229D-33 * '* * '* The number of steps was: 23 * '* * '* * '* Quartile: The calculated zero is X = -1.000000089406967 * '* * '* The associated Y value is Y = -5.712902544996317D-36 * '* * '* The number of steps was: 12 * '* * '***************************************************************** 'Note: The quartile method is slightly faster and more accurate ' than the classical Bisection Method. '----------------------------------------------------------------- 'PROGRAM TQUART DEFDBL A-H, O-Z DEFINT I-N ' real*8 a,b,root, e,x,x0,x1 ' integer m CLS PRINT PRINT " What is the initial range (X0,X1):" INPUT " X0 = ", x0 INPUT " X1 = ", x1 INPUT " Convergence criterion: ", e PRINT a = x0: b = x1 GOSUB 1000 'call Bisect(x0, x1, e, x, m) IF m <> 0 THEN PRINT PRINT " Bisection: The calculated zero is X = "; x PRINT GOSUB 500 'calculate Y(x) PRINT " The associated Y value is Y = "; Y PRINT PRINT " The number of steps was: "; m END IF GOSUB 2000 'call Quartile(a, b, e, root, istep) IF m <> 0 THEN PRINT PRINT " Quartile: The calculated zero is X = "; root PRINT x = root: GOSUB 500 PRINT " The associated Y value is Y = "; Y PRINT PRINT " The number of steps was: "; istep END IF PRINT END '*********************************************** 500 'function Y(x) Y = (x + 1) ^ 5 RETURN '*********************************************** '*********************************************** '* Bisection method subroutine * '* ------------------------------------------- * '* This routine iteratively seeks the zero of * '* function Y(x) using the method of interval * '* halving until the interval is less than e * '* in width. It is assumed that the function * '* Y(x) is available from a function routine. * '* ------------------------------------------- * '* Input values: range (x0,x1), and e. * '* Output values: root x, Y(x) and number of * '* steps m. * '*********************************************** 1000 'subroutine Bisect(x0,x1,e,x,m) 'Label: 10 m = 0 x = x0: GOSUB 500: Y1 = Y x = x1: GOSUB 500 IF Y1 * Y > 0# THEN PRINT " Bisection: No root found in given interval." RETURN ELSE 10 x = x0: GOSUB 500: y0 = Y x = (x0 + x1) / 2: GOSUB 500: yy = Y m = m + 1 IF yy * y0 = 0# THEN RETURN IF yy * y0 < 0# THEN x1 = x IF yy * y0 > 0# THEN x0 = x IF ABS(x1 - x0) > e THEN GOTO 10 END IF RETURN '*********************************************** ' Quartile method subroutine * ' -------------------------------------------- * ' This routine iteratively seeks the zero of * ' function Y(x) using the method of interval * ' quarting until the interval is less than tol * ' in width. It is assumed that the function * ' Y(x) is available from a function routine. * ' -------------------------------------------- * ' Input values: range (a, b), and tol. * ' Output values: found root and number of used * ' steps. * '*********************************************** 2000 'subroutine Quartile(a, b, tol, root, istep) ' integer istep ' double xm, co tol = e co = .25# 'can be 0.3 istep = 0 x = a: GOSUB 500: Y1 = Y x = b: GOSUB 500 IF Y1 * Y > 0# THEN PRINT " Quartile: No root found in given interval." RETURN ELSE 20 x = a: GOSUB 500: Y1 = Y x = b: GOSUB 500 IF ABS(Y1) < ABS(Y) THEN xm = a + co * (b - a) ELSE xm = a + (1# - co) * (b - a) END IF x = xm: GOSUB 500: Y1 = Y x = a: GOSUB 500 IF Y1 * Y > 0# THEN a = xm ELSE b = xm END IF istep = istep + 1 IF ABS(a - b) > tol THEN GOTO 20 root = (a + b) / 2# END IF RETURN 'End of file tquart.bas