Primer.html

MATH 371 -- Numerical Analysis

A MAPLE7 PROGRAMMING PRIMER

John Little

September 6, 2001

Introduction

In addition to the built-in graphical and symbolic computational functions

we have discussed so far, Maple 7 also contains a built-in procedural programming

language (similar in some ways to the Pascal language , in case you have seen that).

This makes it possible to "automate" computations and create new functions to extend

the capabilities of the basic system.

In this presentation we will give an overview of this language including:

1) Basic control structures (loops, selections)

2) Maple data structures

3) Structure of Maple procedures

4) An extended example (the Bisect procedure)

Basic Control Structures

Repetition structures (loops)

for

As a part of big jobs, we often want to do a command or sequence of commands once for each

value of a " counter" variable or once for each item in a list. Hence Maple contains a

for loop structure. There are two basic forms, one very similar to the Pascal count-controlled

for loop, one slightly different.

count-controlled loops

Here the completely general syntax is

for <counter> from <start> by <change> to <finish> do

end do;

The indentation and formatting here is not important, but it

makes the code more readable, I think!

The <statement sequence> is any group of Maple commands.

The end do; line at the end marks the end of this sequence of commands.

What this will do is to set the counter equal to <start> to begin

with, perform the <statement sequence> once, then change the

value of the counter by adding <change> and repeat the

statement sequence. The repetitions continue until counter

reaches the value <finish>.

The <start> and/or <change> can be omitted if they are both 1.

For example--

> s := 0:

> for i to 10 do s := s + i^2 end do:

> s;

>

> Alist := [1,4,9,16,25,36]:

> for i from 2 by 2 to 6 do print(Alist[i]) end do;

list-controlled loops

The for loop can also be used to execute one or more statements for each item in

a list or a set. For more information on these data structures, see the next section.

The syntax here is

For <element> in <list or set> do

od;

> Alist := [3,4,5,6,7]:

> for x in Alist do print(x) end do;

> Aset := {1,1,2,3,5,8,11}:

> for x in Aset do isprime(x) end do;

while

The syntax of the while loop is very similar to the Pascal while loop:

while <condition> do

od;

The line breaks and indentation are irrelevant but a format like this sometimes improves

readability. The "od" matches the "do" at the start of the statement sequence

and marks the end of the body of the loop.

The <condition> is tested once before the statement sequence is executed.

If it is true then the body of the loop is executed once, the condition is retested, and the

process is repeated. The loop is terminated, on completion of the statement sequence, the

first time the condition becomes false. Here are a few simple examples:

> x := 5: y := 0:

> while x > 0 do

> y:= y + x^2;

> x := x - 2;

> end do:

> y;

(Note: The colon on the "od" at the end of the statement sequence suppresses the

intermediate output generated by iterations of the loop. Change to a semicolon

if you want to trace the execution step by step.)

Here is a while loop that performs Euclid's algorithm for integer gcd's, using

the built-in Maple irem (integer remainder) function. The gcd is the last value

of the variable a.

> a := 105: b := 60:

> while b <> 0 do

> d := irem(a,b);

> a := b;

> b := d;

> end do:

> a;

more general forms

There are even more general loop forms that can be constructed using combinations

of for and while and ways to terminate execution of a loop prematurely if that is necessary.

Follow this link for more information: more on loops

Selection structure (if-then)

The basic form of the selection structure in Maple is as follows

if <condition> then

else

fi;

If the <condition> is true, then <statement sequence 1> is executed; if not, then

<statement sequence 2> is. For instance:

> for n from 0 to 5 do

> q := 2^(2^n) + 1;

> if isprime(q) then

> print(q, `is prime`)

> else

> print(q, `is not prime`)

> end if;

> end do;

For a multi-way selection of one case out of a list, you can also use a

more general form:

if <condition 1> then

elif <condition 2> then

...

else

<catch-all statement sequence>

fi;

For more information on selections, follow this link: more on selections

Maple Data Structures

Maple has a rich collection of data structures and operations on them, including most

of the standard mathematical constructs. Here are some of the most common ones:

lists

A list is an ordered collection of objects (constants, expressions, other lists, plots ...)

separated by commas, and enclosed by square brackets: [ ]. The items in a list

can be selected as follows. First we define a list of integers called

ListA, then select the third item in it (in either one of two ways).

> ListA:=[2,5,7,3,-10];

> ListA[3];

> op(3,ListA);

The built-in Maple function called op takes a list (or a set) and

returns one of the sequence of objects in it. If the first argument is omitted, then

the entire sequence of elements of the list is returned by op. This is useful, for instance, if

we want to build up a list one element at a time :

> primes:=[];

> for n to 20 do

> if isprime(n) then primes := [op(primes),n] fi;

> end do;

> primes;

or concatenate lists:

> ListB:=[4,1,6];

> ListC:=[op(ListA),op(ListB)];

Another useful method for constructing lists where the entries are elements

of some sequence defined by a definite pattern is to use the built-in Maple seq function:

> evenlist:=[seq(2*n,n=1..10)];

Note: The seq function creates the sequence of even numbers, then that is enclosed in

square brackets for the list.

sets

As in mathematical notation, a set in Maple is an unordered collection of objects, in

which an object can appear at most once. In a set, the objects are listed separated by

commas, and enclosed by braces: { }.

> SetA := {3,4,5,4,2};

(Note that the elements here have been reordered and the duplicate 4 has been removed.)

The op function discussed above is also defined for sets. In addition, the standard

set operations union, intersection, set difference are available:

> SetB := {5,6,7};

> SetA union SetB;

> SetA intersect SetB;

> SetA minus SetB;

Note that many of the familiar Maple operations take inputs that are either lists

or sets. For instance the standard form for parametric curve plotting in the

plane takes as input either a single list [x(t),y(t),t=range], or a set of these lists.

arrays (matrices)

While a rectangular array can be viewed conceptually as a list of lists (all sublists

the same size), Maple distinguishes between such lists and arrays (which have specified

subscript ranges). Here is a definition of a two-by-two matrix using the array function:

> with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> A:=matrix([[1,2],[3,4]]);

The indexing of a matrix A gives rows numbered 1..rowdim(A) and columns numbered

1..coldim(A). Individual entries can be selected like this:

> A[1,2];

If we wanted to use different ranges of indices, we could also use the array

definition:

> B:=array(0..3,0..4,[[1,2,3,4,5],[2,3,4,5,6],[3,4,5,6,7],[4,5,6,7,8]]):

> B[0,0]; B[3,4];

Matrix multiplication, addition, etc are performed by commands from the linalg package.

linalg package

conversions

The Maple convert function can be used to convert from one data type to another

For more information, follow these links:

more on conversions

Structure of Maple Procedures

A Maple procedure can be thought of as a user-written Maple program, or

equivalently as a specifications for a new command added to the basic Maple system. The

statements making up a procedure can be entered interactively during a session. But this can

be somewhat confusing, especially if you need to correct typing or syntax errors. I

HIGHLY RECOMMEND creating the statements for the procedure definition in a separate file

(using a text editor), reading them into Maple, seeing if the procedure "compiles"

correctly, and then making any necessary changes in the text editor.

syntax

Procedure definitions look like this:

<procedure name> := proc(<par-seq>)

local <local-varseq>;

global <global-varseq>;

end proc:

The <par-seq> is a sequence of variable names which are "inputs" to the procedure

They are analogous to value parameters in Pascal; their values cannot be changed

within the procedure -- you cannot even put one of them on the left of an assignment.

Each of them can be declared with a type if you want -- this is sometimes

useful as a way to catch mistakes in the information supplied to a procedure.

For instance, the procedure heading

MyProc := proc(a::integer,x::array)

says that the procedure will expect exactly two parameters, the first one an integer,

and the second one an array.

The local list consists of local variables that are used while the procedure is

executed, then discarded. Global variables are variables "outside" the procedure

whose values will be used (and possibly modified) by the procedure. The statements

can be any legal Maple commands, including calls to procedures, including the one

being defined (that is -- recursion is OK).

First examples

The simplest Maple procedures are essentially just groups of statements you

want to execute as a unit by entering a single command. For instance, to

compute gcd's of pairs of integers several times, we might want to make

a procedure out of the while loop implementation of the Euclidean algorithm

given above. (There is actually a built-in Maple function that does this too; this

is just an example!) I entered that in a file using a text editor, then read it into Maple

as follows.

> read `/home/fac/little/PuertoRico/Programs/ProcEx.map`;

The definition of the procedure can be viewed like this:

> interface(verboseproc=2);

> print(Euclid);

The procedure can be called by supplying values for the arguments (either constants

or other expressions) for instance as follows:

> Euclid(3423,2345);

The "output" from this procedure is just the output line printed here. If we wanted

to make the value of the gcd the output, so that it would be available for use in later

computations, we can do that by using a return statement in the statement

part of the procedure. Here is the modified version, and an illustration of using a value

returned by a procedure:

> print(Euclid2);

> zz:=5*Euclid2(3452,2994) + 7;

Here's an example of a recursive procedure for computing the nth Hermite polynomial

H_n(x) -- a solution of y'' - 2x y' + 2n y = 0 -- via the recurrence relation

H_n(x) = 2x H_n-1(x) - 2(n-1) H_n-2(x) and initial cases H_0(x) = 1, H_1(x) = 2 x.

> print(HermitePoly);

> HermitePoly(7,t);

>

For more on structure and syntax of Maple procedures, follow this link: more on procedures

An Extended Example

The Bisect Procedure

Here is Maple code for the bisection root-finding method

we discussed in class on September 7:

> read "/home/fac/little/public_html/Num01/Bisection.map";

> print(Bisect);

To call this procedure to compute an approximate root of an equation,

we could proceed like this.

> f:=x->x^5 - 14*x + 3;

> plot(f,-2.1..2);

> Bisect(f,0.0,1.0,.00001);

> fsolve(f(x),x,x=0..1);

If it cannot be determined that the function changes sign, an error statement

is generated. For instance, if we try running Bisect on the interval [0.4,0.6]:

> Bisect(f,0.4,0.6,.01);

Error, (in Bisect) Can't tell whether f changes sign on this interval or not -- try a different interval

The procedure still works on intervals that contain more than one root,

but it will find only one of them -- which one depends on the spacing of the points

> Bisect(f,-3.0,3.0,.000001);

> Bisect(f,-3.0,10.0,.000001);

>

An Exercise on Procedures

Exercise

A derangement of the integers {1, 2, ... , n} is permutation that leaves no integer

fixed. Write a procedure that generates a list of all derangements of {1,2,...,n},

given n as input. This can be done either "from scratch" (in which case you

need to be a bit clever!) or else by using the permute procedure from the

combinat package. Follow this link for some info: more on permute