# Riemann Sum with Fortran 90

• Fortran
[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

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. 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:
!riemann sums from the left

program riemann_left

implicit none

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

print *, "how many 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:
!riemann sums from the left

program riemann_left

implicit none

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

print *, "how many 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

rcgldr
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:
rcgldr
Homework Helper
Another method is to only add numbers of equal magnitudes ... This consumes more time than the "binary" algorithm above, but requires less space.
Correction, this consumes more time and space (assuming that log2(n) is less than 2048), but is more accurate (less truncation).

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.

So what exactly should I use for that first approximation, since it won't have anything prior to compare itself to?

rcgldr
Homework Helper
So what exactly should I use for that first approximation, since it won't have anything prior to compare itself to?
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.