# Bisection formula fortran90/95

Tags:
1. Jun 25, 2015

### 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:

Code (Fortran):
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.

Code (Fortran):
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
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
tot_old=total_calka
enddo
write(1,*)tempred,-2*pi*total_calka
close(1)
end subroutine integration2
end module int2
There is also a function inte()
Code (Text):
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.

2. Jun 25, 2015

### Staff: Mentor

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)

3. Jun 25, 2015

### 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

Code (Text):
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),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.

4. Jun 25, 2015

### Staff: Mentor

That doesn't explain anything to me.

5. Jun 25, 2015

### D H

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

6. Jun 26, 2015

### Staff: Mentor

x = a works

7. Jun 26, 2015

### 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.