Functions and π

Math 242: Modern Computational Math

More Functions

As in other programming languages, functions can combine all sorts of loops, control structures, etc.
To do this in Mathematica, we use a concept called a Module.

As an example, let’s add up the first n values of the inverse squares of positive integers. Note that we already did this on the first assignment. Here we are revisiting the problem as an example of how to build more advanced functions.

A first attempt

computeSum[n_]:=Module[{sum=0},         (* sum is a local variable *)
Do[sum += 1/k^2,{k,1,n}];                 (* here is our loop *)
sum                                                                            (* return the sum *)
]

Try it out:

computeSum[1]
computeSum[2]
computeSum[30]

Maybe we want a decimal answer, with (say) 20 digits:

N[computeSum[30],20]

Another attempt

Here is another way. How is this different than what we did before?

computeSumFP[n_]:=Module[{sum=0.0,k},
Do[sum += 1/k^2,{k,1,n}];
sum
]

computeSumFP[30]

Comparing our two functions: accuracy and runtime

Let’s compare the output of our two functions when we ask for 100 terms of the sums. Run the following:

results={computeSum[100],computeSumFP[100]}

N[results,25]

Do these functions compute the same thing? They agree to six digits, which seems OK. However, it’s important to note that the two answers are not the same.

How fast are are these calculations? Use the Timing function to find out. For this, we will ask for a lot of terms to be computed:

Timing[N[computeSum[100000],20]]

The Timing function returns two numbers. The first number is the runtime, the second is the sum.

How about the other function? Let’s see:

Timing[N[computeSumFP[100000],20]]

That was a lot faster! The lesson is that computers are much faster when working with floating-point (i.e., fixed-precision) numbers than with arbitrary-precision (in this case, rational) numbers.

Where is π?

It is know that the sum of inverse squares converges to ((π^2)/6). Thus, we can use our sum to approximate π. Let’s define a function for this:

approxPi[n_]:=Sqrt[ computeSum[n]*6]

Try it out:

N[approxPi[1000],20]

Looks like we only have 2 digits of agreement to the right of the decimal point.

Here is one way of having Mathematica tell us how many digits of accuracy we have. First, we will compute our approximation using 100,000 terms:

approx = N[approxPi[100000],50]

Now ask Mathematica for π to 50 decimal digits:

exact = N[Pi,50]

The absolute value of the difference gives us the error:

error = Abs[approx - exact]

Now we can use -(Subscript[log, 10])(error) to “count” the number of leading zeros of the error (plus/minus one). The Floor function rounds down the the closest integer.

Floor[-Log[10,error]]