Functions and π

Math 242 Modern Computational Math, spring 2021

Simple Functions

Here is how to code the function f(x)=(x^2) in Mathematica.

f[x_]:=x^2

Note the use of delayed assignment. You almost always want to use delayed assignment for functions.
Also note the underscore following the x inside of the square brackets: this tells Mathematica that the x is a placeholder for anything value that you put into the function.
Let’s use our function to make a list of squares:

squares=Table[f[n],{n,3,20}]

We can also plot the function like this:

Plot[f[x],{x,-1,1}]

We can also define functions with more than one variable. Here is g(x,n)=(x^n).

g[x_,n_]:=x^n

To try this out, let’s make a table (that is, a list) of functions:

functions=Table[g[x,k],{k,1,10}]

Now plot our list of functions:

Plot[functions,{x,0,1},PlotRange->{{0,1},{0,1}}]

More Complicated Functions (that is, Modules)

What if we want to make a more complicated function that encapsulates multiple lines of code? To do this in Mathematica, we use a concept called a module.

Here is the basic syntax for a module:

Module[{local variables}, expressions]

Local variables are variables that exist only within the module. A local variable in a module may have the same name as another variable outside the module, but it is a different variable. Expressions are the statements that run when you call the module. Multiple expressions in a module must be separated by semicolons.

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 Intro Mathematica assignment. Here we are revisiting the problem as an example of how to build more advanced functions.

A first module

computeSum[n_]:=Module[{sum=0},         (* sum is a local variable *)
Do[sum += 1/i^2,{i,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 module that does almost the same thing

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

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

computeSumFP[30]

Comparing our two modules: accuracy and runtime

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

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

N[results,25]

Are these two results the same? They agree to six digits, but they 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[100],20]

Looks like we only have 1 correct digit to the right of the decimal point.

We can use Table to compute many approximations at once:

Table[ N[approxPi[k],20], {k,100,1000,100}]

Even with 1000 terms of the sum, we only get 2 correct digits to the right of the decimal point.

Determining the number of correct digits

Here is a module, digitsOfAgreement, that takes two numbers and returns the number of digits in which they agree.

digitsOfAgreement[a_,b_]:=Module[{i=0},
While[RealDigits[a,10,1,-i][[1,1]] == RealDigits[b,10,1,-i][[1,1]],
i++];
i (* return value: the number of correct digits after the decimal point *)
]

We can ignore the internal details of this module for now, but here is a simple demo:

digitsOfAgreement[1.23456,1.23056]

We can use digitsOfAgreement to find the number of correct digits in an approximation of π. Here is an example:

approx = approxPi[200];
N[approx]
digitsOfAgreement[approx,Pi]

How does digitsOfAgreement work?

The module digitsOfAgreement relies on Mathematica’s RealDigits function, which converts a number to a list of digits. More specifically, it returns a pair {digitList, decPt}, where digitList is the list of digits and decPt is the number of digits to the left of the decimal point. For example:

RealDigits[Pi,10,20]

If we simply want the list of digits, use double brackets (the Part operator) to extract the list of digits:

RealDigits[Pi,10,20][[1]]

If we want just a single digit, we can add a fourth argument to RealDigits specifying the place value at which the digit list should start. Digits to the right of the decimal point are specified by negative place values. For example, -2 gives the hundredths place — that is, the (10^(-2)) place. For example, here is how we can get the hundredths-place digit of π:

RealDigits[Pi, 10, 1, -2][[1]]

Note the curly brackets in the output above. The previous result is a list containing 4, not just the number 4. We can further use the Part operator to extract the digit 4 from its list, like this:

RealDigits[Pi, 10, 1, -2][[1,1]]

Now suppose we want to compare digits at the same place value in two numbers a and b. We can use an equality comparison like this:

a=3.1415;
b=3.13;
RealDigits[a,10,1,-2][[1,1]] == RealDigits[b,10,1,-2][[1,1]]

In digitsOfAgreement, the equality comparison is placed in a While loop to compare digits until we find digits that don’t agree. For each pair of digits that agree, we increment a counter, i, that starts at zero. When we find digits that don’t agree, we stop looping, and the value of i is the number of digits of agreement.

a=3.1415;
b=3.13;
i=0;
While[RealDigits[a,10,1,-i][[1,1]] == RealDigits[b,10,1,-i][[1,1]],i++];
i

Packaging this code in a module, we have the digitsOfAgreement module presented above.