# Pseudo code for reimann sum

1. Nov 28, 2008

### ncp1044

So I'm writing a program to do a riemann sum. The final answer should be Pi.

My function: f(x) = 4/(1+x^2) a = 0 b = 1 n = 40,000,000

For a riemann sum, delta x = (b-a)/n so (1-0)/40,000,000 = 1/40,000,000.

Then the riemann sum looks like:

(1/40,000,000) [f(1/40,000,000) + f(2/40,000,000) + f(3/40,000,000) + ...until you get f(40,000,000/40,000,000)]

That's the math behind the problem, now onto the pseudo code:

Enter limits (lower, upper)
lower = i upper = j
s = 0
while i =< j
s = s + 4/(1 + (i/j)^2)
i = i + 1

The limits mentioned in the code above are i = 1 and j = 40,000,000

I'm not very confident on my pseudo code; any feedback would be helpful.

2. Nov 28, 2008

### mathmate

If I am not mistaken, Riemann sum is the same as the trapezoidal rule, where
I=(0.5T0+T1+T2+....+Tn-1+0.5Tn)/n
In your example, T0 would be correspond to i=0, and tn to j=40,000,000
The number of terms is n+1, so the first and last should be divided by two to get the equivalent of n terms.
You can adjust the pseudocode according to the above. The line
lower=i, upper=j
i=lower, j=upper
Otherwise it looks OK to me.
If you are ready, you can start coding it in Fortran after the adjustments.
Be sure to define i and j as REAL*8, because if they happen to be integers, i/j will give zero as the answer.
Also, I suggest you start with n=10 or 20 just to test the program before looking for precision.

3. Nov 29, 2008

### ncp1044

Here is my code. It seems to work except two things. When the number of subintervals is around 2,000,000 or greater, it starts to slowly diverge from pi. Also, at 15,000,000 the sum equals 4, and at around 17,000,000 it doesn't return an answer.

I'm trying to use 40,000,000 subintervals and parallel computing when the final product is complete.

code:

PROGRAM SUM

REAL S,Pi,I,J
PRINT*, 'This program approximates Pi using a reimann sum.'
PRINT*, 'Enter the number of sub intervals you wish to use in'
PRINT*, 'the reimann sum'
I=0
S=0
5 IF(I.LE.J)THEN
S=S+4/(1+(I/J)**2)
I=I+1
GO TO 5
ENDIF
Pi=S*(1/J)
PRINT*,'The reimann sum is',Pi
END

4. Nov 29, 2008

### mathmate

I do not know which compiler you use, but you can improve by using REAL*8 instead of simply REAL, and also using integer for I and J. This can be done by
IMPLICIT REAL*8(A-H,O-Z), INTEGER*4(I-N)
There are a couple of places you'd have to make a conversion from integer to real.
The code is shown below, and it works well till 200000000 cycles.
Code (Text):
IMPLICIT REAL*8(A-H,O-Z), INTEGER*4(I-N)
C     REAL*8 S,Pi
C     INTEGER I,J
PRINT*, 'This program approximates Pi using a reimann sum.'
PRINT*, 'Enter the number of sub intervals you wish to use in'
PRINT*, 'the reimann sum'
I=0
S=0
5 IF(I.LE.J)THEN
S=S+4/(1+(I*1.0/J)**2)
I=I+1
GO TO 5
ENDIF
Pi=S*(1.0/J)
PRINT*,'The reimann sum is',Pi
END
Also, if you are interested in calculating pi, there are more convergent algorithms. For computing testing, this is OK.

5. Nov 29, 2008

### Hurkyl

Staff Emeritus
That's because there are two sources of error.

The first source of error is that the Riemann sum is only an approximation to the answer.
The second source of error is that you are not doing exact arithmetic. (e.g. you are only computing approximate values for + and *)

The second error gets larger as you increase the number of subintervals, because you are doing more operations and so there is more error to accumulate. Apparently, with the specific amount of precision you're using, 2,000,000 is around where the second error starts getting significant.

P.S. you're misspelling "Riemann"

6. Nov 29, 2008

### ncp1044

To Mathmate: Thanks! That's much better, and it's all accurate. And I'm using Force2.0, a freebie that I found on download.com.

To Hurkyl: The reason for the error makes sense now. And thanks for the heads-up, I feel stupid now haha. rIEmann.

Edit: What would it take to do this in parallel?

Last edited: Nov 29, 2008
7. Nov 29, 2008

### mathmate

I am not familiar with parallel processing. But how effective it would get depends on how many processors you have, and whether you have access to the control of threads.
In the back of my mind, you won't bother with special programming for two or three processors. The extra processors will take care of the backend, like file-IO, system maintenance, networks overhead, etc.
If you have control of the threads, you can always program to split up the work into two parts, one for each thread, for i=0 to J in steps of two, and for i=1 to J in steps of two.
It keeps both threads busy simultaneously, and you only have to add them together (when they are both ready) at the very end. I suppose you could do something similar if you have more processors.

8. Nov 30, 2008

### ncp1044

Can you expound on "If you have control for the threads..."?
I can't say I know anything about parallel processing either...

9. Nov 30, 2008