### Introduction to Scientific ProgrammingComputational Problem Solving Using:Maple and CMathematica and C

Author:
Joseph L. Zachary
Online Resources:
Maple/C Version
Mathematica/C Version

# Cantilevered Blocks Worksheet

Click below to download a Maple V worksheet. You can look at the appended non-interactive HTML version of the worksheet to learn what the worksheet covers.

This worksheet is designed to accompany Chapter 4 of Introduction to Scientific Programming: Computational Problem Solving Using Maple and C by Joseph L. Zachary. In it, we will use Maple to explore the block-stacking problem. (08Oct96)

# Getting Started

To use this worksheet you will need to use some Maple extensions that we have created. Read in our blocks package by evaluating the Maple command below. (You will need to have first installed our custom Maple library and configured Maple to use it.)

> with(blocks);

# Calculating Extensions

The three functions defined in the library are described in Chapter 4 of the text. Each takes as its 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 three different ways) the extension that can be achieved with 100 blocks.

The "blockFloat" function uses floating-point arithmetic to explicitly sum the series . 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 "evalf" to convert the result into a floating-point number. The difference between this sum, and the one calculated with "blockRat", can be attributed to roundoff error.

> evalf(blockRat(100));

The "blockFast" function computes the sum of the series, not by adding up each term, but by using the formula . Notice that we get the same answer as with "blockFast".

> blockFast(100);

"blockFloat" can be used to study how roundoff error can accumulate in a series of computations. Let's begin by setting Digits to 4, which means that we will be computing with four-digit floating-point mantissas.

> Digits := 4;

We can now compare the results of calling the three functions for increasingly large numbers of blocks.

> blockFloat(100);

> evalf(blockRat(100));

> blockFast(100);

> blockFloat(200);

> evalf(blockRat(200));

> blockFast(200);

> blockFloat(400);

> evalf(blockRat(400));

> blockFast(400);

> blockFloat(800);

> evalf(blockRat(800));

> blockFast(800);

> blockFloat(1600);

> evalf(blockRat(1600));

> blockFast(1600);

Notice what happens:

• blockFloat, which operates by directly summing the series using 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.

• blockRat is much slower than blockFloat for large series, but neither is practical for dealing with extremely large series. blockFast, on the other hand, is extremely fast.

## Exercises

Let's go back to 10-digit mantissas and explore how roundoff error accumulates when doing a long series of floating-point calculations.

> Digits := 10;

Using blockFast, we can find the extension possible with ten thousand blocks. (We save the result in the variable "truth".)

> truth := blockFast(10000);

Let's use blockFloat to do the same calculation using mantissa lengths ranging from 1 to 10. For example, with a five-digit mantissa,

> Digits := 5;

the result is as follows. (We save the result in the variable "result".)

> result := blockFloat(10000);

and after going back to ten digits

> Digits := 10;

we find the that relative error is approximately 2.8%.

> 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 blockFloat when calculating the extension for ten 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, .0277]. (Rounding the relative errors to three decimal places is OK.)

> plot([[1, 88], [2, 22], [3, 14], [4, 13], [5, 30],

> [6, 0], [7, 2], [8, 99], [9, 17], [10, 88]],

> 0..10, style=POINT, symbol=BOX,

> title = `Relative error when computing blockFloat(100)`,

> labels = [`Mantissa length`, `Error`]);

It is also interesting to compare the amount of time required by blockFloat and blockRat to do their calculations. Let's begin by going back to ten-digit mantissas.

> Digits := 10;

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 second number displayed by this sequence of commands. (The first is the result of the call to blockFloat.

> start := time(): blockFloat(100); time() - start;

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.

> start := time(): blockFloat(100): blockFloat(100): blockFloat(100): blockFloat(100): blockFloat(100): blockFloat(100): blockFloat(100): blockFloat(100); (time() - start) / 8;

For long calculations such as blockFloat(3200), of course, this would be unnecessary and would take an extremely long time.

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.

> plot([[100, 0], [200, 0], [400, 0], [800, 0], [1600, 0], [3200, 0]],

> 0..3200, style=POINT, symbol=BOX,

> title = `Time in seconds required by blockFloat`,

> labels = [`Stack size`, `Time`]);

You should now repeat the timing experiments using blockRat instead of blockFloat. Plot your results by completing the command below.

> plot([[100, 0], [200, 0], [400, 0], [800, 0], [1600, 0], [3200, 0]],

> 0..3200, style=POINT, symbol=BOX,

> title = `Time in seconds required by blockRat`,

> labels = [`Stack size`, `Time`]);

Compare the shapes of the two curves that you create and try to explain the difference.