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

Click For Summary

Discussion Overview

The discussion revolves around a user's difficulty in obtaining the expected results from a Fortran implementation of the bisection method to find the root of an integrand, specifically where the variable total_calka equals zero. The scope includes programming issues, algorithm implementation, and debugging techniques.

Discussion Character

  • Technical explanation
  • Debate/contested
  • Mathematical reasoning

Main Points Raised

  • The user describes their bisection code and the variables involved, including tl, tup, tempred, and s, and expresses confusion about why the expected root is not found.
  • Some participants question the clarity of the user's code, noting the lack of parameters in the integration2 subroutine and the use of unclear variable names.
  • One participant clarifies that the user is looking for a point where total_calka = 0, while another suggests that the user is seeking a point x such that \int_a^x f(t)\,dt = 0.
  • General programming advice is provided, emphasizing the importance of comments, meaningful variable names, and debugging techniques such as inserting write statements to track variable values.
  • Participants express uncertainty about the user's explanation of the root of the integrand and the overall intent of the code.

Areas of Agreement / Disagreement

Participants generally agree on the need for clearer code and debugging practices, but there is no consensus on the specific issues causing the user's code to fail in finding the root.

Contextual Notes

The discussion highlights potential limitations in the user's code, such as unclear variable definitions and the absence of comments, which may hinder understanding and debugging. There is also a lack of clarity regarding the integration algorithm's verification.

Galizius
Messages
14
Reaction score
0
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:

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.

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
    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()
Code: 
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.
 
Technology news on Phys.org
Galizius said:
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:

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.

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
    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()
Code: 
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)
 
  • Like
Likes   Reactions: FactChecker
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: 
 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.
 
Galizius said:
By a root of an integrand i mean the place where the total_calka = 0.
That doesn't explain anything to me.
 
Mark44 said:
That doesn't explain anything to me.
He's looking for a point x where \int_a^x f(t)\,dt = 0.
 
D H said:
He's looking for a point x where \int_a^x f(t)\,dt = 0.
x = a works :oldbiggrin:
 
  • Like
Likes   Reactions: 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.
 
  • Like
Likes   Reactions: Mark44