!******************************************************* !* Program to demonstrate the zero searching algorithm * !* --------------------------------------------------- * !* Reference; BASIC Scientific Subroutines, Vol. II * !* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * !* * !* F90 version by J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !* --------------------------------------------------- * !* 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.000576 + 0.999857 I * !* * !* Number of iterations: 10 * !* * !******************************************************* PROGRAM ZCIRCLE1 real*8 e,w,x0,y0,xx,yy integer k,m,n character*4 z print *,' ' print *,'What is the initial guess for x and y (x0,y0):' print *,' ' write(*,"(' X0 = ')",advance='no') read *, x0 write(*,"(' Y0 = ')",advance='no') read *, y0 print *,' ' write(*,"(' What is the radius of the first search circle: ')",advance='no') read *, w write(*,"(' By what fraction this circle will be reduced on each iteration: ')", & advance='no'); read *, e write(*,"(' How many evaluation points per quadrant: ')",advance='no') read *, m write(*,"(' Maximum number of iterations: ')",advance='no') read *, n print *,' ' call Zcircle(k,m,n,x0,y0,w,e,xx,yy) ! Call ZCircle subroutine if (yy>=0) z=' + ' if (yy< 0) z=' - ' print *,'The approximate solution is:' write (*,50) xx,z,yy write (*,60) k stop 50 format(' Z = ',F9.6,a3,F9.6,' I'/) 60 format(' Number of iterations: ',i2//) end !*************************************************** ! Complex function(x,y) subroutine ! -------------------------------------------------- Subroutine Eval(x,y,u,v) real*8 x,y,u,v u=x*x-y*y+1.0 v=2.0*x*y return end !*************************************************** ! Local functions used by Zcircle() Subroutine Set_to_zero(x,y) real*8 x,y x=0.0 y=0.0 return end Subroutine Auxilliary(m,i,j) integer m, i,j i=i+1; j=j-1 if (i>m) i=i-m if (j<1) j=j+m return end Subroutine Interpolation(m,i,j,X,Y,U,V,xx,yy) integer m,i,j real*8 xx,yy real*8 X(80),Y(80),U(80),V(80) j=i-1 if (j<1) 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 ! End local functions !************************************************* !* 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), and the * !* number of iterations, k. * !************************************************* Subroutine Zcircle(k,m,n,x0,y0,w,e,xx,yy) !Labels: 50,100,200,250,300,310,320,330,340,350,400,450 ! 460,470,500,510,520,521 integer k,m,n,i,j,m3,m4 real*8 x0,y0,w,e,xx,yy real*8 X(80),Y(80),U(80),V(80) real*8 a,b1,b2,PI,x1,x2,xm1,xm2 real*8 x3,y3,x4,y1,y2,y4,uu,vv m=m*4; k=1 PI = 4.d0*datan(1.d0) a=2.d0*PI/m 50 do i=1,m X(i)=w*dcos(a*(i-1))+x0 Y(i)=w*dsin(a*(i-1))+y0 end do ! Determine the corresponding U[i] and V[i] do i=1,m xx=X(i); yy=Y(i) call Eval(xx,yy,uu,vv) U(i)=uu V(i)=vv end do ! Find the position at which uu changes sign ! in the counterclockwise direction i=1; uu=U(i) 100 call Auxilliary(m,i,j) if (uu*U(i)<0) goto 200 ! Guard against infinite loop if (i.eq.1) call Set_to_zero(xx,yy) goto 100 ! Transition found 200 xm1=i ! Search for the other transition, starting ! on the other side of the circle i=xm1+m/2 if (i>m) i=i-m j=i uu=U(i) ! Flip directions alternately 250 call Auxilliary(m,i,j) if (uu*U(i)<0) goto 300 if (uu*U(j)<0) goto 310 ! Guard against infinite loop if (i.eq.xm1+(xm2/2)) call Set_to_zero(xx,yy) if (j.eq.xm1+(xm2/2)) call Set_to_zero(xx,yy) goto 250 ! Transition found 300 m3=i goto 320 ! Transition found 310 if (j.eq.m) j=0 m3=j+1 ! xm1 and m3 have been determined, now for xm2 and m4 ! now for the vv transitions 320 i=xm1+m/4 if (i>m) i=i-m j=i; vv=V(i) 330 call Auxilliary(m,i,j) if (vv*V(i)<0) goto 340 if (vv*V(j)<0) goto 350 ! Guard against infinite loop if (i.eq.xm1+m/4) call Set_to_zero(xx,yy) if (j.eq.xm1+m/4) call Set_to_zero(xx,yy) goto 330 ! xm2 has been found 340 xm2=i goto 400 350 if (j.eq.m) j=0 xm2=j+1 ! xm2 has been found, now for m4 400 i=xm2+m/2 if (i>m) i=i-m j=i; vv=V(i) 450 call Auxilliary(m,i,j) if (uu*V(i)<0) goto 460 if (vv*V(j)<0) goto 470 ! Guard against infinite loop if (i.eq.xm2+m/2) call Set_to_zero(xx,yy) if (j.eq.xm2+m/2) call Set_to_zero(xx,yy) goto 450 460 m4=i goto 500 470 if (j.eq.m) j=0; m4=j+1; ! All the intersections have been determine ! Interpolate to find the four (x,y) coordinates 500 i=xm1 call Interpolation(m,i,j,X,Y,U,V,xx,yy) x1=xx ; y1=yy ; i=xm2 call Interpolation(m,i,j,X,Y,U,V,xx,yy) x2=xx ; y2=yy ; i=m3 call Interpolation(m,i,j,X,Y,U,V,xx,yy) x3=xx ; y3=yy ; i=m4 call Interpolation(m,i,j,X,Y,U,V,xx,yy) x4=xx ; y4=yy ! Calculate the intersection of the lines ! Guard against a divide by zero if (x1.ne.x3) goto 510; xx=x1; yy=(y1+y3)/2.d0 goto 520 510 xm1=(y3-y1)/(x3-x1); if (x2.ne.x4) goto 520 xm2=1e8; goto 521 520 xm2=(y2-y4)/(x2-x4) 521 b1=y1-xm1*x1 b2=y2-xm2*x2 xx=-(b1-b2)/(xm1-xm2) yy=(xm1*b2+xm2*b1)/(xm1+xm2) ! is another iteration in order ? if (k.eq.n) return x0=xx ; y0=yy ; k=k+1 ; w=w*e goto 50 end ! Zcircle ! End of file zcircle.f90