Introduction to Maple

We introduce Maple in three sections:

1) Basic calculations

2) Calculus

3) Linear algebra.

It is meant as a crash course, and not to replace the usual texts. There is also an intro.mws file in the tutorial section to Maple which one calls by issuing ?intro at the command line. We expect the reader to have taken introductory courses in mathematics and physics.

Standard references are the two books by Michael B. Monagan et al. published by Springer:

1) Maple V Learning Guide

2) Maple V Programming Guide

Of course, some of the topics covered in those texts are dealt with in our worksheets in the physics context, which makes for more interesting reading.

>

Basic calculations and graphs

Maple as a Calculator:

Maple is case sensitive, distinguish carefully upper-case from lower-case typing.We begin by manipulating numbers:

> 1+17;

> answer:=%;

> answer;

Note that := acts as an assignment operator. The equal sign has a different meaning. A common beginner's mistake is to use the equal sign instead of the assignment operator and to get thoroughly confused.

> answ=10;

> answ;

The statement answ=10 is an equation, which was entered above, but is not usable. The variable answ remains unassigned.

We could assign

> ans:=3/4;

> evalf(ans);

The Digits variable controls how many digits are carried in floating-point calculations. It can be changed to any integer number.

> Digits;

> Digits:=20;

We can evaluate Pi now to 20 digits:

> evalf(Pi);

The symbols for Greek letters are obtained by typing their English names.

An awkward matter is that lower-case pi is just a symbol, and not the mathematical pi (that is given as Pi). But as symbols they look the same:

> Pi,pi;

> evalf(%);

Exponentiation:

> ans:=3^4;

What happens when we combine symbols into a formula?

> a:=Pi;

> b:=r^2;

> A:=a*b;

> evalf(A,5);

What if we need the answer for a specific value of r ? Several approaches are possible:

> eval(A,r=5);

> evalf(eval(A,r=5),4);

> evalf(subs(r=5,A));

But the following cannot work, as the second argument in evalf is reserved for the precision (number of digits):

> evalf(A,r=5);

Error, wrong number (or type) of parameters in function evalf

If we don't match brackets, we get some help, but it takes some practice to understand it.

> evalf(subs(r=5,A,4);

Error, `;` unexpected

We are physicists, and often we wish to use units. One way to use them in Maple explicitly is to enter them as variables, and to set them apart from the actual variables by a distinctive naming convention. For example:

> area:=evalf(subs(r=5*_cm,A),4);

It is possible to assign variables:

> r:=5*_cm;

> evalf(A,3);

> A;

and then to unassign them:

> r:='r':

> A;

The best practice for work is to keep variables unassigned (or general) until required, and to use the subs or eval command to evaluate for specific numbers (which is needed for plotting). The subs command is used more often when several variables are to be substituted; note that they will be substituted in a certain order. This ordering becomes important when one substituted variable depends on the other.

> subs(_cm=1,Pi=3.14,r=4*_cm,A);

> subs(Pi=3.14,r=4*_cm,_cm=1,A);

As a final step in this initial session we point out the dangers in assigning/unassigning variables before the final steps of forming an expression:

We look at the volume of a cone, and begin with the area of the circle.

> A;

Suppose we played around with particular values:

> r:=sqrt(7);

> A;

Now we define the volume:

> V:=1/3*A*h;

Now we unassign r :

> r:='r':

> A;

This looks as before. However:

> V;

The expression for V was formed when r was assigned a specific value. This value was used in the definition of V , or more precisely, the actual expression for A was substituted when V was defined, and thus, V has no connection with r anymore. Unassigning r return A to the general expression (it was defined when r was unassigned, but cannot turn back the wheel for V !

We need to assign the volume again, and since A has been turned back to the original general expression, it works now:

> V:=1/3*A*h;

Now we know how to use Maple as a smart calculator. We can perform graphs of expressions, i.e., we have a graphing calculator at this point:

> plot(A,r=0..4);

Combining several plots into one graph is easy, and we use a so-called list construct to pass two elements into the first argument of the plot procedure:

> plot([sin(x),cos(2*x)],x=0..2*Pi,color=[blue,green]);

If you leave a variable unassigned beyond the one used on the x -axis, be ready for an empty plot (or the 'empty plot' error message):

> plot(V,h=0..1);

For a 3-dimensional graph we use plot3d

> plot3d(V,r=0..4,h=0..1,axes=boxed,style=patchcontour,shading=zhue,orientation=[-130,75]);

>

>

Maple knows Calculus

We begin by resetting Maple. This comes close to quitting the program and restarting it (we lose all variables, but we do not recover all the memory used up).

Many worksheets have this command at their beginning, as you may either want to clear memory from a previous calculation within a different worksheet, or re-run the worksheet independently.

> restart;

> y:=exp(3*x);

Note that y is an expression in Maple. We will learn how to use mappings in Maple further below.

> yp:=diff(y,x);

The anti-derivative (indefinite integral) is calculated by

> Y:=int(y,x);

A definite integral between constant boundaries results in a number:

> int(y,x=0..2);

> evalf(%,4);

Euler's number is printed as a non-italicized letter e.

> exp(1);

This is not to be confused with

> e;

which is the unassigned variable e . This sounds a bit like Pi and pi, except that one enters Euler's number via exp(1).

> E;

Note that the subs command introduced earlier does not force evaluation:

> subs(x=0,y);

> simplify(%);

Maple knows something about limits.

> limit(y,x=infinity);

> limit(y,x=-infinity);

Let us see why we need simplify quite often.

> s1:=sin(x);

> s2:=cos(x);

> s3:=s1^2+s2^2;

> s3;

> simplify(s3);

Keep in mind that the latter command did not change the content of the variable s3. For that purpose an assignment is required.

> s3;

> s3:=simplify(s3);

The above simplification can actually be obtained in another way.

> s3:=s1^2+s2^2;

> combine(s3);

or more specifically

> combine(s3,trig);

The combine command is useful in other contexts where we wish to control the display of expressions.

> y:=exp(s)*exp(t);

> simplify(y);

> combine(y);

> expand(%);

The combine and expand commands allow us to switch hence and forth.

Often it is necessary to pull out parts of expressions. While this is tedious at times, and one tries to avoid it, here is the trick.

An expression has a certain number of operands:

> nops(y);

One can pull out the operands one-by-one:

> op(1,y);

> op(2,y);

Suppose we want the argument of the exp function in the second operand:

> nops(op(2,y));

> op(1,op(2,y));

There is no harm in practicing this a bit on complicated expressions. It shows how Maple assembles expressions in a tree-like hierarchical fashion.

Now we introduce mappings. They are not strictly required for simple function assignments, as we could do everything using expressions. Nevertheless they allow more flexibility in substituting arguments, and provide an easy first step towards subroutines (called procedures in Maple).

> g:=t->exp(-t);

> g;

> g(x);

In the latter step we used the mapping g to produce the expression in x . We can plot expressions (as before), and mappings as well, but the latter work differently. Be careful to observe what works, and what doesn't below:

> plot(g(x),x=-1..5);

> plot(g,-2..3);

> plot(g,x=-1..5); # this used to give an empty plot

> plot(g(x),-1..5);

Plotting error, empty plot

Many beginners get very confused by the above. There is always a logical reason when an empty plot is produced, or when a plot fails!

It is safer to work with expressions in Maple. Often we manipulate them, and then want to use them as mappings in the end. There is a command to turn expressions into mappings:

> y:=exp(-2*x)*sin(x)+exp(x)*cos(x);

> ym:=unapply(y,x);

> evalf(subs(x=3,y));

> evalf(ym(3));

Mappings can be differentiated directly using the D differential operator. We do not need to specify the independent variable, as it is just a dummy placeholder.

> ypm:=D(ym);

> plot([ym,ypm],-0.5..1.5,color=[blue,green],axes=boxed);

Note how the function y (blue) has an extremum point (maximum) around x =0.75, and the derivative (green) vanishes there. Make other observations about the function and its derivative.

Note that the display of differentiated mapping may appear confusing, as mappings are referred to by just the function name:

> D(sin);

Mappings can be combined in Maple.

> g1:=t->exp(-t);

> g2:=s->sin(2*s);

> g3:=g1@g2;

> g3(x);

Of course, the order is important as we are nesting functions (mappings) here. Also note that we need brackets if we want to evaluate directly the combined mapping at some value of the independent variable.

> (g2@g1)(x);

We can test out our knowledge of functions and their inverses:

> (ln@exp)(x);

> simplify(%);

> combine(%);

> expand(%);

Now, Maple can be very fussy, which can drive physicists up some trees... In those cases use carefully the following (which simply assumes nice things about the variables, which might be true - or sometimes not!)

> simplify(%,symbolic);

We will have to learn about the assume system, which allows to make assumptions about variable in a controlled way.

> assume(x,real);

> ln(exp(x));

> ln(exp(s));

If we wish to drop the assumption on x , we unassign it.

> x:='x':

> exp(ln(x));

> ln(exp(x));

Mappings can also be iterated.

> y:=sin@@3;

> y(x);

> y(1);

> evalf(%);

> evalf(sin(sin(sin(1))));

Note that Maple reserves the notation sin^ n for the case of the repeated application of the sine function. It will never display sine raised to the n th power as is usually done in Math texts, as this leaves room for ambiguity. sin^2( x ) in Maple means sin(sin( x )), not (sin( x ))^2.

The system of assumptions is important in integral evaluations. Suppose we wish to obtain:

> int(exp(-alpha*r),r=0..infinity);

Definite integration: Can't determine if the integral is convergent.

Need to know the sign of --> alpha

Will now try indefinite integration and then take limits.

> assume(alpha>0);

> int(exp(-alpha*r),r=0..infinity);

Maple has smartened up compared to previous versions, and tells us that assumptions are needed to obtain the result. It normally displays variables with assumptions on them by a tilde attached to the variable, but this feature can be turned off [easiest through the menu, or via interface(showassumed)=0 ]. However, we can always find out about assumptions:

> about(alpha);

Originally alpha, renamed alpha~:

  is assumed to be: RealRange(Open(0),infinity)

We finish this section by a real-life physics example. We solve the Newton equation for free fall at the surface of a planet (constant acceleration).

> restart;

> NE:=m*diff(v(t),t)=m*g;

We begin by writing the integrals in unevaluated form (upper-case version of int ). We have to introduce a new dummy variable s for time, in order not to confuse the integration boundary with the integration variable. We make use of the lhs and rhs functions to pull out the sides of the equation . The mass m is divided out of the equation.

> NE1:=Int(lhs(NE/m),t=0..s)=Int(rhs(NE/m),t=0..s);

We can carry out the integration using the value command.

> NE1:=value(NE1);

Now we introduce the relationship between velocity and position (keeping s as the time variable)

> NE2:=diff(x(s),s)=solve(NE1,v(s));

> NE2:=Int(lhs(NE2),s=0..t)=Int(rhs(NE2),s=0..t);

> sol:=value(NE2);

> solve(sol,x(t));

We used the solve command to isolate x ( t ).

Maple has a built-in differential equation solver, which can crack a large percentage of solvable linear differential equations.

First we re-define NE as the differential equation for position as a function of time:

> NE:=m*diff(x(t),t$2)=m*g;

> sol:=dsolve(NE,x(t));

No boundary (initial) conditions were specified, and therefore we have the general solution with two integration constants.

x ( t ) has not been assigned. There are two options. The more convenient one is actually to not assign x(t), but to extract the solution as an expression:

> rhs(sol);

> x(t);

> assign(sol);

> x(t);

If one wants to use x(t) again (even for re-defining a differential equation), one has to unassign (this is why never making the assignment is more convenient):

> x(t):='x(t)':

One can also find the special solution by specifying initial conditions (either as numbers or as symbols). One needs the D differential operator to specify the condition for the derivative.

> IC:=x(0)=x0,D(x)(0)=v0;

One solves the set of differential equation(s) and initial conditions:

> sol:=dsolve({NE,IC},x(t));

This used to work in previous versions, and does again in Maple 6. In Maple 5 one needs to give Maple a hint:

> sol:=dsolve({NE,IC},x(t),method=laplace);

It is useful to assign dsolve to a variable. In that case when the empty set is returned (which means: "I couldn't come up with a solution", not "There is no solution") one sees it right away (otherwise an empty line is produced).

>

>

and some linear algebra too!

In this part we deal with Maple's handling of linear algebra problems. Vectors play an important role in physics. We learn that Maple has (at least) two ways of dealing with them.

The first way is through objects called lists. Lists can be ordered (tables) or unordered. They allow to append objects, similarly as sets represent accumulations of objects.

Vectors are something more, there are mathematical operations defined on them which are performed by arithmetic on the components. Maple has an array construct to deal with vectors and matrices, and a special package called linalg that contains many useful procedures to compute operations such as determinant, matrix inverse, eigenvalues, etc. In Maple 6 a new linear algebra package has been added; it is called LinearAlgebra , and allows superior handling of numerical linear algebra problems. For compatibility with Maple 5 we will refrain from using the new package for the time being.

We first show the pitfalls of using simple lists:

> restart;

> xv:=[1,4,9];

> xv[1]; xv[2]; xv[3];

> print(xv);

> eval(xv);

> xv;

This looks easy: we defined a quantity xv , and were able to access its components using square brackets. Maple refers to this construct as a list.

> whattype(xv);

There was no need to declare the object. However, it is clear from the first line that Maple had enough information to deal with the construct, namely an ordered list of 3 objects. We cannot trivially add another, e.g., xv[4]:an error message is produced by referring to the object:

> xv[4];

Error, invalid subscript selector

The reference to an element of the list using the subscript in square brackets is equivalent to the op() construct:

> nops(xv);

> op(2,xv);

Suppose we wish to expand the list by adding another element. The following doesn't work, as it generates a list in a list.

> xv1:=[xv,16];

We need to first extract the elements from the previous list. This can be performed the long way or with a shortcut (see ?list for details):

> op(1..nops(xv),xv);

> op(xv);

We are ready to expand our list and re-use its name (be careful with recursive definitions).

> xv:=[op(xv),16];

> xv[4];

For many purposes the work with lists in Maple is the most efficient way to deal with multi-component objects. It is possible for Maple to run out of steam when lists contain 100 entries or more.

We note that the following practice does something different:

> yv[1]:=1;

> yv[2]:=8;

> yv[3]:=27;

> print(yv);

> whattype(yv);

> yv[2];

Tables indexed by integers are ordered, but the way they are stored and displayed is not always predictable.

There are functions to convert between tables and lists; in older Maple versions one had to watch what happened when the entries were out of order.

Tables are extremely flexible: they provide a very general pointer structure, whereby the indexing does not have to be continuous. In fact the index can almost be anything. The table simply provides the link.

> yv[5]:=101;

> yv[4]:=87;

> yv[i]:=2*i;

> yv[j]:=2*j-1;

> print(yv);

> convert(yv,list);

It appears that the numerically indexed entries are coming out in proper order, and that the order of generally defined entries is not quite predictable.

Arrays in Maple.

To deal with vectors and matrices where the indexing is in terms of integer numbers Maple has the array construct. In the linalg package there are special procedures called vector and matrix which define the required arrays. It is possible to define arrays outside the linalg package. Since we are interested in the functions contained in the linalg package, we use the vector and matrix constructs from it to define the arrays.

> with(linalg);

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

We displayed the commands that become available with this package.

> xv:=vector([1,4,9]);

> yv:=vector([1,8,27]);

These are column vectors displayed as rows to save space.

One can refer to vector elements using an index in square brackets, which is analogous to list elements, but the op construct works differently:

> xv[2];

> nops(xv);

Note that this is not the dimension of the vector! There is a special command for that:

> vectdim(xv);

> op(xv);

> op(op(xv));

> whattype(xv);

If one had to access the elements via the op construct, it would be a real nuisance (as compared to lists).

The inner product can be calculated as follows:

> dotprod(xv,xv);

We can write a small loop to calculate this:

> dp:=0:

> for i from 1 to 3 do:

> dp:=dp+xv[i]*xv[i]; od:

> dp;

> dotprod(xv,yv);

Edit the above loop and verify the calculation.

Note that we have referred to the elements of the vector in the same way as list elements were addressed. In fact, we could write our own set of procedures, and always work with lists. We turn the do-loop above into a procedure.

We pass two arguments (the vectors to be multiplied) and check with an if-then-fi construct whether they match.

> dprod:=proc(arg1,arg2) local i,n1,n2,dp;

> n1:=vectdim(arg1):

> n2:=vectdim(arg2):

> if n1 <> n2 then RETURN(`Mismatch of sizes in dprod procedure`); fi;

> dp:=0:

> for i from 1 to n1 do:

> dp:=dp+arg1[i]*arg2[i]: od:

> dp;

> end:

> dprod(xv,xv);

Note how a procedure returns the last executed statement.

It also has its error checking (which could be improved to check the type of arguments passed into the procedure).

> dprod(xv,[1,1,1]);

> evalm(xv);

This worked even though [1,1,1] is a list and not a vector. vectdim must be working properly on lists:

> vectdim([1,2,3,4]);

> dprod(xv,[1,2,3,4]);

> dotprod(xv,[1,2,3,4]);

Error, (in dotprod) arguments not compatible

Maple's dotprod procedure works also with lists:

> dotprod([1,2,3],[1,2,3]);

We now define a matrix:

> A:=matrix([[1,2,3],[4,5,6],[7,8,9]]);

We can check some properties:

> trace(A);

We obtained the sum of the diagonal elements. Next we calculate the determinant:

> det(A);

Obviously the rows/columns are not independent of each other. Thus, no inverse exists.

We try a small modification.

Observe how elements of a matrix are addressed, and how vectors and matrices are not displayed by default, but evalm has to be invoked.

> B:=A;

> B[3,3]:=10;

> B;

> A;

> evalm(B);

> det(B);

> A; evalm(A);

> evalm(A-B);

> B;

> A;

Something apparently absurd happened:

We assigned B to pick up the content of A . Then we changed one element in B . That should not have affected A , but obviously it did.

> Binv:=inverse(B);

Note that since the inverse exists, the system of equations

B x = b

has a unique solution. It can be obtained as follows: we pick some right-hand-side b :

> bv:=vector([1,3,5]);

> xsol:=evalm(Binv &* bv);

Observe that the &* operator works inside the evalm call.

We verify the solution (a mathematician accepts a solution only after verification):

> evalm(B &* xsol - bv);

Thus we have found a vector xsol such that B&*xsol = bv .

Of course, in practice one would not perform a matrix inverse to solve a single system of equations, as it requires n times more operations to find an inverse of an n -by- n matrix than to solve the linear system.

We demonstrate now that a homogeneous system of linear equations has interesting solutions if the coefficient matrix is singular (otherwise the solution is unique, and only the trivial solution x = [0,0,0] is allowed).

We redefine our matrix A :

> A:=matrix([[1,2,3],[4,5,6],[7,8,9]]);

> zerov:=vector([0,0,0]);

> xsol:=vector([x1,x2,x3]);

> evalm(A&*xsol-zerov);

> solve(%,xsol);

Error, (in solve) invalid arguments

This error message means that the solve procedure does not accept vectors as arguments.

> eqs:=evalm(A&*xsol-zerov);

> whattype(eqs);

> eqs[1];

> sol:=solve(eqs,xsol);

This didn't work. solve does not work on lists. We need to convert them to sets (it escapes me why one should need to do this):

> convert(eqs,set);

> sol:=solve(convert(eqs,set),convert(xsol,set));

Note that x3 is a parameter (picking a number makes the solution unique), while x1 and x2 are determined from it.

We now proceed with the linear algebra eigenvalue problem. In classical mechanics it is needed in the calculation of the principal axes of a rigid body (diagonalization of the moments-of-inertia tensor).

Given a square ( n -by- n ) matrix A we are seeking a solution of the problem

A x = lambda x

where x are n -by-1 column vectors, and lambda is a scalar that can be complex-valued.

One can prove that if the matrix A is real and symmetric, the eigenvalues lambda are real-valued. There are always n of them, although the same answer can appear with a multiplicity higher than 1.

The proof that a solution to the problem can be found is by construction. Consider the matrix

B = A - lambda I

where lambda is an unknown scalar, and I the unit matrix of size n -by- n .

The above problem then is a homogeneous system of equations with coefficient matrix B .

A homogeneous system has non-trivial solutions only if its coefficient matrix is singular.

Thus, we can find solutions to the eigenvalue problem when det( B )=0. This becomes a condition on lambda, which can then be found in a straightforward way as shown below.

The solution vectors x to be found separately for each eigenvalue lambda(i) are determined from the solution of the homogeneous system of equations B x = 0. They will not be unique in the sense that the length of the vector is arbitrary.

The geometric interpretation of the problem is as follows:

For an arbitrary vector x we have the property that y = A x points in a direction different from x .

For a given matrix A we are looking for solution vectors such that A does not rotate them. The only effect of A on these vectors is a stretch operation. The amount of stretching is given by lambda.

We first define a matrix that is real and symmetric. We can simply add the transpose of a non-symmetric matrix to it:

> Am:=evalm(transpose(A)+A);

> I3:=array(identity,1..3,1..3);

> evalm(I3);

> Bm:=evalm(Am-lambda*I3);

> cpoly:=det(Bm);

The characteristic polynomial is of order n , i.e., in our case of order 3. The fundamental theorem of algebra asserts that n roots can be found over the field of complex numbers to a polynomial of order n .

> solve(cpoly,lambda);

> evalf(%);

The eigenvalues are real, as we constructed Am to be symmetric. The fact that 0 is an eigenvalue implies that the original matrix Am is singular:

> det(Am);

The solution of roots of polynomials with higher n is difficult (impossible in closed form in general for n >5) and one has to resort to numerical techniques.

For low n they can be found as explained above by using the following routines:

> charpoly(Am,mu);

> eigenvals(Am);

To find the eigenvectors requires more work, but for small n the following works (if the roots are simple):

> eigenvects(Am);

The structure of the answer is: eigenvalue, multiplicity, set of eigenvectors (even if the multiplicity is 1);

The answer for lambda=0 is easy to read.

To obtain numerically calculated eigenvalues (and eigenvectors) consult the following help page:

> ?Eigenvals

>

>