What is the Pseudo Code for a Riemann Sum Calculation?

[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

answer = s * 1/j



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.

 

[mathmate]

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

answer = s * 1/j

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
should read
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.

 

[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'
READ*, J
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

 

[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.

      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'
      READ*, J
      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.

 

[Hurkyl]

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.

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"

 

[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?

 

[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.

 

[ncp1044]

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

 

[mathmate]

Basically, you would start your programme. Inside the single programme, you would spilt up the work in, say, two parts, each of which should be independent of each other to a certain extent. In your case, one part would calculate the sum of the odd values of i, and the other part for the even values. You would initiate threads, or processes for each part, each of which will work independently of each other. The threads may or may not occupy the same processor. It's usually up to the operating system to decide. Most of the time, if you have 4 processors, you can be quite sure that each thread will occupy a processor. All is well until they finish the work, report to your main program with a subtotal to be added to the other one. Your main program will take control and assure that all the threads will have returned the answers before the addition is done. This coordination requires programming to synchronize work, as well as assuring that certain part of the work is done after the prerequisite work is finished. This is achieved through signals.
As I said, I do not know much about parallel processing, so this is about as much as I can tell you. If you google parallel processing or something like that, you will probably know more in an hour than I could ever tell you.
As I said, if you have less than 4 processors, don't even bother thinking about it, becaue the overhead will eat away most of your benefits. Parallel is very useful for problems which are massively parallel, like 1000 processors or more. Then the additional programming is well worth it, especially if you have to run the same programme all the time.

 

[ncp1044]

Perhaps I should buy a book on it. It sounds like a pretty dense topic. Thanks for pointing me in the right direction.