This will be a two-week unit.
(A) Random Number Generators.
In Maple we can generate a 12-digit non-negative random integer using the command
> rand();
To check if the 12-digit random number is prime or not we could do the following:
> n := rand();
> isprime(n);
In fact, it would be interesting to compute the percentage of all 12-digit integers that are actually prime. (This would eventually lead us to the famous Prime Number Theorem, a very deep result in classical Number Theory.) Later we will explore a Monte Carlo strategy to estimate this percentage.
More modestly, suppose that we wish to simulate the throw of a fair die (having faces numbered 1 through 6). One way to achieve this is by defining a Maple procedure (that we will call die) as follows:
> die := rand(1..6):
> die();
> die();
> die();
Notice that each time you enter die(), another random integer between 1 and 6 is "selected."
Suppose that we wish to simulate the throwing of a pair of dice and the resulting sum of numbers on the two faces.
> die1 := rand(1..6):
> die2 := rand(1..6):
> sumOfFaces := die1() + die2();
Of course, to simulate anything interesting we would need to generate
a large number of random numbers. Thus we will explore loops
and branching statements next.
(B) Branching.
The if-then or if-then-else statements allow us to make decisions in Maple. For example, consider the following:
> a := 3;
> b := 5;
> if (a > b) then a
>
else b
> end if;
(Observe the syntax. Indentation is provided only to improve readability.
Maple is oblivious to indentation.)
Notice that since 3>5 is false,
the value of b is printed. Change a
to a number larger than 5 and see what happens.
Here is an example involving prime numbers:
> number := rand(); # we generate
a 12 digit random integer.
> if (isprime(number)) then print
(`The number is prime.`)
>
else print (`The number is composite.`)
> end if;
Repeat this code several times. How many prime numbers do you obtain?
Relational and Logical operators: The expressions a>b and isprime(number) in the examples above are called Boolean expressions because they evaluate to either true or false. The familiar arithmetic relations <, <=, >, >=, =, and <> are called relational operators. (The last symbol <> means not equal.) Sometimes we may wish to combine two or more Boolean expressions:
a>b and b<cThe operators and, or, not are called logical operators.
a=5 or a=3
not (a >= 4)
(a=b) and ((c=d) or (a<>d))
(C) Loops.
To simulate the repetition of an experiment we need a loop command. In Maple the loop works as follows:
> for i from 6 to 25 do
>
print(i)
> end do;
Next let us add all the integers between 1 and 1000 using a loop construct:
> sumOfNumbers := 0; #
Here we initialize the value of our counter to 0.
> for i from 1 to 1000 do
>
sumOfNumbers := sumOfNumbers + i # Here we increment
the value of the counter by the number i.
> end do:
> print (` The sum of the integers between
1 and 1000 is:`);
> sumOfNumbers;
(Again, observe the syntax. The spacing is only to aid the human reader.) What happens if the colon after the end do is replaced by a semi-colon?
(D) Generating sequences of random numbers and Monte Carlo methods.
Let us imagine that we wish to compute the area of a quarter circle x2 + y2 < 1 in the first quadrant. (Let us assume that x and y are measured in cm.) Of course we know in advance the answer, but let us employ a Monte Carlo method. We will generate a sequence of N random points in the unit square, 0 < x < 1 and 0 < y < 1, and calculate how many of them lie inside our quarter circle. We then use this to estimate the area of the quarter circle: For example, if 34 points of the 50 in the square lie inside the circle, then 68% of the random points land inside the quarter circle. Thus we would estimate the area of the quarter circle to be 0.68 times the area of the square (which of course is 1), or 0.68 sq cm.
To generate a sequence of random numbers we will use the stats library in Maple.
> with(stats[random]):
# invoking the stats library.
> N:= 50;
# here we choose a relatively small number of points.
> xrand:=uniform[0,1](N): #
we generate a sequence of 50 random numbers between 0 and 1 that will represent
x-coordinates of our 50 points.
> yrand:=uniform[0,1](N): #
we generate a sequence of 50 random numbers between 0 and 1 that will represent
y-coordinates of our 50 points.
> with(plots): # invoking
the plots library.
> points:=[seq([xrand[i], yrand[i]], i = 1..N)]:
# Here we define a collection of 50 random points.
> f := sqrt(1 - x^2);
> fGraph := plot(f, x = 0..1):
> pointGraph := plot(points, style = POINT):
> display({fGraph, pointGraph});
Now let's count the number of points that land inside the circle:
> count:=0;
> for i from 1 to N
> do
>
if yrand[i] < subs(x=xrand[i],sqrt(1-x^2)) then count:=count+1
>
end if
end do;
> count;
> print (`The fraction of points
that landed inside the quarter circle is:`);
> evalf(count / N);
> print (`Thus the area of the quarter
circle is approximately equal to:`);
> evalf( count/N * 1);
# Note that 1 represents the area of the square.
Exercise: Change the number of points to 500 and see what happens.
(E) Distance function.
In some of our geometry problems (circles, spheres, ellipses, etc.) it is efficient to employ Maple's Euclidean distance function. For example:
> with(student):
# here we invoke the student library.
> distance ([3,4], [7,11]);
# Two-dimensional distance
> distance ([1, 3, 4.5], [19.1,
3, 5.55]); # Three-dimensional distance
(F) More examples.
(1) The area of the annulus: 4 < x2 + y2 < 9. Hint: A point Q is in the annulus if and only if the distance from Q to (0,0) is between 2 and 3.
(2) The volume of the unit sphere: x2 +y2 + z2 < 1. Hint: A point P is in the unit sphere if and only if its distance to (0,0,0) is less than 1.
(3) Throw a pair of fair dice 1000 times. Find the percentage of times that the outcome is: "seven or eleven".
(4) The area of the crescent formed by the region outside the circle centered at the origin of radius 2 and inside the circle centered at (0,1) with radius sqrt(3) .