November 8, 2007 --- Class 21---Tar; Testing Trapeziod, Simpson and
Romberg Methods; Hit or Miss Method
Activities:
Tar
I did my work on the tests of integration methods on my linux
computers. I am going to show you a convenient way to move
a group of files to another computer using the unix tar command.
Tar stands for tape archive. I don't have a tape drive anymore,
but tar can also create files on disk that can be moved with scp.
To create a tar file, use
tar -cvf myarchive.tar files
or
tar -cvf myarchive.tar directory
c is for create, v is for verbose and f indicated that the next
argument is the name of the tape archive file to be created. Then
comes a list of files and/or directories that are to be include in the
archive. To extract the files from the archive, you may use the
command:
tar -xvf myarchive.tar
This unpacks the files from myarchive.tar and puts them in the current
directory. If you used the second form above, then a directory will
be created and the files placed there.
Testing Trapeziod, Simpson and Romberg Methods
In ~sg/integration_test/double_precision, you will find source code
and output files. The source code has been modified so we can see
how the method converges. We went over the source code in detail.
I also showed how to compile the code as the executables in the
tar file were for linux, not Mac OS.
The integral of cos(x) from x=0 to x=pi is calculated by xqsimp2.
Here is the output:
1 7.853981633974e-01 3.333333333333e+29
2 9.480594489685e-01 1.002279877492e+00
3 9.871158009728e-01 1.000134584974e+00
4 9.967851718862e-01 1.000008295524e+00
5 9.991966804851e-01 1.000000516685e+00
6 9.997991943200e-01 1.000000032265e+00
7 9.999498000921e-01 1.000000002016e+00
8 9.999874501175e-01 1.000000000126e+00
9 9.999968625353e-01 1.000000000008e+00
10 9.999992156342e-01 1.000000000000e+00
The first column is the level at which we call the routine, i.e., with
2^(j-1) intervals in the integration range. The second column is the
trapzoid rule and the third column is for Simpson's rule. Note that
the j=1 value for Simpson's rule is bogus since we need to go at least
to j=2 to use Simpson's rule. In cos.ax, you will find the error
plotted vs. j. The crosses are for the trapezoid rule, the diamonds
are for Simpson's rule. You will also find a single point from
the Romberg routine.
We found that from j=0 to j=10 the value for the trapezoid rule
fall by 6 orders of magnitude. This corresponds to increasing the
number of intervals by 2^10 or 1024. Since the error for the
trapzoid rule is supposed to fall like n^{-2}, and 1024^2 is about
10^6, this makes perfect sense.
Simpson's rule is seen to converge much more quickly and falls
12 orders of magnitude. This is the expected behavior. With the
Romberg method, the error is about 10^{-12} for j=5. Quite an
impressive improvement!
A second example integrates sqrt(1-x^2) from 0 to 1. Curiously, it
looks very different. This is because of the singularity in the
derivative of the integrand at x=1. Simpson's rule does not converge
with a higher power of j. Plot sqrt_1minus_xsqr.ax or plot both graphs
together. We recalled that the expected error for the trapezoid rule
depends on the average value of the second derivative of the function.
Testing Hit-or-Miss Method
Program ~sg/src/pi3.c implements the hit-or-miss method using a
feedback shift register random number generator.
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.
I ran the pi3 code a few times with this command
pi3 10
The output ranged from 2.4 to 3.6. Clearly the output is a random
variable. We have to answer several questions.
Do we get the right answer on average?
Can we estimate the error?
How does the error depend on n, the number of trials? We expect
the error to fall like 1/sqrt(n), but can we show that?
I asked everyone to create datafiles with 100 batches of various
trial sizes: 10 20 40 80 160. I did not use a shell script, and used
a interactive foreach command instead. I then used another foreach
loop to make histograms from each datafile. In our next class
we will study the statistics further.