• 5.3. More loop examples
  • 5.3.1. Example 1: Riemann sums
  • 5.3.2. Example 2: Pascal's triangle
  • 5.3.3. Example 3: Periodic extensions
  • 5.3.4. Example 4: Drawing graphs
  • 5.3.5. Example 5: Butterfly curve
  • 5.3.6. Example 6: Animations

Programming in Maple

Roger Kraft
Department of Mathematics, Computer Science, and Statistics
Purdue University Calumet

roger@calumet.purdue.edu

5.3. More loop examples

In this section we work on more involved examples using loops.

>

5.3.1. Example 1: Riemann sums

A common use for a for-loop is to compute a sum. In this example we show how to use a loop to compute Riemann sums from calculus. But first let us look at a couple of simple sums written as for-loops.

A well known result about sums of integers is that the sum of the first  positive integers is /2. Let us verify this for  using a for-loop. We want a for-loop that will compute the following sum, written in sigma notation.

The basic idea of using a for-loop to compute a sum is that we compute a running total. If you want to add 1+2+3+4+5+6+7, we start with 0, then we add 1 to 0, then we take that result and add 2 to it, then we take that result and add 3 to it, then we take that result and add 4 to it, etc. In the following for-loop, we will let the running total be called s  (for s um) and we will initialize s  with the value 0. Notice how the single command in the body of the for-loop adds another integer to the sum for each iterate of the loop. Compare this for-loop with the sigma notation above. In what ways are they alike and in what ways do they differ?

>

s := 0:

>	for i from 1 to 1000 do s := s + i od:

Notice that we have a colon at the end of the loop so that we do not see 1000 lines of output. Let us check the value of s , the sum.

>

s;

Now verify this result against the formula /2.

>	1000*(1000+1)/2;

By the way, Maple can verify the above (symbolic) formula for us by using the sum  command.

>	i:='i': n:='n':

>	sum( i, i=1..n );

>	factor( % );

>

Exercise : Write a procedure called add_list  that takes one input, a list of numbers, and computes the sum of the numbers in the list. Do not use Maple's add  or sum  commands.

>

Here is another example of computing a sum. We know from calculus that the number  is the sum of the reciprocals of all the factorials. So we can compute an approximation of  by summing the reciprocals of the factorials from say 0! to 10!. So we want a for-loop that will compute the following sum, written in sigma notation.

Compare the next for-loop with this sigma notation.

>

e := 0;

>	for n from 0 to 10 do e := e + 1/(n!) od;

That is not what we really want. Let us modify the loop so that it computes with decimal numbers.

>

e := 0;

>	for i from 0 to 10 do e := e + 1.0/(i!) od;

In this example we put a semi colon at the end of the loop so that we could see the answers converge to the correct value of , which is given by the next command (to ten decimal places).

>	evalf( exp(1) );

Now let us turn to Riemann sums. The following execution group will compute a (left hand) Riemann sum for the function f  over the interval from a  to b  using n  rectangles. If you want to, you can change any of the values for f , a , b  or n .

>	f := x -> sin(x);      # The function.

>	a := 0;                # Left hand endpoint.

>	b := 3/2*Pi;           # Right hand endpoint.

>	n := 75;               # How many rectangles.

>	L||n := 0:             # Set the running total to 0.

>	for i from 0 to n-1 do

>	x||i := a+i*(b-a)/n;          # Compute a partition point.

>	L||n := L||n+f(x||i)*(b-a)/n: # Add area of rectangle to sum.

>

od:

>	'L[n]' = %;            # Display the result symbolically,

>	'L[n]' = evalf(%%);    # and numerically

You should think of the for-loop in this execution group as the Maple version of the following formula from calculus.

                                                              

This sigma-notation represents the sum and the for-loop in the execution group computes the sum. An interesting feature of actually computing the sum, instead of just representing it abstractly, is that the computation must be very precise about what it is doing. For example, notice how  in the sigma-notation is a short hand for the Maple expression a+i*(b-a)/n . In calculus courses where this sigma-notation is used, students often do not realize just what the   represents. But if you have to write a Maple for-loop to implement the sum, then you are forced to explicitly state what the symbol  means. Maple will not be able to figure out for you what  represents if you use it in a for-loop without defining it.

>

Exercise : Modify the execution group for computing  to compute , right hand sums. Then try  nd , the trapezoid sums and Simpson's rule.

>

5.3.2. Example 2: Pascal's triangle

Our next example shows how we can develop a procedure for printing out Pascal's triangle. The next for-loop prints out a table of expansions of .

>	x:='x': y:='y':

>	for i from 0 to 8 do

>	expand( (x+y)^i );

>

od;

The coefficients in these polynomials have a name (the binomial coefficients) and they make up what is called Pascal's triangle. Pascal's triangle and the binomial coefficients have a lot of interesting properties, so let us see if we can extract Pascal's triangle out of the last example.

The Maple command coeffs  returns the coefficients of a polynomial as an expression sequence. The coefficients of the above polynomials make up Pascal's triangle. So let us see if the next for-loop will give us Pascal's triangle.

>	for i from 0 to 8 do

>	coeffs( expand( (x+y)^i ) );

>

od;

That did not work. The output from coeffs  gives the coefficients of each polynomial as an expression sequence, but the coefficients are not in the same order as in the polynomial (see the next two commands).  

>	expand( (x+y)^8 );

>	coeffs( % );

So we do not yet have Pascal's triangle. Let us try another approach. Since polynomials are data structures, let us use our knowledge of Maple's data structure commands to get the coefficients of each polynomial in the order we want them. Since we want the coefficients to be in an expression sequence, we will use the seq  command. We can use an op  command to pick off the individual terms of the polynomial, and then we can use the coeffs  command to return the coefficient of just one term at a time.

>	expand( (x+y)^8 );

>	seq( coeffs(op(j, %)), j=1..nops(%) );

That worked. Here is a for-loop built around the last two commands.

>	for i from 0 to 8 do

>	expand( (x+y)^i );

>	seq( coeffs(op(j, %)), j=1..nops(%) );

>

od;

Oops, too much output. Put a colon at the end of the final od  and use a print  command in the body of the loop.

>	for i from 0 to 8 do

>	expand( (x+y)^i );

>	seq( coeffs(op(j, %)), j=1..nops(%) );

>	print( % );

>

od:

There we go. Notice that every number in Pascal's triangle is the sum of the two numbers just above it.

Of course, if binomial coefficients are important in mathematics, then we should expect Maple to have a command for computing them directly instead of picking them out of the expansion of .  The command binomial(n,j)  returns the coefficient of the j'th term of .

>	seq( binomial(8,j), j=0..8 );

So we can get Pascal's triangle with the following for-loop.

>	for i from 0 to 8 do

>	seq( binomial(i,j), j=0..i );

>

od;

Notice how this loop uses two index variables, one for the loop itself and one for the seq  command. The i  index variable is counting the rows of our output and the j  index variable is counting the "columns". And the i  index variable from the "outer loop" is the final value for the j  index in the "inner loop".

>

Here is a procedure that allows us to conveniently print out as many lines of Pascal's triangle as we want.  (Notice the multiple levels of indentation to make this easier to read.)

>	pascal := proc(m, n)

>	local i, j, x, y;

>	for i from m to n do

>	expand( (x+y)^i );

>	seq( coeffs(op(j, %)), j=1..nops(%) );

>	print( % );

>

od;

>

end;

Let us try it out. (On my computer screen I can fit up to 18 lines of Pascal's triangle.)

>	pascal(0,8);

Notice how this procedure uses the print  command to print out one line of Pascal's triangle for each loop through the for-loop. Normally a procedure only has one output, which is the last line executed by the procedure (the "return value" of the procedure). But this procedure has n  lines of output. These extra lines of output have a name, they are called side effects.  Anything else that a procedure does besides returning its "return value" is called a side effect. (You may wonder just what the return value of this procedure is. We will look at that question later in this chapter.)

>

Exercise : Modify the procedure so that it uses the binomial  command instead of the expand  and coeff  commands.

>

5.3.3. Example 3: Periodic extensions

Our third example will use a while-loop. We show how to take an arbitrary function g defined on an interval between two numbers  and , , and produce a function f that is periodic on the whole real line, with period , and is equal to the original function  between  and . This new function is called the periodic extension  of the original function.

>

Before showing how to use Maple to define the periodic extension, let us try to describe it in words. We start with what we might think of as a segment of a function g, defined just between  and . The periodic extension f will have copies of this segment repeated over and over again on the real line. There will be one copy between  and  (where  is the length of the interval that g is originally defined on) and another copy between  and , etc. How should we define the extension f on the "first" interval from  to ? Visually, we would just imagine sliding the graph of g from the interval  over onto the interval . We would express this mathematically by saying that, if  is in , then  (make sure you understand this step). For the interval from  to , we would again just imagine sliding the graph of g from  over to . And if  is in , then we would let . In general, for any positive integer , if  is in , then . In short, to define  for  we need to keep subtracting  from  until we get a number that is between  and , and then use that number to evaluate g. This process of subtracting 's from  until we get a number between  and  is exactly what a while-loop can do for us. Since we do not know the value of  ahead of time, we do not know how many 's we need to subtract, but while-loops are exactly what we should use when we want to iterate something an unknown number of times.

Exercise : Figure out how we would define  for .

>

Let g  be a Maple function and let a  and b  be two numbers with a<b . Here is how we define the periodic extension f  of g .

>	f := proc(x)

>

local y;

>

y := x;

>	while y >= b do y := y-(b-a) od;

>	while y < a  do y := y+(b-a) od;

>

g(y);

>

end;

This definition of f  works for any g , a  and b . Let us define specific g , a , and b .

>	g := x -> x^2;

>	a := -2; b:= 3;

Now graph the periodic extension of g .

>	plot( f, -7..18, scaling=constrained, discont=true, color=red );

We should note here that if we try to graph f  as an expression, then we get an error message.

>	plot( f(x), x=-7..18, scaling=constrained, discont=true, color=red );

Error, (in f) cannot determine if this expression is true or false: 3-x <= 0

Similarly, if we try to look at the formula for f  as an expression, we get the same error message.

>

f(x);

Error, (in f) cannot determine if this expression is true or false: 3-x <= 0

We will find out how to fix this later in this chapter. For now, here is a way to graph f  as an expression if you really want to.

>	plot('f(x)', x=-7..18, scaling=constrained, discont=true, color=red);

>

Exercise : With the same g  as in the last example, try changing the values of a  and b . How does this change the periodic extension?

>

Exercise : Try defining periodic extensions for some other functions g .

>

Exercise : How are we defining f(a+k*p)  for any integer k ? What other reasonable choices are there for the value of f(a+k*p) ?

>

Recall that a function f is even  if for all  it is true that . A function is even if its graph is symmetric with respect to the -axis. Recall that a function is odd  if for all  it is true that . A function is odd if its graph is symmetric with respect to the origin, meaning that if we rotate the graph of f by 180 degrees, then the graph remains unchanged.

If  is a positive number and we define the periodic extension of a function g on the interval , then the periodic extension will be an odd function for some g, an even function for some other g, and neither even nor odd for most g. For example, the periodic extension of  on the interval  is an even function.

>	g := x -> sin(x);

>	a:=0; b:=Pi;

>	plot( f, -3*Pi..3*Pi );

And the periodic extension of  on the interval  is an odd function.

>	g := x -> cos(x);

>	a:=0; b:=Pi;

>	plot( f, -2*Pi..2*Pi, discont=true, color=red );

There is a way to define a periodic extension so that it is always even or always odd, no matter what the function g on the interval  looks like. The following procedure defines the even periodic extension of a function g on the interval from 0 to b. The even extension defined by this procedure has period . (Notice the use of an anonymous function in the last line of the procedure.)

>	f_even := proc(x)

>

local y;

>

y := x;

>	while y >= b do y := y-2*b od;

>	while y < -b do y := y+2*b od;

>	( z -> piecewise(z<0, g(-z), g(z)) )(y);

>

end;

The next procedure defines the odd,  periodic extension of a function g on the interval from 0 to .

>	f_odd := proc(x)

>

local y;

>

y := x;

>	while y >= b do y := y-2*b od;

>	while y < -b do y := y+2*b od;

>	( z -> piecewise(z<0, -g(-z), g(z)) )(y);

>

end;

Let us try some examples using these two procedures. Here are the even and odd  periodic extensions of  on the interval .

>

g := sin;

>

b := Pi;

>	plot( f_even, -3*Pi..3*Pi );

>	plot( f_odd,  -3*Pi..3*Pi );

Here are the even and odd  periodic extensions of  on the interval .

>

g := cos;

>

b := Pi;

>	plot( f_even, -3*Pi..3*Pi );

>	plot( f_odd,  -3*Pi..3*Pi );

Let us try the identity function,  on the interval .

>	g := x -> x;

>

b := 1;

>	plot( f_even, -3..3 );

>	plot( f_odd,  -3..3 );

>	plot( f, -3..3 );

Let us try a non symmetric piece of a parabola on the interval .

>	g := x -> (x-1/4)^2;

>

b := 1;

>	plot( f_even, -3..3 );

>	plot( f_odd,  -3..3, discont=true, color=red );

>	plot( f, -3..3, discont=true, color=red );

>

Exercise : Try defining even and odd periodic extensions for some other functions g .

>

Exercise : Part (a) Under what conditions on the function g  will the function f_odd  be continuous at 0 and b ? Under what conditions on the function g  will the function f_even  be continuous at 0 and b ?

>

Part (b) Under what conditions on the function g  will the functions f_odd  and f  be the same function (where f  is the b  periodic extension of g  on the interval from 0 to b ). Under what conditions on the function g  will the functions f_even  and f  be the same function.

>

Here is an interesting example that uses a parabola.

>	g := x -> 4*x*(1-x);

>

b := 1;

>	plot( f_odd, -4..4 );

Let us compare this periodic function to a sine curve.

>	plot( [ f_odd, x->sin(Pi*x) ], -2..2 );

Notice how amazingly similar the two functions are. The following parametric graph uses the odd periodic extension of g . Which is the real circle?

>	plot( [ [f_odd, t->f_odd(t-1/2), 0..2],

>	[cos, sin, 0..2*Pi] ],

>	scaling=constrained );

>

Exercise : Given a function g defined on an interval , write a procedure f_oddeven  that defines a  periodic extension of g that is odd and such that the horizontal shift of the extension by the amount  is even. Similarly, write a procedure f_evenodd  that defines a  periodic extension of g that is even and such that the horizontal shift of the extension by the amount  is odd. (Hint: Think of how sin  is built up from just one quarter of sin 's graph, the part from 0  to Pi/2 , and then think of how cos  is built up from just one quarter of cos 's graph.)

>

Exercise : Make up your own periodic function by extending some function g , and then use your periodic function in place of the trig functions sin and cos in several parametric curves and surfaces like cardiods, spirals, lemniscates, roses, sphere, torus, etc. Try to come up with a really unusual curve or surface.

>

>

5.3.4. Example 4: Drawing graphs

Our fourth example uses a for-loop and a list data structure to draw a graph. Each iterate of the for-loop will compute the coordinates of a point and put the point in a list of points that are to be plotted. Then the plot  command will be used to draw a graph of all the points in the list. This is a fairly common way to create a graph in Maple and this is fundamentally how all graphs are drawn in Maple.

The following for-loop computes 9 equally spaced sample points in the interval from 0 to . Each iterate of the loop computes one sample point, evaluates the sin  function at the sample point, and puts the ordered pair of the sample point and the value of sin  into a list called data . Then the plot  command is used to graph the list data , and we get a graph of the sin  function. Notice that data  starts off as the empty list, [] , and as each new point [x,sin(x)]  is computed, it is put into the list using a command of the form data:=[op(data),[x,sin(x)]] . Notice how this is very much like computing a running sum s , where we start off with s  equal to the "empty" sum, 0, and then as each new term in the sum is computed it is added to the sum using a command like s:=s+something .

>

N := 8;

>	b := 2*Pi;

>	data := [];

>	for n from 0 to N do

>	data := [ op(data), [n*b/N, sin(n*b/N)] ];

>

od;

Notice that the output from the for-loop is in symbolic form. Since we are going to just plot the contents of the list, we do not need exact symbolic results. So let us modify the for-loop so that it computes with decimal numbers. When working with loops that might create thousands of points, this could help speed things up. (If we want, we could even use the evalhf  command and hardware floating points for even more speed.) First, we need to reinitialize list  to be the empty list again.

>	data := [];

>	for n from 0 to N do

>	data := [ op(data), [evalf(n*b/N), sin(evalf(n*b/N))] ];

>

od;

Now use plot  to graph the list data .

>	plot( data );

>

Exercise : Modify the above example so that it can be used to graph the sin  function over any interval from a  to b . (You should suppress the output from the for-loop in your example, so that you do not fill up the worksheet with an unwieldy amount of output.)

>

Exercise : Modify the above example so that it is easy to change the function that is being graphed. Try graphing the function  over the interval .

>

Exercise : Consider the following three execution groups. Each one graphs the sin  function and each one uses a slight variation on the last for-loop. Suppose N  is a very large number. How do you think these three execution groups would compare, in terms of speed, with the last for-loop and with each other?

>

N := 8;

>	b := 2*Pi;

>	data := []:

>	for n from 0 to N do

>	data := [ op(data), [evalf(n*b/N), evalf(sin(n*b/N))] ];

>

od:

>	plot( data );

>

>

N := 8;

>	b := 2*Pi;

>	data := []:

>

n := 'n':

>	x := evalf( n*b/N );

>	for n from 0 to N do

>	data := [ op(data), [x, sin(x)] ];

>

od:

>	plot( data );

>

>

N := 8;

>	b := 2*Pi;

>	data := []:

>

n := 'n':

>	x := evalf( b/N );

>	for n from 0 to N do

>

p := n*x;

>	data := [ op(data), [p, sin(p)] ];

>

od:

>	plot( data );

>

The following for-loop computes five equally space points on the circumference of the unit circle and puts these point in a list called data . Notice that once again data  starts off as the empty list, [] , and as each new point [x,y]  is computed it is put into the list using a command of the form data:=[op(data),[x,y]] .

>

N := 5;

>	data := [];

>	for n from 0 to N do

>	data := [ op(data), [cos(evalf(n*2*Pi/N)), sin(evalf(n*2*Pi/N))] ];

>

od;

Notice how the list of points grew by one point with each iteration of the loop. Now use the plot  command to plot the list.

>	plot( data, style=point, scaling=constrained );

Now plot it with lines connecting the dots.

>	plot( data, style=line, scaling=constrained );

In the next for-loop we use the names x  and y  to represent the calculation of the coordinates of each of the points. Convince yourself that the next execution group does exactly the same calculation as the previous execution group. The reason for rewriting it like this is to try and make the execution group easier to read. (And this time we are suppressing the output from the for-loop.)

>

N := 5;

>	data := [];

>	n := 'n';  # be sure that n is unassigned

>	theta := evalf( n*2*Pi/N );

>	x := cos(theta);

>	y := sin(theta);

>	for n from 0 to N do

>	data := [ op(data), [x, y] ];

>

od:

Plot the points again, to verify that it was the same calculation.

>	plot( data, style=line, scaling=constrained );

Now go back to the last execution group and change the value of N  from 5 to 7 or 9. The execution group will plot 7 or 9 equally spaced points on the circumference of the unit circle. Try any positive integer for N .

>

Here is a variation on this example. The next execution group plots 5 equally spaced points on the circumference of the unit circle, but they are not computed in sequential order around the circle. Try to figure out exactly what this version computes and how it does it.

>

N := 5;

>

J := 2;

>	data := []:

>

n := 'n':

>	theta := J*2*Pi/N;

>	x := cos(n*theta);

>	y := sin(n*theta);

>	for n from 0 to N do

>	data := [ op(data), [x, y] ];

>

od:

>	plot( data, style=line, scaling=constrained );

Try changing N  to 7 and J  to 3. Try several different values for N  and J .

>

Exercise : Convert the last execution group into a procedure that takes two positive integer parameters, the N  and J , and draws the appropriate graph. This will make it a lot easier to try out different parameter values.

>

One more variation on this example. This version has three parameters, N , J , and K . Try to figure out exactly what this execution group is doing. (Hint: It does pretty much the same thing as the previous version, but it does it twice per iterate of the loop.)

>	N, J, K := 36, 21, 9;

>	data := []:

>

n := 'n':

>	theta1 := J*2*Pi/N;

>	theta2 := K*2*Pi/N;

>	x1, y1 := cos(n*(theta1+theta2)),

>	sin(n*(theta1+theta2));

>	x2, y2 := cos(n*(theta1+theta2)+theta1),

>	sin(n*(theta1+theta2)+theta1);

>	for n from 0 to N do

>	data := [ op(data), [x1, y1], [x2, y2] ];

>

od:

>	plot( data, style=line, scaling=constrained, axes=none );

>

Here are some values for the parameters N , J  and K  that produce nice graphs.

15, 8, 13

28,19,15

39,33,27

19,13,11

>

Exercise:  Part (a) Convert the last execution group into a procedure that has three positive integer parameters and draws the appropriate graph. Call your procedure with a number of different parameters.

>

Part (b): Convert your procedure from Part(a) into a procedure that takes no input parameters and generates the three integers it needs randomly. Your procedure should use the Maple function rand  to generate the random integers (the expression rand(a..b)() , with both  sets of parentheses,  generates a random integer between a  and b ). Have the procedure print out the values of the randomly chosen integers and then draw the appropriate graph. Run this procedure many times. You should get some very elegant graphs.

>

>

5.3.5. Example 5: Butterfly curve

Our fifth example uses a for-loop and a data structure (a list) to draw a fairly complex graph, a version of a butterfly curve.

The curve that we draw is not in fact a curve. Instead we will be plotting points, thousands of points, and not connecting them together with line segments. The points we will be plotting will be computed by a for-loop and placed in a list (as in the last example). Then the plot  command will be used to draw a graph of all the points in the list (without connecting them together with line segments). Almost all of the work in this example is in computing the coordinates of all the points that we want to plot.

The loop in this example is fairly computationally intensive. You should save all of your work before executing it, in case it takes too long for it to compute on your computer (do not try this on a Pentium, use at least a Pentium II).  If your computer is really fast, you can try changing N . Try  N:=21000 , or N:=41900  (which should run for a pretty long time). Different values of N  give different butterfly curves.

>	r := phi -> ( exp(cos(phi))-2*cos(4*phi) )*sin(99999999*phi)^4;

>	N := 11500;              # Number of points to compute.

>	h := evalf(2*Pi/N);      # Step size between points.

>	n := 'n';                # Just to be safe.

>	x := r(n*h)*sin(n*h);    # x-coord of a point on the curve

>	y := r(n*h)*cos(n*h);    # y-coord of a point on the curve

>	data := [];              # Start with an empty list.

>	for n from 1 to N do     # This do-loop computes the butterfly.

>	data := [ op(data), [evalhf( x ), evalhf( y )] ]

>

od:

Now that we have computed our list of points, let us graph it.

>	plot( data, style=point, symbol=point, color=black );

The original reference for butterfly curves is The Butterfly Curve , by Temple H. Fay, in The American Mathematical Monthly, Volume 96, Issue 5 (May, 1989), pages 442-443. The version of the butterfly curve in this example is from A Study in Step Size , by Temple H. Fay, in Mathematics Magazine, Volume 70, No. 2, April 1997, pages 116-117.

>

>

5.3.6. Example 6: Animations

Our last example is the use of a for-loop to compute the frames of an animation. Maple's animate  command is limited in the kind of animations it can draw. For example, it cannot animate the graphs of equations. There is another way to create animations in which we use a for-loop to create and label a sequence of graphs (the frames) and then we use a special form of the display  command to turn the sequence of graphs into an animation. Here is an example of this technique.

First, a for-loop is used to create and name, using concatenated names, 51 two-dimensional graphs, each with a parameter slightly changed. Then the display  command is used to "sequence" the 51 graphs into a short movie. To view the movie, after the first graph is visible, place the cursor on the graph. Some VCR type buttons will appear at the top of the Maple window. Click on the "play" button. (The for-loop takes a little while to complete.)

>	x := 'x': y := 'y':

>	for i from -20 to 30 do

>	p||i := plots[implicitplot]( x^3+y^3-5*x*y = 1-i/8,

>	x=-3..3, y=-3..3, numpoints=800, tickmarks=[2,2] )

>

od:

>	plots[display]( p||(-20..30), insequence=true );

>

>

>