*=======================================================================
* program main
* test program for root finding routine.
program main
implicit none
double precision func0, func1, func2, func3, func4, func5, func6,
& a0, b0, xend
integer nmbf
external findroot, func0, func1, func2, func3, func4,
& func5, func6
write (*,'(2x,a)') 'Bisection method'
write (*,'(2x,a8,2x,a16,2x,a16,2x,a6)')
& 'function', 'xend', 'error', 'nmbf'
a0 = 1
b0 = 2
call findroot_bisect (func0, a0, b0, xend, nmbf)
write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
& 'func0', xend, xend-sqrt(2.0d0), nmbf
a0 = 1
b0 = 2
call findroot_bisect (func1, a0, b0, xend, nmbf)
write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
& 'func1', xend, xend-sqrt(2.0d0), nmbf
a0 = 1
b0 = 2
call findroot_bisect (func2, a0, b0, xend, nmbf)
write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
& 'func2', xend, xend-sqrt(2.0d0), nmbf
a0 = 1
b0 = 2
call findroot_bisect (func3, a0, b0, xend, nmbf)
write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
& 'func3', xend, xend-sqrt(2.0d0), nmbf
a0 = 1
b0 = 2
call findroot_bisect (func4, a0, b0, xend, nmbf)
write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
& 'func4', xend, xend-sqrt(2.0d0), nmbf
a0 = 1
b0 = 2
call findroot_bisect (func5, a0, b0, xend, nmbf)
write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
& 'func5', xend, xend-sqrt(2.0d0), nmbf
a0 = -1
b0 = 2
call findroot_bisect (func6, a0, b0, xend, nmbf)
write (*,'(2x,a8,2x,1p,d16.7,2x,d16.7,2x,i6)')
& 'func6', xend, xend-0.0d0, nmbf
stop
end
*=======================================================================
* subroutine findroot_bisect
subroutine findroot_bisect (f, a0, b0, xend, nmbf)
implicit none
double precision f, a0, b0, xend
integer nmbf
external f
* ------------------------------------------------------------------
double precision signum, bisect, a, b, c, fa, fb, fc
external signum, bisect
a = a0
b = b0
fa = f(a)
fb = f(b)
nmbf = 2
if (fa.eq.0) then
c = a
fc = fa
else if (fb.eq.0) then
c = b
fc = fb
else
if (signum(fa).eq.signum(fb)) then
stop 'error in findroot_bisect -- no bracket for root'
endif
c = bisect(a,b)
fc = f(c)
nmbf = nmbf+1
endif
do while (c.ne.a.and.c.ne.b)
if (fc.eq.0) then
a = c
fa = fc
b = c
fb = fc
else if (signum(fc).eq.signum(fa)) then
a = c
fa = fc
c = bisect(a,b)
fc = f(c)
nmbf = nmbf+1
else if (signum(fc).eq.signum(fb)) then
b = c
fb = fc
c = bisect(a,b)
fc = f(c)
nmbf = nmbf+1
else
stop 'error in findroot_bisect -- f generates NaN?'
endif
enddo
xend = c
return
end
*=======================================================================
double precision function bisect(a,b)
implicit none
double precision a, b
* ------------------------------------------------------------------
* Modified bisection routine. If a and b are both non-zero, of the
* same sign, and not too far apart on a logarithmic scale, returns
* c=(a+b)/2. Other cases ... see the code. The intent is that
* bisect(a,b) divides the interval so that there are about equal
* number of floating point values in either half. Some preference
* towards zero.
double precision signum, eps, tiny, t0, t1
external signum
* let eps be IEEE machine precision
eps = 2.0d0**(-52)
* set tiny to IEEE smallest positive normal number
tiny = 2.0d0**(-1022)
if (signum(a)*signum(b).lt.0) then
bisect = 0
else if (signum(a)*signum(b).gt.0) then
t0 = min(abs(a),abs(b))
t1 = max(abs(a),abs(b))
if (t1.lt.4*t0.or.t1.lt.4*tiny) then
bisect = (a+b)/2
else
bisect = signum(a)*sqrt(t0)*sqrt(t1)
endif
else if (a.eq.0.or.b.eq.0) then
t0 = a+b
if (abs(t0).ge.1/sqrt(eps)) then
bisect = signum(t0)
else if (abs(t0).ge.sqrt(eps)) then
bisect = t0*sqrt(eps)
else if (abs(t0).ge.sqrt(tiny)) then
bisect = t0*abs(t0)
else if (abs(t0).ge.4*tiny) then
bisect = signum(t0)*tiny
else if (abs(t0).gt.eps*tiny) then
bisect = signum(t0)*eps*tiny
else
bisect = 0
endif
else
* must be NaN
bisect = (a+b)/2
endif
return
end
*=======================================================================
* double precision function signum(x)
* returns the sign of x; returns 0 if x.eq.0.
double precision function signum(x)
implicit none
double precision x
if (x.lt.0) then
signum = -1
else if (x.eq.0) then
signum = 0
else if (x.gt.0) then
signum = 1
else
* (must be that x is Not-a-Number)
signum = x
endif
return
end
*=======================================================================
* double precision function func0(x)
* test function for rootfinding routine.
* func0 is smooth and has a simple root at x=sqrt(2.0d0).
double precision function func0(x)
implicit none
double precision x
func0 = (x-sqrt(2.0d0))*sqrt(1+x**2)
return
end
*=======================================================================
* double precision function func1(x)
* test function for rootfinding routine.
* func1 is smooth and has a triple root at x=sqrt(2.0d0).
double precision function func1(x)
implicit none
double precision x
func1 = (x-sqrt(2.0d0))**3
return
end
*=======================================================================
* double precision function func2(x)
* test function for rootfinding routine.
* func2 is smooth and has a simple root at x=sqrt(2.0d0); the
* asymptotic behavior may pose problems for a Newton-like method.
double precision function func2(x)
implicit none
double precision x
func2 = atan(x-sqrt(2.0d0))
return
end
*=======================================================================
* double precision function func3(x)
* test function for rootfinding routine.
* func3 is smooth and has a simple root at x=sqrt(2.0d0); the
* asymptotic behavior may pose problems for a Newton-like method.
double precision function func3(x)
implicit none
double precision x
func3 = atan(256*(x-sqrt(2.0d0)))
return
end
*=======================================================================
* double precision function func4(x)
* test function for rootfinding routine.
* func4 is discontinuous with sign change at x=sqrt(2.0d0).
double precision function func4(x)
implicit none
double precision x
if (x.lt.0.or.x*x.lt.2) then
func4 = -1
else
func4 = 1
endif
return
end
*=======================================================================
* double precision function func5(x)
* test function for rootfinding routine.
* func5 has a vertical asymptote and sign change at x=sqrt(2.0d0).
double precision function func5(x)
implicit none
double precision x
double precision t
t = (x-sqrt(2.0d0))*sqrt(1+x**2)
if (t.ne.0) then
func5 = 1/t
else
func5 = 0
endif
return
end
*=======================================================================
* double precision function func6(x)
* test function for rootfinding routine.
* func6 is discontinous with sign change at 0.
double precision function func6(x)
implicit none
double precision x
if (x.lt.0) then
func6 = -1
else
func6 = 1
endif
return
end