November 8, 2005 --- Class 20 --- Simpson's Rule Rederived, Hit or Miss Method
Activities:
Simpson's Rule Rederived
In the last class, starting from the trapezoid rule, we used
the formulae below for correcting the leading 1/N^2 errors.
S_N = A + E/N^2
S_2N = A + E/(2N)^2
so A = 4/3 S_2N -1/3 S_N
Applying this to a one-dimensional itegral, we easily found that
the weighting of the sum of two expressions on the right hand side
exactly coincides with the Simpson's rule formula. We previously
found that Simpson's Rule has a leading error of order 1/N^4. Since
the approach above removes the leading 1/N^2 error and since the
error of the trapezoid rule only has even powers of 1/N (asserted, but
not proved), this is not surprising.
Hit-or-Miss Method
The Hit-or-Miss method and Sample Mean methods are described in
Section 11.2. Let's try the first method and see how well it does.
Program ~sg/src/pi1.c implements the hit-or-miss method using a
random number generator that was in the first edition of Numerical
Recipes. It is called ran1 and the source is in ~sg/src/ran1.c.
You have to tell the program how many "darts you want to throw."
This is the first command line argument. You may also have it
repeat the experiment multiple times. This is the optional second
argument. Look at the code to see how this is done if you are
unfamiliar with this type of feature. Also note how the random
number generator gets its seed through a call to pid(). This is
the "process id" of the calling routine. It is essentially a
random number. If you like, you can modify the code so that you
specify the seed yourself. This can be very useful when you are
trying to debug the code and want to see if results are
reproducible. Sometimes I set up my code so that it asks for a
seed, but if the seed I give is negative, it calls pid() to get a
random seed. That way, I can do many runs without having to think
up new seeds, and when I want to reproduce runs with a known seed I
can also do that.
Call pi1 with the argument 50 and 500, and send the output to a
file. Each line will have 50 (the number of darts) and the
estimate of pi/4 (the value of the integral).
Hist
It is interesting to look at the distribution of the estimates. One
way to do this is by making a histogram. You break the x-interval
up into different ranges and then count the number of times a
value falls in each range. hist is a program for doing this. I
find it very useful. You are welcome to copy the source code from
~sg/src/hist.c. If you have a file with a bunch of numbers to
histogram, just type hist.c