'***************************************************** '* Program to demonstrate the complex domain * '* Mueller's subroutine * '* ------------------------------------------------- * '* Reference: BASIC Scientific Subroutines, Vol. II * '* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * '* ------------------------------------------------- * '* Example: Find a complex root of z^2 + 1 = 0 * '* * '* SAMPLE RUN: * '* * '* Write the initial guesses and their bounds: * '* * '* X0 = 0 * '* Bond on X0 = 3 * '* Y0 = 0 * '* Bond on Y0 = 3 * '* * '* Convergence criterion: 0 * '* * '* Maximum number of iterations: 100 * '* * '* The calculated zero is (X,Y) = * '* -0.00001400 -1.00003402 * '* * '* The associated Z value is (U,V) = * '* -0.00006805 0.00002801 * '* * '* The number of steps was: 100 * '***************************************************** defint i-n defdbl a-h,o-z cls 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 gosub 2000 'Call complex domain Mueller subroutine print : print print " The calculated zero is (X,Y) = " print using " ##.######## ##.########";x;y print gosub 1000 print " The associated Z value is (U,V) = " print using " ##.######## ##.########";u;v print print " The number of steps was: "; k print end '************************************* 1000 ' Function subroutine u = x*x-y*y+1 v = 2#*x*y 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, * * '************************************************* 2000 'Start calculations k=1 x3=x0 : y3=y0 : x1=x3-b1 : y1=y3-b2 : x2=x3+b1 : y2=y3+b2 2100 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 x=x1 : y=y1 : gosub 1000 : u1=u : v1=v x=x2 : y=y2 : gosub 1000 : u2=u : v2=v x=x3 : y=y3 : gosub 1000 : 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 2200 a1=e1-b1 : a2=e2-b2 : goto 2300 2200 a1=e1+b1 : a2=e2+b2 2300 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 x=x3+xl1*(x3-x2)-xl2*(y3-y2) y=y3+xl2*(x3-x2)+xl1*(y3-y2) 'Test for convergence if ABS(x-x0)+ABS(y-y0)=n then return 'Continue k=k+1 : x0=x : y0=y : x1=x2 : y1=y2 : x2=x3 : y2=y3 : x3=x : y3=y goto 2100 'end of subroutine 2000 'End of file zMueller.bas