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

FORTRAN program for harmonic series sum

  1. Jul 20, 2012 #1
    Hi

    Here is my code to get the sum of harmonic series. Harmonic series is

    [tex]\sum_{i=1}^{\infty}\; \frac{1}{i} [/tex]

    here is the code

    Code (Text):
    program harmonic

    implicit none

    integer :: i,n
    real :: sum=0.0

    write(*,*)'How many terms you want to sum ?'
    read(*,*) n


    do i= 1 ,n
       

          sum=sum +(1.0/real(i))
         

    end do

    write(*,*) 'The sum of series is = ', sum


    end program harmonic
    Now if I give n= 10000000 , it gives me the answer 15.403683. But Mathematica gives me answer 16.6953. So there is something wrong happening in this simple code.

    Can anybody point out the mistake ?
     
  2. jcsd
  3. Jul 20, 2012 #2

    Mark44

    Staff: Mentor

    I suspect that the problem is that your sum variable is of type real (or single precision), which is a four-byte type. See if changing the type of this variable to double precision, and write 1.0 as 1.0D0. You'll also need to convert your loop counter, i, to double precision, which you can do with the DBLE intrinsic function.
     
  4. Jul 20, 2012 #3

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    I think Mark44 is probably right.

    But you might try the experiment of summing the terms smallest first, not largest, i.e. change
    do i = 1,n
    to
    do i = n,1,-1
    (and stay with real variables, not double precision).

    If this "works", figuring out why it works is left as an exercise for you (and it's something important to know about floating-point calculations in general, not just for this situation).

    FWIW A good approximation should be ##\ln(10000000) + \gamma## where ##\gamma## is the http://en.wikipedia.org/wiki/Euler–Mascheroni_constant which is the same as the Mathematica answer.
     
    Last edited: Jul 20, 2012
  5. Jul 20, 2012 #4
    I would definitely start the loop with the largest value; otherwise, the little numbers are simply going to fall off.

    The typical fortran integer is 4 bytes, so, no problem there for i or n to accomodate 10 million.

    Actually, the problem is with the reals...they are also typically 4 bytes and maybe not large enough to accomodate little numbers; so, for the real, I would go double precision, like real*8.

    Then, I would cast the loop variable right away and use the declared real variable in the division, like this
    Code (Text):

    program harmonic
    implicit none
    integer :: i,n
    real*8 :: s,sum=0.0

    write(*,*)'How many terms you want to sum ?'
    read(*,*) n

    do i= n,1,-1
        s = i
        sum = sum + 1.0/s
    end do
    write(*,*) 'The sum of series is = ', sum

    end program harmonic
     
    This way there is no need to be writing 1.0D0; instead, 1.0 would automatically be promoted to the same type as 's' in the denominator...which we already declared to be REAL*8
     
  6. Jul 21, 2012 #5
    thanks everybody....

    when I use the following code as suggested by alephzero ,

    Code (Text):
    program harmonic

    implicit none

    integer :: i,n
    real :: sum=0.0

    write(*,*)'How many terms you want to sum ?'
    read(*,*) n


    do i= n ,1, -1
       

          sum=sum +(1.0/real(i))
         

    end do

    write(*,*) 'The sum of series is = ', sum


    end program harmonic
    For n=10000000, it gives now the answer 16.686031.

    Then I tried the solution suggested by gsal

    Code (Text):
    program harmonic
    implicit none
    integer :: i,n
    real*8 :: s,sum=0.0

    write(*,*)'How many terms you want to sum ?'
    read(*,*) n

    do i= n,1,-1
        s = i
        sum = sum + 1.0/s
    end do
    write(*,*) 'The sum of series is = ', sum

    end program harmonic
    now for n=10000000, I get the sum as 16.695311365859890 , which is very close
    to the Mathematica answer, 16.6953. I don't know how to get more precise answers
    in Mathematica. I used the command

    Code (Text):
    Sum[1.0/i, {i,10000000}]
    Next, I tried to alter gsal's code as

    Code (Text):
    program harmonic
    implicit none
    integer :: i,n
    real*8 :: s,sum=0.0

    write(*,*)'How many terms you want to sum ?'
    read(*,*) n

    do i= 1, n
        s = i
        sum = sum + 1.0/s
    end do
    write(*,*) 'The sum of series is = ', sum

    end program harmonic
    Now here, for n=10000000, I get the sum as 16.695311365856710 , which is different in
    last 4 digits from the answer from gsal's code.

    So using double precision variables certainly improves the precision , but the order in which
    we sum the series seems to do some funny things......As commented by alephzero,
    summing smallest terms first is probably more accurate. He has asked me to think over it
    as an exercise. I found something on yahoo answers

    Here, its suggested that

    So how would I do this kind of splitting. Also I want to time my code to get the time the compiler takes to compute the result. I have Win XP as OS and I am using Ubuntu inside
    VirtualBox. Compiler is gfortran......

    thanks
     
  7. Jul 21, 2012 #6
    If you are going to get into timing the program, that's another matter which I personally am not that good at...I am talking about setting things up according to what task takes longer, moving things in and out of the register, vectorizing loops, etc.

    By the way, not every decimal number can be represented in binary with a given number of bits and so, just because you keep getting a larger and larger total sum, it does not mean it is correct...the correct number could very well be smaller than those large numbers, it is just that some of the resulting 1.0/i quantities may have had to be represented in binary in a number that when converted back to decimal is larger than 1.0/i .....do you get this?

    For example, if you have a 3-bit binary system, you can exactly represent the numbers 0 thru 7...but how can you represent 5.5? You can't, you would have to chose either 101 or 110 !!! Neither one of which converts back to 5.5

    Also, you are going to have to choose between speed and accuracy. In Fortran, you can work with various precisions upon request...read up on "selected_real_kind".

    Then, as mentioned above, there is the more elemental thing about moving things in and out of the registry for addition, leaving it there during the entire loop or what.


    Example 1:

    Code (Text):

    program harmonic
    implicit none
    integer :: i,n,d(5)
    real*8 :: s,rsum(5),total

    !write(*,*)'How many terms you want to sum ?'
    !read(*,*) n
    n=10000000

    d(1) = 1*n/5
    d(2) = 2*n/5
    d(3) = 3*n/5
    d(4) = 4*n/5
    d(5) = 5*n/5

    rsum=0.0
    do i= 0,n/5-1
        rsum(1) = rsum(1) + 1.0/(d(1)-i)
        rsum(2) = rsum(2) + 1.0/(d(2)-i)
        rsum(3) = rsum(3) + 1.0/(d(3)-i)
        rsum(4) = rsum(4) + 1.0/(d(4)-i)
        rsum(5) = rsum(5) + 1.0/(d(5)-i)
    end do

    total = sum(rsum)

    write(*,*) 'The sum of series is = ', total

    end program harmonic
     
    needless to say, the code above assumes that the number given can be nicely divided by 5...besides that, though:
    • it has a single loop
    • the loop needs to be executed only n/5 times
    • BUT there are 5 different operations in the loops which require the various variables to be moved in and out of the registry.

    Then, there is:
    Code (Text):

    program harmonic
    implicit none
    integer :: i,n,d(5)
    real*8 :: s,rsum(5),total

    !write(*,*)'How many terms you want to sum ?'
    !read(*,*) n
    n=10000000

    d(1) = 1*n/5
    d(2) = 2*n/5
    d(3) = 3*n/5
    d(4) = 4*n/5
    d(5) = 5*n/5

    rsum=0.0
    do i= 0,n/5-1
        rsum(1) = rsum(1) + 1.0/(d(1)-i)
    end do
    do i= 0,n/5-1
        rsum(2) = rsum(2) + 1.0/(d(2)-i)
    end do
    do i= 0,n/5-1
        rsum(3) = rsum(3) + 1.0/(d(3)-i)
    end do
    do i= 0,n/5-1
        rsum(4) = rsum(4) + 1.0/(d(4)-i)
    end do
    do i= 0,n/5-1
        rsum(5) = rsum(5) + 1.0/(d(5)-i)
    end do

    total = sum(rsum)

    write(*,*) 'The sum of series is = ', total

    end program harmonic
     
    where:
    • you have several loops
    • but each loop is only executed n/5 times
    • and the quantities inside the loop do not have to keep getting moved in and out of the re
      gistry during the entire execution of the loop

    Then, of course, you could start taking advantage of Fortran 90 features about vectorization and increased precision; the following code:

    Code (Text):

    program harmonic
    implicit none
    integer :: i,n,counter
    integer, parameter :: rkind = selected_real_kind(20)
    real(kind=rkind) :: total, rsum(10000000)

    n=10000000
    forall (i=1:n) rsum(i) = 1.0_rkind/i
    total = sum( rsum((0*n/5+1):(1*n/5)) ) + &
            sum( rsum((1*n/5+1):(2*n/5)) ) + &
            sum( rsum((2*n/5+1):(3*n/5)) ) + &
            sum( rsum((3*n/5+1):(4*n/5)) ) + &
            sum( rsum((4*n/5+1):(5*n/5)) )
    write(*,*) 'The sum of series is = ', total

    end program harmonic
     
    for example, takes quite a while to run, but produces the following number:

    16.69531136585985181539911893954098

    ...then, again, I don't quite know what I am doing...using your post as an excuse, I went ahead and tried to use both "forall" and "selected_real_kind" for my first time! Be warned.
     
  8. Jul 21, 2012 #7

    uart

    User Avatar
    Science Advisor

    Hi gsal. Did you compare the accuracy of the results obtained by splitting the code into 5 separate sums, compared with just using just the one (summing terms from smallest to largest as suggested by AlphZero of course). In my own tests I couldn't notice any improvement by splitting the sum. Increasing the real precision of course makes a big difference.

    For reference, to about 80 dp the results of f(10000000) is:
    16.69531136585985181539911893954045188424986975237308046278513595435628869217425467

    Hi issacn.
    If you're interested in both speed and precision then it is worth considering that the highest precision that is supported natively by your CPU is probably 80 bits (10 bytes).

    The syntax seems to be compiler dependent, in ftn95 use "real (kind=3) :: sum=0" for example. And in gnu g95 use "real (kind = 10) :: sum=0" for example.

    As far as I can tell it seems that ftn95 treats "kind" in a simple ordinal way, with kinds 1,2,3 representing the supported precisions in order from lowest to highest (single, double, extended). Whereas the gnu "g95" compiler appears to use kind to specify the number of bytes used.
     
    Last edited: Jul 21, 2012
  9. Jul 21, 2012 #8
    Thank you gsal and uart for the explanation. I am following Chapman's "Fortran 95/2003
    for scientists and engineers" third edition. And I am currently on the 4th chapter about loops.
    So I wanted to apply that immediately to do the sum of this famous series. As a result of this
    thread , I have learned that there are various issues involved here. So may be I need to wait till I finish this book. Author doesn't discuss "kind" till , I think , chapter 11. But all the posts in this thread were very informative. The world of high computing seems very exciting......

    Thanks
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: FORTRAN program for harmonic series sum
  1. Fortran Riemann Sum (Replies: 2)

  2. Programming in fortran (Replies: 5)

Loading...