Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Bisection formula fortran90/95

  1. Jun 25, 2015 #1
    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
        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 (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. jcsd
  3. Jun 25, 2015 #2

    Mark44

    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)
     
  4. Jun 25, 2015 #3
    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):: 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.
     
  5. Jun 25, 2015 #4

    Mark44

    Staff: Mentor

    That doesn't explain anything to me.
     
  6. Jun 25, 2015 #5

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    He's looking for a point [itex]x[/itex] where [itex]\int_a^x f(t)\,dt = 0 [/itex].
     
  7. Jun 26, 2015 #6

    Mark44

    Staff: Mentor

    x = a works :oldbiggrin:
     
  8. Jun 26, 2015 #7

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Bisection formula fortran90/95
  1. Fortran90 help (Replies: 4)

  2. Bisection Method (Replies: 2)

  3. Fortran 95 (Replies: 12)

  4. Fortran90 help (Replies: 5)

  5. Function problem (Replies: 8)

Loading...