Fortran Why Isn't My Fortran Bisection Code Finding the Root?

[Galizius]

Hello, I was told to find the root of an integrand with using the bisection method. I am kind of new to such a stuff so I am kind of confused why I am not getting the expected result.

Here is my bisection code:

tl=3.0
tup=4.0
tempred=0.0
s=-0.234726128709272990292
  do while (abs(total_calka) > eps)
    tempred=(tl+tup)/2.0
    beta=4.0/tempred
    if (total_calka * s < 0) then
      tup=tempred
    else
      tl=tempred
    endif
    call integration2
  end do

And also there is necessary to explain what other things mean. Eps is the approximation; tl and tup are values of an interval where the root is placed in. S is the value of total_calka for the "tl" temperature value. The total_calka statement is a value which I am getting inside subroutine integration2. And here is the code for that.

module int2
use variables;
use fx; use int1
implicit none
contains
subroutine integration2
open(1,file='calka.dat', access='append')
total_calka=0.0
    division=0.5
    tot_old=10000.0
    blad=10000.0
    do while (blad>error)
    irange=int(intrange/division)
    total_calka=-(0.5**3)/3.0
    do j=1,irange
    r1=0.5+(j-1)*division
    r2=r1+division
    calka=0.0
    do i=1,ngauss
      integrand(i)=inte(x(i),beta,r2,r1)
     calka=calka+integrand(i)*w(i)
    enddo
    total_calka=total_calka+calka*(r2-r1)
    enddo
        write(*,*) irange,division,tot_old,total_calka
    division=division/2.0
    blad=abs(tot_old-total_calka)
    tot_old=total_calka
       write(*,*) blad
    enddo
      write(1,*)tempred,-2*pi*total_calka
   close(1)
   end subroutine integration2
end module int2

There is also a function inte()

module fx
  implicit none
  contains
     real(10) function inte(y,beta,r2,r1)
     real(kind=10)::r,beta,r6,r2,r1,u,y
     r=(r2-r1)*y+r1
     r6=(1.0/r)**6
     u=beta*r6*(r6-1.0d0)
     if (u>100.d0) then
     inte=-1.0d0
     else
     inte=exp(-u)-1.d0
     endif
     inte=r*r*inte
     end function inte
end module fx

I don't know why it's not getting me a value of the root, while when I was doing a simple loop for temperature interval 0.1<T<5.0 I was getting the values of every total_calka for each temperature. I would be grateful for any hints that are helpful.

 

[Mark44]

Hello, I was told to find the root of an integrand with using the bisection method. I am kind of new to such a stuff so I am kind of confused why I am not getting the expected result.

Here is my bisection code:

tl=3.0
tup=4.0
tempred=0.0
s=-0.234726128709272990292
  do while (abs(total_calka) > eps)
    tempred=(tl+tup)/2.0
    beta=4.0/tempred
    if (total_calka * s < 0) then
      tup=tempred
    else
      tl=tempred
    endif
    call integration2
  end do

And also there is necessary to explain what other things mean. Eps is the approximation; tl and tup are values of an interval where the root is placed in. S is the value of total_calka for the "tl" temperature value. The total_calka statement is a value which I am getting inside subroutine integration2. And here is the code for that.

module int2
use variables;
use fx; use int1
implicit none
contains
subroutine integration2
open(1,file='calka.dat', access='append')
total_calka=0.0
    division=0.5
    tot_old=10000.0
    blad=10000.0
    do while (blad>error)
    irange=int(intrange/division)
    total_calka=-(0.5**3)/3.0
    do j=1,irange
    r1=0.5+(j-1)*division
    r2=r1+division
    calka=0.0
    do i=1,ngauss
      integrand(i)=inte(x(i),beta,r2,r1)
     calka=calka+integrand(i)*w(i)
    enddo
    total_calka=total_calka+calka*(r2-r1)
    enddo
        write(*,*) irange,division,tot_old,total_calka
    division=division/2.0
    blad=abs(tot_old-total_calka)
    tot_old=total_calka
       write(*,*) blad
    enddo
      write(1,*)tempred,-2*pi*total_calka
   close(1)
   end subroutine integration2
end module int2

There is also a function inte()

module fx
  implicit none
  contains
     real(10) function inte(y,beta,r2,r1)
     real(kind=10)::r,beta,r6,r2,r1,u,y
     r=(r2-r1)*y+r1
     r6=(1.0/r)**6
     u=beta*r6*(r6-1.0d0)
     if (u>100.d0) then
     inte=-1.0d0
     else
     inte=exp(-u)-1.d0
     endif
     inte=r*r*inte
     end function inte
end module fx

I don't know why it's not getting me a value of the root, while when I was doing a simple loop for temperature interval 0.1<T<5.0 I was getting the values of every total_calka for each temperature. I would be grateful for any hints that are helpful.

What do you mean by "root of an integrand"? Bisection is typically used to find a solution to the equation f(x) = 0 in some interval [a, b], and for which f(a) and f(b) have different signs. If f is continuous on the interval, there has to be a root somewhere in the interval. By bisecting the interval you're looking for smaller and smaller subintervals on which the function still has different signs.

I took a look at your code, but found that it was too hard to follow easily. Some of the things that make your code difficult to comprehend are the following:
1. Your integration2 subroutine has no parameters.
2. Your variable names don't make much sense to me (e.g., r, r1, r2, u, calka, blad, etc.).
3. There is not a single comment anywhere in your code.
4. One of your modules uses several other modules, but it's not clear to me whether the snippets you showed are from these other modules (variables, fx, int1)

 

[Galizius]

By a root of an integrand i mean the place where the total_calka = 0.

There is the variable module. I just forgotten to add it, because i thought it would be unnecessary

 module variables
implicit none
real(10) :: r, u, r6, tempred, f, r2, r1, calka,beta
real :: start, finish
integer:: i,k
integer, parameter :: Ngauss = 8
real(10),dimension(ngauss),parameter::xx=(/-0.9602898565d0,&
-0.7966664774d0,-0.5255324099d0,-0.1834346425d0,&
0.1834346425d0,0.5255324099d0,0.7966664774d0,0.9602898565d0/)
real(10),Dimension(ngauss),parameter::ww=(/0.1012285363d0,&
0.2223810345d0,0.3137066459d0,0.3626837834d0,&
0.3626837834d0,0.3137066459d0,0.2223810345d0,0.1012285363d0/)
real(10),dimension(ngauss)::x,w,integrand
integer, parameter :: out_unit=1000
integer, parameter :: out_unit1=1001
integer, parameter :: out_unit2=1002
real(10), parameter :: error=0.000001
real(10):: total_calka, division,tot_old,blad
real(10),parameter:: intrange=7.0
integer::j,irange
real, parameter :: pi = 4.0d0*datan(1.0d0)
end module variables

The whole code is about to count the value of second virial coefficient vs temperature. And now I need to find the Boltzmann temperature which is total_calka=0.

 

[Mark44]

By a root of an integrand i mean the place where the total_calka = 0.

That doesn't explain anything to me.

 

[D H]

That doesn't explain anything to me.

He's looking for a point x where \int_a^x f(t)\,dt = 0.

 

[Mark44]

He's looking for a point x where \int_a^x f(t)\,dt = 0.

x = a works

 

[FactChecker]

Some general programming notes:
1) Comments in code are your friends
2) So are variable names that mean something and don't have to be explained later. (Examples: What is beta and why is it set to 4/tempred? What is tempred?)
3) Debugging code can be done by putting write statements at significant points that print the values of variables there.
4) You should check the lower level procedures before you worry about the top level results. Are there parts of the code that you are confident of? Has the integration algorithm been verified?

If I was debugging this, I would insert more write statements in and run it. Then you can follow the progress of the calculation. The binary search algorithm should be easy to follow as it narrows in on a solution.