{**************************************************** * Program to demonstrate the complex domain * * Mueller's subroutine * * ------------------------------------------------- * * Reference: BASIC Scientific Subroutines, Vol. II * * By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * * * * Pascal Version By J-P Moreau, Paris. * * (www.jpmoreau.fr) * * ------------------------------------------------- * * Example: Find a complex root of zē + 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) = * * -1.43254447933788E-0005 -1.00003709578020E+0000 * * * * The associated Z value is (U,V) = * * -7.41927312872814E-0005 2.86519524138603E-0005 * * * * The number of steps was: 100 * ****************************************************} PROGRAM DEMO_ZMUELLER; Uses WinCrt; Var k,n : INTEGER; b1,b2,e,u,v,x,y,x0,y0 : DOUBLE; {************************************* Function subroutine } PROCEDURE Z(x,y:DOUBLE; VAR u,v:DOUBLE); Begin u := x*x-y*y+1; v := 2.0*x*y End; {************************************} {************************************************ * 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. * * --------------------------------------------- * * writeS: * * 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, * * ************************************************} PROCEDURE ZMueller; Label 100,200,300,fin; Var d,d1,d2,x1,x2,x3,y1,y2,y3 : DOUBLE; u1,u2,u3,v1,v2,v3,xl1,xl2 : DOUBLE; a,a1,a2,b,c1,c2,e1,e2,e3 : DOUBLE; Begin {Start calculations} k:=1; x3:=x0; y3:=y0; x1:=x3-b1; y1:=y3-b2; x2:=x3+b1; y2:=y3+b2; 100: 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} Z(x1,y1,u1,v1); Z(x2,y2,u2,v2); Z(x3,y3,u3,v3); {Calculate Mueller parameters} e1:=u1*(xl1*xl1-xl2*xl2)-2.0*v1*xl1*xl2-u2*(d1*d1-d2*d2); e1:=e1+2.0*v2*d1*d2+u3*(xl1+d1)-v3*(xl2+d2); e2:=2.0*xl1*xl2*u1+v1*(xl1*xl1-xl2*xl2)-2.0*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.0*(u3*d1*c1-u3*d2*c2-v3*d2*c1-v3*d1*c2); b2:=2.0*e1*e2-4.0*(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:=ARCTAN(b2/b1); a:=a/2.0; b:=SQRT(SQRT(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 200; a1:=e1-b1; a2:=e2-b2; goto 300; 200: a1:=e1+b1; a2:=e2+b2; 300: 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.0*xl1/a; xl2:=-d1*u3*a2+d2*v3*a2+a1*u3*d2+a1*v3*d1; xl2:=-2.0*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 goto fin; {Continue} k:=k+1; x0:=x; y0:=y; x1:=x2; y1:=y2; x2:=x3; y2:=y3; x3:=x; y3:=y; goto 100; fin: End; {of ZMueller} {main program} BEGIN clrscr; writeln; writeln(' Write the initial guesses and their bounds:'); writeln; write(' X0 = '); read(x0); write(' Bond on X0 = '); read(b1); writeln; write(' Y0 = '); read(y0); write(' Bond on Y0 = '); read(b2); writeln; write(' Convergence criterion: '); read(e); writeln; write(' Maximum number of iterations: '); read(n); ZMueller; {Call complex domain Mueller subroutine} writeln; writeln; writeln(' The calculated zero is (X,Y) ='); writeln(' ',x,' ',y); writeln; Z(x,y,u,v); writeln(' The associated Z value is (U,V) ='); writeln(' ',u,' ',v); writeln; writeln(' The number of steps was: ',k); writeln; ReadKey; DoneWinCrt END. {End of file zmueller.pas}