Hi everyone!

I want to solve f(x)=0 by using Bisection’s method. In this way, my function f(x) is any function such as f(x)=A.x^3+Bx^2+Cx+D. So I set an array (x,f(x)) then I use the subroutine having Bisectional method. My expectation is the solution having three roots. If the equation has 2 roots, I set one of the roots is 0.1 for my purpose. My interval is from 0 to 100.

The problem is all the equation I input the roots always appear 0.0 or no roots found.

Here is my code

```
program any_equation
implicit none
double precision A,B,C,D,E,dminroot,dmaxroot
double precision r,dr,r1,r2,r3,ff
integer i,nroots
double precision x(3000),f(3000),roots(10)
! input
print*, "Input A,B,C,D,E:"
read *, A,B,C,D,E
write(*,12) A,B,C,D,E
12 FORMAT('Your equation: ',F7.2,'X^3+',F7.2,'X^2+',F7.2,'X+',F7.2,'=',F7.2)
dr=0.1d0
! creating an array
open(20,file='array_x_fx',status='unknown')
do i=0,1000
r=dr*dfloat(i)
call find_f(A,B,C,D,r,ff)
x(i)=dr*dfloat(i)
f(i)=ff
write(20,*) x(i),f(i)
end do
close (20)
! find roots
call find_roots(roots,nroots,x,f,E)
!----
if (nroots ==0) then
write(*,*) 'no roots found'
stop
end if
!-----
if (nroots.ge.4) stop
!-----
if (nroots.eq.2) then
if(roots(1).lt.roots(2)) then
r2=roots(1)
r3=roots(2)
else
r2=roots(2)
r3=roots(1)
end if
r1=0.1d0
end if
!-----
if (nroots.eq.3) then
dmaxroot=roots(1)
dminroot=roots(1)
do i =1,nroots
if (dmaxroot.lt.roots(i)) dmaxroot = roots(i)
if (dminroot.gt.roots(i)) dminroot = roots(i)
enddo
do i =1,nroots
if((roots(i).lt.dmaxroot).and.(roots(i).gt.dminroot)) then
r2= roots(i)
end if
enddo
r1 =dminroot
r3 =dmaxroot
end if
write(*,131) r1,r2,r3
131 format('Solution are:',/,F7.2,/,F7.2,/,F7.2)
end program
! tinh f(x)
Subroutine find_f(A,B,C,D,r,f)
implicit none
double precision A,B,C,D,f,r
f=A*r**3.0d0+B*r**2.0d0+C*r+D
end subroutine timf
Subroutine find_roots(Roots,nroots,x,f,E)
!=====================================================================
! 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|<eps
! K - number of subintervals for [x1,x2]
! output ...
! Root(n) - roots of the equation f(x)=0 on [x1,x2]
! nroots - number of roots (nroots<=n)
! Comments:
! The program divide [x1,x2] into n subintervals
! Max number of iterations for every subinterval - 200
!======================================================================
implicit none
integer, parameter:: iter=200
integer nroots,i,j
double precision x1,x2,dx,e
double precision xmin,xmax,fmin,fmax,deriv3
double precision fa,fae,fb,fbe,fc,fce,froot,froote
double precision a,b,c,root,eps
double precision roots(10), y2(3000), x(3000), f(3000)
! initializayion
eps = 1.0e-6
x1=0.0d0
x2=100.0d0
dx = (x2-x1)/60. !the number of interval
nroots = 0
xmin=x(1)
xmax=x(1000)
fmin= deriv3(xmin,x,f,1000,1)
fmax= deriv3(xmax,x,f,1000,1)
call SPLINE(X,F,1000,fmin,fmax,y2)
! loop over subintervals
do j=1,60
a = x1 + real(j-1)*dx
b = a + dx
CALL SPLINT(x,f,y2,1000,a,fa)
fae=fa-e
CALL SPLINT(x,f,y2,1000,b,fb)
fbe=fb-e
!======================================================================
! check the closed domain condition f(a)*f(b)<0
if(fae*fbe>0) cycle
! Iterative refining the solution
do i=1,iter
c=(b+a)/2.0
! find c
CALL SPLINT(x,f,y2,1000,c,fc)
fce=fc-e
if(fae*fce.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.0d0
CALL SPLINT(x,f,y2,1000,root,froot)
froote=froot-e
if (abs(froote) < 1.0) then
nroots = nroots+1
Roots(nroots)=root
end if
end do
end subroutine find_roots
c------------------------------------------------------spline
SUBROUTINE SPLINE(x,y,n,yp1,ypn,y2)
INTEGER n,NMAX
REAL*8 yp1,ypn,x(n),y(n),y2(n)
PARAMETER (NMAX=20000)
INTEGER i,k
REAL*8 p,qn,sig,un,u(NMAX)
if (yp1.gt.0.99D40) then
y2(1)=0.0D0
u(1)=0.0D0
else
y2(1)=-0.50D0
u(1)=(3.0D0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
endif
do 11 i=2,n-1
sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
p=sig*y2(i-1)+2.0D0
y2(i)=(sig-1.0D0)/p
u(i)=(6.0D0*((y(i+1)-y(i))/(x(i+1)-x(i))-
& (y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
11 continue
if (ypn.gt.0.99D40) then
qn=0.0D0
un=0.0D0
else
qn=0.5D0
un=(3.0D0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
endif
y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.0D0)
do 12 k=n-1,1,-1
y2(k)=y2(k)*y2(k+1)+u(k)
12 continue
return
END
c------------------------------------------------------splint
SUBROUTINE SPLINT(xa,ya,y2a,n,x,y)
INTEGER n
REAL*8 x,y,xa(n),y2a(n),ya(n)
INTEGER k,khi,klo
REAL*8 a,b,h
klo=1
khi=n
1 if (khi-klo.gt.1) then
k=(khi+klo)/2
if(xa(k).gt.x)then
khi=k
else
klo=k
endif
goto 1
endif
h=xa(khi)-xa(klo)
if(h.eq.0.00D0) stop
a=(xa(khi)-x)/h
b=(x-xa(klo))/h
y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*
&(h**2)/6.0D0
return
END
c------------------------------------------------------deriv3
function deriv3(xx, xi, yi, ni, m)
!====================================================================
! Evaluate first- or second-order derivatives
! using three-point Lagrange interpolation
! written by: Alex Godunov (October 2009)
!--------------------------------------------------------------------
! input ...
! xx - the abscissa at which the interpolation is to be evaluated
! xi() - the arrays of data abscissas
! yi() - the arrays of data ordinates
! ni - size of the arrays xi() and yi()
! m - order of a derivative (1 or 2)
! output ...
! deriv3 - interpolated value
!============================================================================*/
implicit none
integer, parameter :: n=3
double precision deriv3, xx
integer ni, m
double precision xi(ni), yi(ni)
double precision x(n), f(n)
integer i, j, k, ix
! exit if too high-order derivative was needed,
if (m > 2) then
deriv3 = 0.0
return
end if
! if x is ouside the xi(1)-xi(ni) interval set deriv3=0.0
if (xx < xi(1) .or. xx > xi(ni)) then
deriv3 = 0.0
return
end if
! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i)
i = 1
j = ni
do while (j > i+1)
k = (i+j)/2
if (xx < xi(k)) then
j = k
else
i = k
end if
end do
! shift i that will correspond to n-th order of interpolation
! the search point will be in the middle in x_i, x_i+1, x_i+2 ...
i = i + 1 - n/2
! check boundaries: if i is ouside of the range [1, ... n] -> shift i
if (i < 1) i=1
if (i + n > ni) i=ni-n+1
! just wanted to use index i
ix = i
! initialization of f(n) and x(n)
do i=1,n
f(i) = yi(ix+i-1)
x(i) = xi(ix+i-1)
end do
! calculate the first-order derivative using Lagrange interpolation
if (m == 1) then
deriv3 =(2.0*xx - (x(2)+x(3)))*f(1)/((x(1)-x(2))*(x(1)-x(3)))
deriv3 = deriv3 + (2.0*xx - (x(1)+x(3)))*f(2)/((x(2)-
1 x(1))*(x(2)-x(3)))
deriv3 = deriv3 + (2.0*xx - (x(1)+x(2)))*f(3)/((x(3)-
2 x(1))*(x(3)-x(2)))
! calculate the second-order derivative using Lagrange interpolation
else
deriv3 = 2.0*f(1)/((x(1)-x(2))*(x(1)-x(3)))
deriv3 = deriv3 + 2.0*f(2)/((x(2)-x(1))*(x(2)-x(3)))
deriv3 = deriv3 + 2.0*f(3)/((x(3)-x(1))*(x(3)-x(2)))
end if
end function deriv3
!-------------------------------------------------------
Function gauss8(xx,fx,dx,a,b)
!==========================================================
! Integration of f(x) on [a,b]
! Method: Gauss 8 points
! written by: Alex Godunov (October 2009)
!----------------------------------------------------------
! IN:
! f - Function to integrate (supplied by a user)
! a - Lower limit of integration
! b - Upper limit of integration
! OUT:
! gauss8 - Result of integration
!==========================================================
implicit none
integer, parameter :: n=4
double precision gauss8, a, b, mi,fmi,mti,fmti
double precision ti(n), ci(n)
data ti/0.1834346424,0.5255324099,0.7966664774,0.9602898564/
data ci/0.3626837833,0.3137066458,0.2223810344,0.1012285362/
double precision r, m, c
integer i,dx
double precision xmin,xmax,fmin,fmax,deriv3
double precision xx(1000),fx(1000),y2(1000)
xmin=xx(1)
xmax=xx(dx)
fmin=deriv3(xmin,xx,fx,dx,1)
fmax=deriv3(xmax,xx,fx,dx,1)
call SPLINE(xx,fx,dx,fmin,fmax,y2)
r = 0.0;
m = (b-a)/2.0;
c = (b+a)/2.0;
do i = 1,n
mi= m*(-1.0)*ti(i) + c
call SPLINT(xx,fx,y2,dx,mi,fmi)
mti=m*ti(i) + c
call SPLINT(xx,fx,y2,dx,mti,fmti)
r = r + ci(i)*(fmi + fmti)
end do
gauss8 = r*m
return
end function gauss8
```