'******************************************************* '* Program to demonstrate the zero searching algorithm * '* --------------------------------------------------- * '* Reference: BASIC Scientific Subroutines, Vol. II * '* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * '* --------------------------------------------------- * * '* SAMPLE RUN: * '* What is the initial guess for x and y (x0,y0): * '* * '* X0 = 0 * '* Y0 = 0 * '* * '* What is the radius of the first search circle: 5 * '* * '* By what fraction this circle will be reduced on * '* each iteration: .5 * '* * '* How many evaluation points per quadrant: 4 * '* * '* Maximum number of iterations: 10 * '* * '* The approximate solution is: * '* Z = 0.000549 + 0.999919 I * '* * '* Number of iterations: 3 * '* * '******************************************************* DEFINT I-N DEFDBL A-H, O-Z CLS PRINT PRINT " What is the initial guess for x and y (x0,y0): " PRINT INPUT " X0 = ", x0 INPUT " Y0 = ", y0 PRINT INPUT " What is the radius of the first search circle: ", w PRINT INPUT " By what fraction this circle will be reduced on each iteration: ", e PRINT INPUT " How many evaluation points per quadrant: ", m PRINT INPUT " Maximum number of iterations: ", n PRINT DIM X(4 * m), Y(4 * m), U(4 * m), V(4 * m) GOSUB 2000 'Call ZCircle subroutine PRINT IF yy >= 0 THEN z$ = " + " IF yy < 0 THEN z$ = " - " PRINT " The approximate solution is:" PRINT USING " Z = ##.###### "; xx; : PRINT z$; PRINT USING " ##.######"; yy; : PRINT " I" PRINT PRINT " Number of iterations: "; k PRINT END '******************************************* 1000 ' Complex function(x,y) subroutine uu = xx * xx - yy * yy + 1# vv = 2# * xx * yy RETURN '******************************************* '************************************************* '* Complex root search subroutine * '* --------------------------------------------- * '* This routine searches for the complex roots * '* of an analytical function by encircling the * '* zero and estimating where it is. The circle * '* is subsequently tightened by a factor e, and * '* a new estimate made. * '* The inputs to the subroutine are: * '* w initial radius of search circle * '* x0,y0 the initial guess * '* e factor by which circle is reduced * '* n maximum number of iterations * '* m evaluation points per quadrant * '* The routine returns Z=X+IY (X,Y), the number * '* of iterations, k. * '************************************************* 2000 m = m * 4: k = 1 a = 2# * 3.1415926535# / m 2050 FOR i = 1 TO m X(i) = w * COS(a * (i - 1)) + x0 Y(i) = w * SIN(a * (i - 1)) + y0 NEXT i 'Determine the corresponding U(i) and V(i) FOR i = 1 TO m xx = X(i): yy = Y(i) GOSUB 1000 U(i) = uu V(i) = vv NEXT i 'Find the position at which uu changes sign 'in the counterclockwise direction i = 1: uu = U(i) 2100 GOSUB 3000 IF uu * U(i) < 0 THEN GOTO 2200 'Guard against infinite loop IF i = 1 THEN GOTO 2900 GOTO 2100 'Transition found 2200 m1 = i 'Search for the other transition, starting 'on the other side of the circle i = m1 + m / 2 IF i > m THEN i = i - m j = i uu = U(i) 'Flip directions alternately 2250 GOSUB 3000 IF uu * U(i) < 0 THEN GOTO 2300 IF uu * U(j) < 0 THEN GOTO 2310 'Guard against infinite loop IF i = m1 + xm2 / 2 THEN GOTO 2900 IF j = m1 + xm2 / 2 THEN GOTO 2900 GOTO 2250 'Transition found 2300 m3 = i GOTO 2320 'Transition found 2310 IF j = m THEN j = 0 m3 = j + 1 'm1 and m3 have been determined, now for xm2 and m4, 'now for the vv transitions 2320 i = m1 + m / 4 IF i > m THEN i = i - m j = i: vv = V(i) 2330 GOSUB 3000 IF vv * V(i) < 0 THEN GOTO 2340 IF vv * V(j) < 0 THEN GOTO 2350 'Guard against infinite loop IF i = m1 + m / 4 THEN GOTO 2900 IF j = m1 + m / 4 THEN GOTO 2900 GOTO 2330 'xm2 has been found 2340 xm2 = i GOTO 2400 2350 IF j = m THEN j = 0 xm2 = j + 1 'xm2 has been found, now for m4 2400 i = xm2 + m / 2 IF i > m THEN i = i - m j = i: vv = V(i) 2450 GOSUB 3000 IF uu * V(i) < 0 THEN GOTO 2460 IF vv * V(j) < 0 THEN GOTO 2470 'Guard against infinite loop IF i = xm2 + m / 2 THEN GOTO 2900 IF j = xm2 + m / 2 THEN GOTO 2900 GOTO 2450 2460 m4 = i GOTO 2500 2470 IF j = m THEN j = 0 m4 = j + 1 'All the intersections have been determined 'Interpolate to find the four (x,y) coordinates 2500 i = m1 GOSUB 4000 x1 = xx: y1 = yy: i = xm2 GOSUB 4000 x2 = xx: y2 = yy: i = m3 GOSUB 4000 x3 = xx: y3 = yy: i = m4 GOSUB 4000 x4 = xx: y4 = yy 'Calculate the intersection of the lines 'Guard against a divide by zero IF x1 <> x3 THEN GOTO 2510 xx = x1: yy = (y1 + y3) / 2# GOTO 2520 2510 m1 = (y3 - y1) / (x3 - x1) IF x2 <> x4 THEN GOTO 2520 xm2 = 1E+08 GOTO 2521 2520 xm2 = (y2 - y4) / (x2 - x4) 2521 b1 = y1 - m1 * x1 b2 = y2 - xm2 * x2 xx = -(b1 - b2) / (m1 - xm2) yy = (m1 * b2 + xm2 * b1) / (m1 + xm2) 'is another iteration in order ? IF k = n THEN RETURN x0 = xx: y0 = yy: k = k + 1: w = w * e GOTO 2050 2900 xx = 0: yy = 0: RETURN 3000 i = i + 1: j = j - 1 IF i > m THEN i = i - m IF j < 1 THEN j = j + m RETURN 4000 'interpolation routine j = i - 1 IF j < 1 THEN j = j + m 'Regula Falsi interpolation for the zero xx = (X(i) * U(j) + X(j) * U(i)) / (U(i) + U(j)) yy = (Y(i) * V(j) + Y(j) * V(i)) / (V(i) + V(j)) RETURN 'End of file Zcircle.bas