'***************************************************** '* Program to demonstrate the complex domain * '* Allroot subroutine * '* ------------------------------------------------- * '* Reference: BASIC Scientific Subroutines, Vol. II * '* by F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * '* ------------------------------------------------- * '* Example: Find the two complex roots of z^2+1 = 0 * '* * '* SAMPLE RUN: * '* * '* Input the initial guesses and their bounds: * '* * '* X0 = 4 * '* Bond on X0 = 1 * '* Y0 = 4 * '* Bond on Y0 = 1 * '* * '* Convergence criterion: 1e-9 * '* Maximum number of iterations: 30 * '* How many roots are to be sought: 2 * '* Is the function defined in 1000 (1) * '* or is it a series (2): 1 * '* * '* The estimated roots are: * '* * '* X = 5E-22 * '* Y = 1 * '* * '* X = -2E-22 * '* Y = -.99999999 * '* * '* The last number of iterations was: 6 * '* * '***************************************************** defint i-n defdbl a-h,o-z DIM X(10),Y(10) cls PI=4#*ATN(1#) print print " Input the initial guesses and their bounds:" print input " X0 = ",x0 input " Bond on X0 = ",b1 print input " Y0 = ",y0 input " Bond on Y0 = ",b2 print input " Convergence criterion: ",e print input " Maximum number of iterations: ",n print input " How many roots are to be sought: ",n2 print print " Is the function defined in 1000 (1)" input " or is it a series (2) ", n3 gosub 2000 'Call complex domain Allroot subroutine print : print print " The estimated roots are:" for i=1 to n2 print print using " X = ##.########"; X(i) print using " Y = ##.########"; Y(i) next i print print " The last number of iterations was: "; k print end 'Rectangular to polar conversion 500 u=SQR(xx*xx+yy*yy) 'Guard against ambiguous vector if yy=0 then yy=1e-16 'Guard against divide by zero if xx=0 then xx=1e-16 v=ATN(yy/xx) 'Check quadrant and adjust if xx<0 then v=v+PI if v<0 then v=v+2#*PI return 'Polar to rectangular conversion 600 xx=u*COS(v) yy=u*SIN(v) return 'Polar power 700 u1=u^n : v1=n*v v1=v1-(2#*PI)*INT(v1/2#/PI) return 'Rectangular complex number power 800 'Rectangular to polar conversion gosub 500 'Polar power gosub 700 'Polar to rectangular conversion u=u1 : v=v1 : gosub 600 return '*** Function subroutine *** 1000 u=xx*xx-yy*yy+1 v=2#*xx*yy return '*********************************************** '* Complex series evaluation subroutine * '* ------------------------------------------- * '* The series coefficients are A(i), assumed * '* real. The order of the polynomial is m. The * '* routine uses repeated calls to the nth * '* power (Z^N) complex number subroutine. * '* INPUTS: x,y,m and A(i) * '* OUTPUTS: z1(real) and z2(imaginary) * '*********************************************** 1500 z1=A(0) : z2=0 'Store x and y a1=xx : a2=yy for n=1 to m 'Recall original x and y xx=a1 : yy=a2 'Call Z^N routine gosub 800 'Form partial sum z1=z1+A(n)*xx : z2=z2+A(n)*yy next n 'Restore x and y xx=a1 : yy=a2 return '**************************************************** '* General root determination subroutine * '* ------------------------------------------------ * '* This routine attempts to calculate the several * '* roots of a given series or function by repea- * '* tedly using the Zmueller subroutine and removing * '* the roots already found by division. * '* ------------------------------------------------ * '* INPUTS: X0, Y0 - The initial guess * '* b1, b2 - The bounds on this guess * '* e - The convergence criteria * '* n - The maximum number of iterations* '* for each root * '* n2 - the number of roots to seek * '* n3 - A flag (1) for a function F(Z) * '* (2) for a polynomial. * '* ------------------------------------------------ * '* OUTPUTS: The n2 roots found X(i), Y(i) * '* The last number of iterations, k. * '**************************************************** 2000 k=0 if n3=1 then goto 2100 if n3<>2 then return 'error in n3 2100 j1=0 'Save the initial guess x4=x0 : y4=y0 'If n3=2 get the series coefficients if n3=2 then gosub 2300 'Test for completion 2200 if j1=n2 then return 'Call Zmueller subroutine gosub 3000 j1=j1+1 X(j1)=xx : Y(j1)=yy : x0=x4 : y0=y4 'Try another pass goto 2200 '*** Coefficients subroutine *** 2300 m=5 : A(0)=0 : A(1)=24 : A(2)=-50 : A(3)=35 : A(4)=-10 : A(5)=1 return '************************************************* '* Mueller's method for complex roots subroutine * '* --------------------------------------------- * '* This routine uses the parabolic fitting tech- * '* nique associated with Mueller's method, but * '* does it in the complex domain. * '* --------------------------------------------- * '* INPUTS: * '* x0, y0 - The initial guess * '* b1, b2 - A bound on the error in this guess * '* e - The convergence criteria * '* n - The maximum number of iterations * '* OUTPUTS: * '* x,y - The value of the complex root found * '* k - The number of iterations performed, * * '************************************************* 3000 'Start calculations k=1 x3=x0 : y3=y0 : x1=x3-b1 : y1=y3-b2 : x2=x3+b1 : y2=y3+b2 3100 d=(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1) 'Avoid divide by zero if d=0 then d=1e-7 xl1=(x3-x2)*(x2-x1)+(y3-y2)*(y2-y1) xl1=xl1/d xl2=(x2-x1)*(y3-y2)+(x3-x2)*(y2-y1) xl2=xl2/d d1=(x3-x1)*(x2-x1)+(y3-y1)*(y2-y1) d1=d1/d d2=(x2-x1)*(y3-y1)+(x3-x1)*(y2-y1) d2=d2/d 'Get function values xx=x1 : yy=y1 : gosub 4000 : u1=u : v1=v xx=x2 : yy=y2 : gosub 4000 : u2=u : v2=v xx=x3 : yy=y3 : gosub 4000 : u3=u : v3=v 'Calculate Mueller parameters e1=u1*(xl1*xl1-xl2*xl2)-2#*v1*xl1*xl2-u2*(d1*d1-d2*d2) e1=e1+2#*v2*d1*d2+u3*(xl1+d1)-v3*(xl2+d2) e2=2#*xl1*xl2*u1+v1*(xl1*xl1-xl2*xl2)-2#*d1*d2*u2-v2*(d1*d1-d2*d2) e2=e2+u3*(xl2+d2)+v3*(xl1+d1) c1=xl1*xl1*u1-xl1*xl2*v1-d1*xl1*u2+xl1*d2*v2+u3*xl1 c1=c1-u1*xl2*xl2-v1*xl1*xl2+u2*xl2*d2+v2*d1*xl2-v3*xl2 c2=xl1*xl2*u1+xl1*xl1*v1-d2*xl1*u2-xl1*d1*v2+v3*xl1 c2=c2+u1*xl1*xl2-v1*xl2*xl2-u2*xl2*d1+v2*d2*xl2+u3*xl2 b1=e1*e1-e2*e2-4#*(u3*d1*c1-u3*d2*c2-v3*d2*c1-v3*d1*c2) b2=2#*e1*e2-4#*(u3*d2*c1+u3*d1*c2+v3*d1*c1-v3*d2*c2) 'Guard against divide by zero if b1=0 then b1=1e-7 a=ATN(b2/b1) : a=a/2# b=SQR(SQR(b1*b1+b2*b2)) b1=b*COS(a) : b2=b*SIN(a) a1=(e1+b1)*(e1+b1)+(e2+b2)*(e2+b2) a2=(e1-b1)*(e1-b1)+(e2-b2)*(e2-b2) if a1>a2 then goto 3200 a1=e1-b1 : a2=e2-b2 : goto 3300 3200 a1=e1+b1 : a2=e2+b2 3300 a=a1*a1+a2*a2 xl1=a1*d1*u3-a1*d2*v3+a2*u3*d2+a2*v3*d1 'Guard against divide by zero if a=0 then a=1e-7 xl1=-2#*xl1/a xl2=-d1*u3*a2+d2*v3*a2+a1*u3*d2+a1*v3*d1 xl2=-2#*xl2/a 'Calculate new estimate xx=x3+xl1*(x3-x2)-xl2*(y3-y2) yy=y3+xl2*(x3-x2)+xl1*(y3-y2) 'Test for convergence if ABS(xx-x0)+ABS(yy-y0)=n then return 'Continue k=k+1 : x0=xx : y0=yy : x1=x2 : y1=y2 : x2=x3 : y2=y3 : x3=xx : y3=yy goto 3100 'end of subroutine 3000 '*** Supervisor subroutine *** 4000 n5=n : u5=u1 : v5=v1 'Function or series ? if n3=1 then gosub 1000 if n3=2 then gosub 1500 if n3=1 then goto 4100 u=z1 : v=z2 'Restore parameters n=n5 : u1=u5 : v1=v5 4100 if j1=0 then return 'Divide by the j1 roots already found for j2=1 to j1 u5=u u=(xx-X(j2))*u+(yy-Y(j2))*v v=(xx-X(j2))*v-(yy-Y(j2))*u5 a4=(xx-X(j2))*(xx-X(j2))+(yy-Y(j2))*(yy-Y(j2)) 'Gard against divide by zero if a4=0 then a4=1e-7 v=v/a4 : u=u/a4 next j2 return 'End Sub 4000 'End of file Allroot.bas