|
|
This notebook is designed to accompany Chapter 4 of "Introduction to Scientific Programming: Computational Problem Solving Using Mathematica and C" by Joseph L. Zachary. In it, we will use Mathematica to explore the block-stacking problem. (26Jul97}
To use this notebook you will need to use some Mathematica extensions that we have created. Load our blocks package by evaluating the Mathematica command below. (You will need to have first installed our custom Mathematica library.)
Needs["ISP`Blocks`"]
The four functions defined in the library are described in Chapter 4 of the text. Each takes as a parameter the number of blocks in a stack of blocks arranged as described in the text and reports back the amount by which the top block extends beyond the table. For example, let's compute (in four different ways) the extension that can be achieved with 100 blocks.
The "BlockFloat" function uses floating-point arithmetic to explicitly sum the series 1/2 + 1/4 + 1/6 + 1/8 + ... + 1/2n. For n = 100, that sum is
BlockFloat[100]
The "BlockRat" function uses rational arithmetic to sum the same series. After it has computed an exact sum, we use "N" to convert the result into a floating-point number.
N[BlockRat[100]]
The "BlockFast" function computes the sum of the series, not by adding up each term, but by using the formula (1/2) (Psi(n+1) + gamma). Notice that we get the same answer as with "BlockFast".
BlockFast[100]
The "BlockPrecision" function takes a second parameter "p". It uses simulated floating-point arithmetic with precision p (if p<=10); it uses arbitrary-precision arithmetic with precision p (if p>$MachinePrecision); othewise, it uses regular floating-point numbers. For example, here we compute the extension possible with 100 blocks using 3-digit floating-point numbers
BlockPrecision[100, 3]
"BlockPrecision" can be used to study how roundoff error can accumulate in a series of computations. Let's compare the results of using 4-digit floating-point numbers, rational numbers, and the formula embodied in "BlockFast" to sum increasingly larger numbers of blocks.
StyleBox[{BlockPrecision[100, 4], N[BlockRat[100]], BlockFast[100]}, ShowStringCharacters -> True]
StyleBox[{BlockPrecision[200, 4], N[BlockRat[200]], BlockFast[200]}, ShowStringCharacters -> True]
StyleBox[{BlockPrecision[400, 4], N[BlockRat[400]], BlockFast[400]}, ShowStringCharacters -> True]
StyleBox[{BlockPrecision[800, 4], N[BlockRat[800]], BlockFast[800]}, ShowStringCharacters -> True]
StyleBox[{BlockPrecision[1600, 4], N[BlockRat[1600]], BlockFast[1600]}, ShowStringCharacters -> True]
Notice what happens:
BlockPrecision, which operates by directly summing the series using 4-digit floating-point arithmetic, suffers more and more from roundoff error as the number of terms increases.
BlockRat, which operates by directly summing the series using rational number arithmetic, always returns an exact answer. We then convert that exact answer into a floating-point number that is exact to the available number of digits.
BlockFast, which operates by exploiting specialized mathematical knowledge about the Harmonic series, returns a number that is exact to the available number of digits. BlockFast is also extremely fast.
Let's go back to 16-digit mantissas and explore how roundoff error accumulates when doing a long series of floating-point calculations.
Using BlockFast, we can find the extension possible with one thousand blocks. (We save the result in the variable "truth".)
truth = BlockFast[1000]
Let's use BlockPrecision to do the same calculation using mantissa lengths ranging from 1 to 10. For example, with a five-digit mantissa, the result is as follows. (We save the result in the variable "result".)
result = BlockPrecision[1000, 5]
The relative error is approximately 0.024%.
Abs[truth - result] / truth
You should do similar calculations using for each mantissa length from 1 to 10. When you are done, you will know the relative error for each mantissa length that is exhibited by BlockPrecision when calculating the extension for one thousand blocks.
You can visualize the effect of mantissa length on relative error by plotting mantissa length on the x-axis against relative error on the y-axis. To do this, use the command below. It plots the numbers 1 through 10 against some arbitrarily chosen numbers. You will have to replace the arbitrarily chosen numbers with the relative errors that you have calculated. For example, [5, 30] would become [5, .00024]. (Rounding the relative errors to two decimal places is OK.)
ListPlot[{{1, 88}, {2, 22}, {3, 14}, {4, 13}, {5, 30}, {6, 0}, {7, 2}, {8, 99}, {9, 17}, {10, 88}}, AxesLabel -> {"Mantissa Length", "Error"}, PlotStyle -> PointSize[0.02], PlotRange -> {{0, 10}, {0, 100}}]
It is also interesting to compare the amount of time required by BlockFloat and BlockRat to do their calculations. For block stack sizes of 100, 200, 400, 800, 1600, and 3200, let's calculate the time required to calculate the extension using both BlockFloat and BlockRat. For example, the time required to compute BlockFloat[100] is the first number displayed by this expression. (The second is the result of the call to BlockFloat.)
Timing[BlockFloat[100]]
The Timing function is so crude that for quick calculations, it is a good idea to do several calculations and calculate an average. To obtain a better timing result for blockFloat[100], for example, we can do eight identical calculations and then divide the overall result by 8. (You'll have to ignore the calculated extension, however.)
Timing[BlockFloat[100] ; BlockFloat[100] ; BlockFloat[100] ; BlockFloat[100] ; BlockFloat[100] ; BlockFloat[100] ; BlockFloat[100] ; BlockFloat[100]]/8
For long calculations such as BlockFloat[3200], of course, this would be unnecessary .
You should time BlockFloat for stack sizes of 100, 200, 400, 800, 1600, and 3200. You can plot your results by completing the following command. You will need to insert the times that you obtain in place of the zeroes.
ListPlot[{{100, 5}, {200, 5}, {400, 5}, {800, 5}, {1600, 5}, {3200, 5}}, AxesLabel -> {"Stack Size", "Time"}, PlotStyle -> PointSize[0.02], PlotRange -> {{0, 3200}, {0, 10}}]
Repeat the timing experiments using BlockRat instead of BlockFloat. Plot your results by completing the command below.
ListPlot[{{100, 5}, {200, 5}, {400, 5}, {800, 5}, {1600, 5}, {3200, 5}}, AxesLabel -> {"Stack Size", "Time"}, PlotStyle -> PointSize[0.02], PlotRange -> {{0, 3200}, {0, 10}}]
Compare the shapes of the two curves that you create and try to explain the difference.