'*******************************************************
'* 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