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

Riemann Sum with Fortran 90

  1. Apr 20, 2008 #1
    [SOLVED] Riemann Sum with Fortran 90

    My assignment: Use Reimann Sums to estimate pi to 6 decimal places (ie: you can stop when successive iterations yield a change of less than 0.000001. For the Reimann Sums solution, an iteration equals 2X the number of segments as the trial before. Print out the results and number of iterations.

    So, I'm trying the best way to do this. We're supposed to use a unit circle, so the equation for a circle will be 1 = (x^2+y^2)^.5

    I'm not that good at Fortran, but I've started part of it. This assignment is really confusing me, I just don't know if the start I have is even going to work. If there's anyone that can point me in the right direction, that'd be great.

    PROGRAM riemann
    IMPLICIT NONE

    REAL,PARAMETER::lower=0.0, upper=1.0
    REAL::dx,x,y
    integer::i

    y=0.0
    x=0.0
    dx=(upper-lower)/i

    DO i=1,i*2

    y=(1-x**2)**.5
    x=x+dx
    ?
    ?
    END DO
    END PROGRAM
     
  2. jcsd
  3. Apr 20, 2008 #2
    Note: I haven't read your code at all, I think that it's for your to figure out the syntax. I just want to help you wrap your brain around what you have to do.

    Let's see... I forgot how exactly to do Riemann sums, but the general rule is to split your area into discrete bins to add up, right?

    So, what I would do is split the unit circle into 2 semi-circles, starting from the Origin and going all the way out to x = 2. That way you have a length, right? Or you could just take 1 semi-circle and double it. :smile:

    So you can continually iterate and say something like:

    binSize = length/numberOfBins

    Your length will likely be 2 if you just use a unit circle and your numberOfBins will be 2^(Iterations) power, since you want to double at every iteration. So now you know how long each bin will be for a given iteration. You have your dx. It should be constant per iteration. You get dy explicitly by setting y= root(1-x^2) or whatever. You have that already. You'll need two of those, or simply doubling your semi-circle.

    Anyway, you then take two values of y, and a value for dx. You can now use the formula for the area of a trapezoid.

    So, in pseudo-code it would look like this:

    x = 0
    y = root(1 - x^2)
    dx = length/numberOfBins
    totalArea = 0

    for (i = 0, i < number of bins, i++)
    {
    y_1 = y
    x += dx
    y_2 = y

    area = x * [(y_1 + y_2)/2]

    totalArea += area
    }
    Something like that.

    Then it's a matter of dividing the area by your radius-squared. So really, it's just the total area you get.
     
    Last edited: Apr 20, 2008
  4. Apr 21, 2008 #3
    !riemann sums from the left

    program riemann_left

    implicit none

    real :: x, approx, del_x
    integer :: i, n_points

    print *, "how many points?"
    read *, n_points

    x = 0

    del_x = 1/n_points
    approx = 0

    do i=1,n_points-1

    approx = approx + del_x*sqrt(1-(x+del_x)**2)
    x = x + del_x

    end do

    print *, "pi=",approx*4

    end program riemann_left


    there's a mistake somewhere i'll find it when i get back from my final
     
    Last edited: Apr 21, 2008
  5. Apr 21, 2008 #4
    !riemann sums from the left

    program riemann_left

    implicit none

    real :: x, approx, del_x
    integer :: i, n_points

    print *, "how many points?"
    read *, n_points

    x = 0

    del_x = 1./n_points
    approx = 0

    do i=1,n_points-1

    approx = approx + del_x*sqrt(1.-x**2)
    x = x + del_x
    x
    end do

    print *, "pi=",approx*4.

    end program riemann_left


    mixed arithmetic isn't good
     
  6. Apr 24, 2008 #5

    rcgldr

    User Avatar
    Homework Helper

    What's missing so far here is a variable for the number of segments. This will double on each iteration as described. You'll need to compare the new approximation with the previous one, and if it differs >= .000001, copy the new appxoximation value into a variable to hold the previous approximation and do another interation. If it differs by < .000001, then the program is done and should print out the approximation and the number of segments.

    In this case, there aren't any actual bins or segements, since you're summing as you calculate each value. This will work fine as long as there isn't a huge difference in magnitude between the sum and each individual value. Once this becomes an issue due to truncation, you'll need a summing algorithm, but that's beyond the scope of the stated problem and shouldn't be an issue if you're using double precision and only going to 7 digits.

    If you're curious, summing algorthims use actual bins. One simple method is for the bins to hold: a number, sum of 2 numbers, sum of 4 numbers, sum of 8 numbers, ... sum of 2^n numbers. This summing function would get a new number as input each time, and by using a static counter, would keep track of what steps to peform in order to update and/or zero out bins. A final call would be made to sum up all the bins to produce the total sum. Note that the array is initialized to zeroes before summing starts.

    Another method is to only add numbers of equal magnitudes, by using the exponent fields of the numbers and intermediate sums to index an array of size 2^e elements, where e is the number of bits in the exponent of the floating points numbers. For double precision, this requires 2048 bins (2 of these are speicals, Nan Infinity use the same exponent, and zero), since the exponent is 11 bits. Upon receiving a number it the array is indexed, if the value in the array is zero, the number is just stored, if not, the number and the array value are added to produce a new intermediate sum, the array entry is zeroed, and the process repeated using the new intermediate sum until a zero entry is found or overflow occurs. This consumes more time than the "binary" algorithm above, but requires less space. Again, a final call is made to sum up all the bins to produce the total sum, and the array is initialized to zeroes before summing starts.
     
    Last edited: Apr 24, 2008
  7. Apr 24, 2008 #6

    rcgldr

    User Avatar
    Homework Helper

    Correction, this consumes more time and space (assuming that log2(n) is less than 2048), but is more accurate (less truncation).
     
  8. Apr 28, 2008 #7
    So what exactly should I use for that first approximation, since it won't have anything prior to compare itself to?
     
  9. Apr 29, 2008 #8

    rcgldr

    User Avatar
    Homework Helper

    Since you know the target number is pi, just start with 0.000000 or any number sufficiently different than pi.

    If this was a more generic program, then you would need check for iteration < 2 before checking if the difference is < .000001.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Riemann Sum with Fortran 90
  1. Fortran 90 (Replies: 7)

  2. Fortran Riemann Sum (Replies: 2)

  3. Fortran 90 (Replies: 1)

  4. Fortran 90 (Replies: 1)

Loading...