program main !==================================================================== ! Multiple roots of an equation f(x)=0 in [x1,x2] interval ! Method: Brute force + Closed Domain (Bisectional or False position) !==================================================================== implicit none integer, parameter:: n=1000 double precision f,x1,x2,eps double precision Roots(n) integer key, nroots, i external f key = 2 x1 = -10.0 x2 = 10.0 eps = 1.0e-7 write (*,*) ' Roots of f(x) by brute force method' if(key == 1) write (*,*) ' Subintervals - Bisectional method' if(key == 2) write (*,*) ' Subintervals - False position method' write (*,100) x1, x2 write (*,101) n write (*,102) eps call BForce(f,x1,x2,eps,Roots,key,n,nroots) ! print solutions if(nroots == 0) then write(*,*)' no roots found' stop end if write(*,*) ' number root f(root)' do i=1,nroots write(*,103) i, Roots(i), f(Roots(i)) end do 100 format(' interval [',f7.3,',',f7.3,']') 101 format(' subintervals= 'i12) 102 format(' tolerance = ',1pe12.3,/) 103 format(i5,3x,f11.6,2(1pe12.3)) end program main Subroutine BForce(f,x1,x2,eps,Roots,key,n,nroots) !============================================================ ! Multiple roots of equation f(x)=0 on [x1,x2] interval ! Method: Brute force with one of closed domain methods ! Close domain methods: bisectional or false position ! Alex G. January 2010 !------------------------------------------------------------ ! input ... ! f - function - evaluates f(x) for any x in [x1,x2] ! x1 - left endpoint of initial interval ! x2 - right endpoint of initial interval ! eps - desired uncertainity of the root as |b-a|0) cycle ! Iterative refining the solution do i=1,iter if(key == 1) then c=(b+a)/2.0 else c = b - f(b)*(b-a)/(f(b)-f(a)) end if if(f(a)*f(c).le.0.0) then b = c else a=c end if ! condition(s) to stop iterations) if(abs(b-a)<= eps) exit end do ! check if it is a root or singularity root = (b+a)/2.0 if (abs(f(root)) < 1.0) then nroots = nroots+1 Roots(nroots)=root end if end do end subroutine BForce !----------------------------------------------------- Function f(x) implicit none double precision f, x double precision, parameter::pi=3.141592654 f = exp(x)-x*sin(pi*x/2.0) end function f