!******************************************************* !* Program to demonstrate the Newton root subroutine * !* in the complex domain * !* --------------------------------------------------- * !* 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: * !* * !* (Example: find a complex root of z^2 + 1 = 0) * !* * !* What is the initial guess: * !* * !* X0 = .2 * !* Y0 = 1.2 * !* * !* Convergence criterion: 1e-8 * !* Maximum number of iterations: 10 * !* * !* The root estimate is: * !* * !* X0 = 0.000000 * !* Y0 = 1.000000 * !* * !* The number of iterations performed was: 5 * !* * !******************************************************* PROGRAM DEMO_ZNEWTON integer k,n real*8 a,e,x0,y0,x,y print *,' ' print *,'What is the initial guess:' print *,' ' write(*,"(' X0 = ')",advance='no') read *, x0 write(*,"(' Y0 = ')",advance='no') read *, y0 print *,' ' write(*,"(' Convergence criterion: ')",advance='no') read *, e write(*,"(' Maximum number of iterations: ')",advance='no') read *, n call ZNewton(n,k,a,e,x0,y0,x,y) ! Call ZNewton subroutine if (a.ne.0.d0) then print *,' ' print *,'The root estimate is:' print *,' ' write(*,50) x write(*,60) y write(*,70) k end if stop 50 format(' X0 = ',f9.6) 60 format(' Y0 = ',f9.6) 70 format(/' The number of iterations performed was: ',I2//) end !********************************************************* ! Functions subroutine ! -------------------------------------------------------- Subroutine Eval(x,y,u,v,u1,v1,u2,v2) real*8 x,y,u,v,u1,v1,u2,v2 u=x*x-y*y+1.d0; v=2.d0*x*y u1=2.d0*x; u2=-2.d0*y v1=2.d0*y; v2=2.d0*x return end !********************************************************* !************************************************* !* Complex root seeking using Newton's method * !* --------------------------------------------- * !* This routine uses the complex domain form of * !* Newton's method for iteratively searching * !* for roots. The complex function and its first * !* partial derivatives must be available in the * !* form: F(Z) = U(X,Y) + I V(X,Y). The required * !* derivatives are DU/DX and DU/DY. * !* --------------------------------------------- * !* INPUTS; initial guess, x0 and y0, convergence * !* criteria, e, maximum number of iterations, n. * !* OUTPUTS; approximation to the root, x and y, * !* number of performed iterations, k. * !************************************************* Subroutine ZNewton(n,k,a,e,x0,y0,x,y) !Label: 100 parameter(TINY=1.d-12) !small number integer n,k real*8 e,x0,y0,x,y real*8 u,v,u1,v1,u2,v2 k=0 100 k=k+1 ! Get u,v and the derivatives u1,v1,u2,v2 call Eval(x0,y0,u,v,u1,v1,u2,v2) a=u1*u1+u2*u2 ! Guard against a=0 if (a < TINY) then print *,' ' print *,'ZERO DIVIDE ERROR - u1*u1+u2*u2 must be <> 0.' print *,' ' return end if x=x0+(v*u2-u*u1)/a y=y0-(v*u1+u*u2)/a ! Check for convergence in euclidean space if ((x0-x)*(x0-x)+(y0-y)*(y0-y)<=e*e) return if (k>=n) return x0=x ; y0=y goto 100 end ! ZNewton() ! End of file znewton.f90