Numbers

M. E. Mauel
Columbia University

Introduction

In this notebook, we'll learn how Mathematica represents numbers.
Mathematica has three basic types of numbers:
1. Integer. These are exact numbers. 
2. Real or "MachinePrecision". These are approximate, floating-point numbers in the representation used by your computer. Usually, desktop computers use 64-bit, or "double-precision" numbers. The results produced by Mathematica using Machine-Precision numbers can change from computer-to-computer.
3. Arbitrary-Precision Numbers. These are approximate, floating-point numbers in a special, representation, internal to Mathematica. With arbitrary-precision numbers of a given precision, p, the first p significant digits are known, and Mathematica considers there to be an "infinite number" of unknown digits afterwards.

The use of these numbers are demonstrated by showing:
1. Computation with Real, MachinePrecision, numbers are very fast. (Although not as quickly as program produced with an optimizing compiler!)
2. Since Mathematica uses your computer's processor when calculating MachinePrecision expressions, there may be considerable numerical error. Beware!
3. Computation with Arbitrary-Precision numbers are very slow. Beware!

In this course, we will be pushing Mathematica to its computational limits. For this reason, most (if not all) of your calculations will be Real, MachinePrecision, numbers. You will work to make your calculations run as quickly as possible--if for no other reason that you will not want to wait to see your results!

Off[General::"spell1"];

Typing Numbers

Integer

11^24

Precision[1]

13^7/14^6 + 1 + Sin[45]

Integers are exact numbers! Mathematica will simplify expressions, but it will not make approximations.

Real (MachinePrecision)

11.0^24

1.23*^-6

Precision[1.0]

?*Machine*

$MachinePrecision

{$MaxMachineNumber, $MinMachineNumber}

$MachineEpsilon

i = 1.0;
While[1.0 + i != 1.0, i = i/2]; i

{1.0 == 1.0 + i, 1.0 == 1.0 + 2i}

Example of "round-off" error

Series[(1 - Cos[x^2])/x^2, {x, 0, 8}]

Plot[{((x^2)/2),((1-Cos[(x^2)])/(x^2))},{x,-0.0005,0.0005},PlotStyle->{{Blue},{Thickness[0.01]}},PlotPoints->100,PlotLabel->"(1 - Cos[(x^2)]/(x^2)"]

Using Numerical "Microscope"

Some functions are not automatically loaded into Mathematica. The "Needs" command tells Mathematica to load a package if a function within is called. The following package helps to illustrate computer precision.

Needs["ComputerArithmetic`"];

Ulp[1.0] (* "Unit in Last Place" *)

MachineError[Sin[x]^2, x -> 1]

Needs["ComputerArithmetic`"];ComputerArithmetic`MicroscopicErrorPlot[((Sin[x])^2),{x,1,15}]

Needs["ComputerArithmetic`"];ComputerArithmetic`MicroscopicErrorPlot[Log[Tan[((Sin[x])^2)]],{x,1,15}]

Subtraction of Machine-Precision numbers significantly reduces accuracy.

Needs["ComputerArithmetic`"];ComputerArithmetic`MicroscopePlot[x-(1/x),{x,1,15}]

{MachineError[x - 1/x, x -> 1 - Ulp[1.0]], MachineError[x - 1/x, x -> 1 ], MachineError[x - 1/x, x -> 1 + Ulp[1.0]]}

{Ulp[0.0], Ulp[1.0]}

Needs["ComputerArithmetic`"];ComputerArithmetic`MicroscopicErrorPlot[x-(1/x),{x,1,15}]

Addition is not so bad:

{MachineError[x + 1/x, x -> 1 - Ulp[1.0]], MachineError[x + 1/x, x -> 1 ], MachineError[x + 1/x, x -> 1 + Ulp[1.0]]}

Arbitrary-Precision

1`32

Precision[1.0`32]

i = 1.0`32;
While[1.0`32 + i != 1.0`32, i = i/2];i

{1.0`32 == 1.0`32 + i, 1.0`32 == 1.0`32 + 2i}

$MachinePrecision gives the Fastest Calculation

Nest[Sin, 0.1, 10^4] // Timing

Nest[Sin, 0.1`32, 10^4] // Timing

Nest[Sin, 0.1, 10^6] // Timing

MachinePrecision is about 500 times faster!

(10^6/%[[1]])/(10^4/%%[[1]])

Summation Example

Since Mathematica can perform algebraic computations, it can often find exact simplifications for expressions. However, numerical computation with Machine-Precision numbers contain errors. As more calculations are performed, these errors increase. Usually, the only way to maintain accuracy for large computations is to increase the precision of the representation of floating-point numbers.

One simple example of numerical (round-off) error is the computation of the following sum:

sumUp[n_] := Sum[(1/i), {i, 1, n}];
sumDown[n_] := Sum[(1/i), {i, n, 1, -1}];

Notice the order of the two sums, and notice that the sums are equal.

sumUp[10] == sumDown[10]

Mathematica knows the exact formula for the sum:

sumUp[n]

In order to compute these sums using Machine-Precision numbers, we will use the numerical value function, N[...].

sumUpN[n_] := Sum[N[(1/i)], {i, 1, n}];
sumDownN[n_] := Sum[N[(1/i)], {i, n, 1, -1}];

As you can see below, these are not equal!

ListPlot[Table[sumUpN[n]-sumDownN[n],{n,1,1000}],PlotStyle->PointSize[0.02]]

sumUpN[10^5] - sumDownN[10^5] // Timing

The computation used for the plot above requires a very long time as n increases. This is because the required number of calculations scale as n!. (Do you know why?) We can make things faster in the up direction by using Nest[...], but each expression working down from large n must be calculated independently.

((sumList = Table[{n, sumUpN[n]}, {n, 1, 1000}]) // Shallow) // Timing

ListPlot[sumList,PlotLabel->"HarmonicNumber[n]"]

nextUp[{n_, sum_}] := {n+1, sum +N[ (1/(n +1))]};

((sumFast = NestList[nextUp, {1,1}, 1000-1]) // Shallow) // Timing

Max[Abs[sumFast - sumList]]

Subtractive Cancellation Example

Another practical and useful example of a major source of numerical error in computational physics is subtractive cancellation (a type of "round-off" error.) When two numbers are subtracted from each other, the result is less precise as the original two numbers!

We will illustrate this by computing a first-order and second-order approximation to a derivative and compare this approximation with the "exact" value.

A first-order approximation to the derivative of the function f[x] is
	(∂_x) f ~ ((f[x + h] - f[x])/h)
It is called "first-order" because it's accuracy improves linearly with decreasing h.

A second-order approximation to the derivative of the function f[x] is
	(∂_x) f ~ ((f[x + h] - f[x-h])/(2h))
It is called "second-order" because it's accuracy improves linearly with decreasing (h^2). This is called a "central difference" because the difference of f[x] on both sides of x are evaluated.

f[x_] := x^3;

First Order Derivative

error1[h_,x_:1] := Abs[3 x^2 - (f[x+h] - f[x])/h]

err1List = Table[{10.0^h,error1[10.0^h]},{h,-14,-1,0.5}];

err1List // Shallow

ListLogLogPlot[err1List,PlotStyle->PointSize[0.02], AxesLabel -> {"Δh", "Error"}]

Often, I set Graphics options to make the plot easier to read.

myPlot[items__]:=ListLogLogPlot[items,Frame->True,GridLines->Automatic,PlotStyle->{PointSize[0.02`],RGBColor[1,0,0]},BaseStyle->{FontWeight->"Bold",FontSize->12}]

p1 = myPlot[err1List,PlotLabel->"1st Order Difference Approximation"];

Second Order Derivative

error2[h_,x_:1] := Abs[3 x^2 - (f[x+h] - f[x-h])/(2h)]

err2List = Table[{10.0^h,error2[10.0^h]},{h,-14,-1,0.5}];

p2 = myPlot[err2List,PlotLabel->"2nd Order Difference Approximation"];

Show[p2,p1,PlotLabel->"Difference Error", FrameLabel -> {"Δh", "Error"}]

Arbitrary-Precision

10.0`32^-16

err3List = Table[{10.0^h,error2[10.0`32^h]},{h,-14,-1, 0.5`32}];

err3List // Shallow

p3 = myPlot[err3List,PlotLabel->"2nd Order Arbitrary-Precision"];

Show[p3,p2,p1,PlotLabel->"Finite  Difference  Error", FrameLabel -> {"Δh", "Error"}]

Summary

Mathematica computes with three types of numbers: Integers, Reals, and Arbitrary-Precision Numbers. 
Integers are exact, but they can not represent a continuous variable. Real numbers use the floating-point processor within your computer. They are approximate numbers, but calculations can be fast. Most of your work in this course will be with Real numbers. Arbitrary-Precision numbers use a special, internal representation within Mathematica.  They have a known (and controlable) precision.
Finally, we examined the finite-difference approximation to the derivative of a function. We found that round-off error limited the accuracy of numerical approximations to derivatives.