C============================================================== C A simple algorithm for computing all the zeros of a nonlinear C function of a variable in a given interval [a,b]. b>a. C C The algorithms tries to determine all the zeros of a nonlinear C function in a given interval by discretizing of this interval C and by find the intervals where the function changes its sign. C C Every subinterval is reduced by its halving for finding the C zero of the function under accuracy epsm. C C Neculai Andrei C April 16, 1975 C============================================================== subroutine func(x,f) real*8 x,f c f = x**2 - 4.d0 c f = x**2 - 8.d0*x - 9.d0 c f = 4.d0*x**3 - 3.d0*x - 1.d0 c f = x**3 + 4.d0*x**2 - 4.d0*x - 16.d0 c f = 1.d0/((x-0.3d0)**2+0.01) + 1.d0/((x-0.9d0)**2+0.04) -6.d0 c f = x**4 - 5.d0*x**2 + 4.d0 c f = x**6 - 14.d0*x**4 + 49.d0*x**2 - 36.d0 c f = x**7 + 7.d0*x**6 - 14.d0*x**5 - 98.d0*x**4 + c + 49.d0*x**3 +343.d0*x**2 - 36.d0*x - 252.d0 c f = x**3 - 2.d0*x**2 + 1.d0 c f = cos(x) - x**2 c f = x*dexp(x) - 7.d0 c f = x - cos(x) c f = cos(x) c f = sin(x) c f = cos(x**2) f = x**3 + 3.d0*x**2 - 18.d0*x return end c=========================================== Main program c ============== real*8 a, b, c, x real*8 fa,fb,fc,fx real*8 h real*8 epsm integer iprint integer i, ni, iint integer iter, itertot integer fgcnt, fgcntot epsm=0.000000001d0 iprint = 0 C Searching interval c ================== a = -7.d0 b = 7.d0 open(unit=4,file='zeros.out',status='unknown') itertot = 0 fgcntot = 0 write(4,1) 1 format(4x,'All zeros of a nonlinear function') write(4,2) a,b 2 format(4x,'in interval [',e20.13,','e20.13,' ]') write(4,3) 3 format(4x,56('='),/) C----------------------------------------------------- call func(a, fa) fgcnt = fgcnt + 1 C Take a discretization of the interval [a,b] C =========================================== ni = 200 h = (b-a)/float(ni) c Start iteratuions c ================= iint= 0 i = 1 400 continue iter = 0 fgcnt = 0 if(i .gt. ni) go to 998 30 b = a + h iter = iter+1 call func(b, fb) fgcnt = fgcnt + 1 if(fa*fb .gt. 0.d0) then a = b fa = fb i = i+1 if(i .gt. ni) go to 998 go to 30 else c write(4,104) c104 format(4x,'fa*fb le 0.d0: Change of sign------') iint = iint+1 if(fa .le. 0.d0) then 10 continue if(dabs(fa) .le. epsm .or. dabs(fb) .le. epsm) go to 999 iter=iter+1 c = a + (b-a)/2.d0 call func(c,fc) fgcnt = fgcnt + 1 if(fc .ge. 0.d0) then b = c fb = fc if(iprint .eq. 1) then write(4,99) iter 99 format(4x,'iter='i5) write(4,110) a,b,c 110 format(4x,' a='e20.13,4x,' b=',e20.13,4x,' c=',e20.13) write(4,111) fa,fb,fc 111 format(4x,'fa='e20.13,4x,'fb=',e20.13,4x,'fc=',e20.13) write(4,112) b-a 112 format(4x,'fa<0. b-a=',e20.13) end if go to 10 else a = c fa = fc if(iprint .eq. 1) then write(4,99) iter write(4,110) a,b,c write(4,111) fa,fb,fc write(4,112) b-a end if go to 10 end if end if if(fa. ge. 0.d0) then 20 continue if(dabs(fa) .le. epsm .or. dabs(fb) .le. epsm) go to 999 iter=iter+1 c = a + (b-a)/2.d0 call func(c,fc) fgcnt = fgcnt + 1 if(fc .ge. 0.d0) then a = c fa = fc if(iprint .eq. 1) then write(4,99) iter write(4,110) a,b,c write(4,111) fa,fb,fc write(4,112) b-a end if go to 20 else b = c fb = fc if(iprint .eq. 1) then write(4,99) iter write(4,110) a,b,c write(4,111) fa,fb,fc write(4,112) b-a end if go to 20 end if end if end if c end do c------------------------------------------------------------ 999 continue if(dabs(fa) .le. epsm .or. dabs(fb) .le. epsm) x=a call func(x, fx) fgcnt = fgcnt + 1 write(4,119) iint 119 format(2x,'Zero #',i2) write(4,120) x, fx 120 format(2x,'ZERO: ',e20.13,4x,'Function value: ',e20.13) write(4,121) iter 121 format(2x,'Number of iterations: ',i7) write(4,122) fgcnt 122 format(2x,'Number of function evaluations:',i7,/) a = x + h call func(a,fa) c fgcnt = fgcnt + 1 i = i+1 itertot = itertot + iter fgcntot = fgcntot + fgcnt go to 400 998 continue write(4,130) 130 format(2x,39('=')) write(4,131) itertot 131 FORMAT(2X,'TOTAL # OF iterations: ',i6) write(4,132) fgcntot 132 format(2x,'TOTAL # of function evaluations:',i6) stop end c======================================================================== c Last line