Introduction to Maple

introMaple.mws

> restart;

A Short Introduction to Maple

>

Arithmetic Operations

Variable Assignement

> a := 100;

> b := 200;

> a+b;

> a mod 3;

> 12*23;

> 12^23;

> 30!;

> 100/20;

> 100/30;

> trunc(100/30);

> round(100/30);

> trunc(-4.3);

> round(-4.3);

> floor(2.3); floor(-2.3);

> ceil(2.3); ceil(-2.3);

> abs(-4);

Quotes, Concatenation operator

> a:=3;

> a:='a';

> a;

> sqrt(s^2);

> assume(s>0); sqrt(s^2);

> assume(s<0); sqrt(s^2);

> s:='s';

> a:=(x+y)^3;

In Maple 5, the concatenation operator was the dot.

In Maple 6, the concatenation operator is the two vertical bars.

> a||1;

> a||1||2;

Definition of strings

> c:="abcd";

Interface, Code Generation

Customize the user interface

> interface(version);

> interface(screenwidth=100);

> interface( quiet=true);

> interface( quiet=false);

> interface( verboseproc = 2 );

> eval(isprime);

LaTeX interface

> expr:=expand((x+y)^10);

> latex(expr);

{x}^{10}+10\,{x}^{9}y+45\,{x}^{8}{y}^{2}+120\,{x}^{7}{y}^{3}+210\,{x}^{6}{y}^{4}+252\,{x}^{

5}{y}^{5}+210\,{x}^{4}{y}^{6}+120\,{x}^{3}{y}^{7}+45\,{x}^{2}{y}^{8}+10\,x{y}^{9}+{y}^{10}

C and Fortran code generation

> with(codegen);

Warning, the protected name MathML has been redefined and unprotected

> C(1-x/2+3*x^2-x^3+x^4);

      t0 = 1.0-x/2.0+3.0*x*x-x*x*x+x*x*x*x;

> fortran(1-x/2+3*x^2-x^3+x^4);

      t0 = 1-x/2+3*x**2-x**3+x**4

> C([a=10.3*x-y,b=z/t*log(z)]);

      a = 0.103E2*x-y;

      b = z/t*log(z);

> fortran([a=10.3*x-y,b=z/t*log(z)]);

      a = 0.103D2*x-y

      b = z/t*log(z)

> cost(expand((x+y)^5));

Types (whattype, type, hastype)

Maple has a type system. Using it is not obligatory though.

> a:=1; whattype(a);

> b:=1.2; whattype(b);

> c:=2-3*I; whattype(c);

> c:=cos(x); whattype(c);

> whattype(1/2);

$fraction$

> whattype(x);

For lazy typists, there is an alias mechanism in Maple.

> alias(w=whattype);

> w(a+b);

> w(x+y); w(x*y); w(x^3);

> w(w);

Besides whattype, there are two more commands that handle types

> type(a,integer); type(b,integer);

> hastype(a,integer); hastype(b,integer);

The hastype(expr,t) returns true if expr has a subexpression of the type t

> pol:=(x+x*y)^2;

> w(pol);

> type(pol,`*`);

> hastype(pol,`*`);

> type(pol,`+`);

> hastype(pol,`+`);

> pol:=expand((x+x*y)^2);

> w(pol);

> type(pol,`*`);

> hastype(pol,`^`);

Access parts of a Maple expression (op, nops, indets)

Tree structure of an expression in Maple.

Internal representation of mathematical expressions.

> a:=x+y;

> nops(a);

> op(1,a);

> op(2,a);

> b:=sqrt(a);

> nops(b);

> op(1,b);

> op(2,b);

> w(b);

> nops(op(2,b));

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

Because fractions are quite common,

there are two special commands to get

the numerator and the denominator.

> numer(1/2); denom(1/2);

There is a more compact way to write this

> op([2,1],b); op([2,2],b);

> op(0,b);

> op(1..3,[5,6,7,8]);

> op(b);

> indets(x+y^2);

> a:=x+z^3: indets(a);

> indets(sqrt(x)); nops(sqrt(x));

${x, x^(1/2)}$

> indets(sqrt(x)+y^(1/3));

${x, y, x^(1/2), y^(1/3)}$

> indets(exp(x+y)); indets(log(x+y));

Sequences, Lists and Sets

A sequence is defined as a succession of elements separated by commas

> 1,2,3;

> whattype(%);

We can create fairly complicated sequences

> seq(i,i=1..10);

> seq(i^2,i=1..10);

> seq(x^i-1,i=1..7);

A list is defined as a sequence enclosed in (square) brackets

> l:=[1,2,3];

> whattype(%);

You can use seq to define new lists

> l:=[seq(ifactor(i),i=1..6)];

To access the i-th element of a list

> l[4];

A set is defined as a sequence enclosed in (curly) braces

> s:={1,2,3};

> whattype(%);

Sets cannot contain multiple elements as opposed to lists

> l:=[1,1,2,3];

> s:={1,1,2,3};

Apply a function to every element of a list (map)

> l1:=[1,9,64];

> map(sqrt,l1);

> l2:=[-1,4,-6];

> map(abs,l2);

> l3:=[10/3,2,3/7];

> map(evalf,l3);

> map(evalf[15],l3);

> f:=x->x mod 3;

> map(f,[1,10,12]);

> map(x->x^2+1,[1,10,12]);

> map(x->x^3-1,[x,y+1,z+2]);

> map(factor,%);

When the function has arguments,

these appear at the end.

> l:=[1+x^2+x^2*y+x-x^2*y^3,x^3-x^3*y+x^2];

> map(collect,l,x);

Operations on sets (union, intersect, member)

> s1 := {1,2,3}; s2 := {1,4,5};

> s1 union s2;

> s1 intersect s2;

> member(2,s1); member(2,s2);

> s1 minus s2;

> s2 minus s1;

> s1 minus {1};

> s2 union {6};

Operations on lists

Add an element to a list

> m1:=[1,2,3];

> m1:=[op(m1),4];

> m1:=[5,op(m1),6];

> m1:=[op(1..3,m1),100,op(4..nops(m1))];

This technique is not efficient when

we are dealing with big lists.

Delete an element from a list

> m2:= subsop(2=NULL,m1);

Note that the initial list has been left unaltered.

m2 is a simply copy of m1, with the second element removed.

> m1;

Modify an element of a list

> m1[2]:=20;

> m1;

> subsop(2=25,m1);

> m1;

Sums and Products (sum, prod)

> sum(i,i=1..10^6);

> add(i,i=1..10^6);

sum looks for a closed form and then substitutes the

numerical values of the bounds of summations

add performs a loop from the lower bound to the upper bound

> sum(i^2,i=1..1000000000000);

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

> factor(%);

> factor(sum(i^3,i=1..n));

> factor(sum(i^3,i=a..b));

> product(i,i=1..10); 10!;

Infinite sums and products are handled

> sum('1/i^2','i=1..infinity');

> sum('1/k!', 'k'=0..infinity);

> product('1-z^2/n^2','n=1..infinity');

product:   "Cannot show that 1-z^2/n^2 has no zeros on [1,infinity]"

Sometimes the results are not immediately readable

> product('1/i^2','i=1..n');

> i:='i':k:='k':
sum(combinat[binomial](i+k-1,k-1),i=0..n);

> sum('(2*z)/(z^2-n^2)','n=1..infinity');

use single quotes to avoid premature evaluation

> i:=2;

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

Error, (in sum) summation variable previously assigned, second argument evaluates to 2 = 1 .. 1000

> sum('i','i'=1..1000);

in the following sum we access the command binomial from the

package combinat in the long form, avoiding the with(combinat)

> n:='n';
sum(combinat[binomial](n,k),k=0..n);

the sum is actually equal to (1+1)^n (binomial theorem)

we can have nested sums and products

> factor(sum(sum('i*j','i=1..n'),'j=1..n'));

> i:='i': j:='j': n:='2':

> sum(sum(n!/(i!*j!*(n-i-j)!),i=0..n),j=0..n-i);

> evalc(%);

the sum is actually equal to (1+1+1)^2 = 9, by the trinomial theorem

now let us try to verify that (1+1+1)^3 = 27:

> i:='i': j:='j': n:='3':

> sum(sum(n!/(i!*j!*(n-i-j)!),i=0..n),j=0..n-i);

(more interesting things happen for n=4)

to actually do this sum we break it in two pieces

> i:='i': j:='j': n:='3':

> sum(n!/(i!*j!*(n-i-j)!),j=0..n-i);

> sum(%,i=0..n);

now we know how to do the general sum (1+1+1)^n = 3^n

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

> sum(n!/(i!*j!*(n-i-j)!),j=0..n-i);

> simplify(sum(%,i=0..n));

Conversions (convert)

> convert([1,2,3],set);

> convert({1,2,3,4},list);

> convert([1,1,2,3],set);

> l:=seq(i^3,i=1..5);

> convert([l],`+`);

> convert([l],`*`);

> series(tan(x),x,9);

> w(%);

> convert(%%,polynom);

> w(%);

> convert(16341,hex);

> convert(16341,binary);

> convert(16341,octal);

For more general bases, use the following syntax

> convert(13,base,6);

> convert(15,base,15);

Partial fraction decomposition of a rational function

> r:=(3*x+1)/(x^2-4*x+4); w(r);

> convert(r,parfrac,x);

Iterative constructs (for, while)

> for i from 1 to 4 do i^2; od;

> l:=[];

> for i from 1 to 4 do l:=[op(l),i^2]; od;

> i:=1;

> while i+2 < 6 do i:=i+1: od;

Functions

> f:=x->x+1;

> f(2);

> type(f,'function');
type(cos(x),'function');

> whattype(f);

The operators @ and @@ are used to compose functions

> (f@f)(2);

> (f@@5)(2);

> g:=y->y+2;

> (f@g)(x);

> ((f@g)@@3)(x);

> h:=(x,y,z)->x^2+y^2-z^2;

> h(3,4,5);

Recursive functions (factorial, Fibonacci)

> fact:=n->if n = 0 then 1; else n*f(n-1); fi;

> fact(3); fact(7);

>    fibo:=n->if n = 0 then 0
    elif n = 1 then 1
    else fibo(n-1)+fibo(n-2);
fi;

> seq(fibo(i),i=0..12);

> phi := (1+sqrt(5))/2;
Phi := (1-sqrt(5))/2;

> fibo1:=n->(phi^n-Phi^n)/(sqrt(5));

> fibo1(100);

> expand(%);

> seq(fibo(i)-expand(fibo1(i)),i=1..10);

> time(fibo(25));

> time(expand(fibo1(500)));

Procedures

> f:=proc(x) x^2; end;

> f(3);

>    ff:=proc(n,g)
   local aux,i;
      aux:=[seq(i^2,i=1..n)];
      map(g,aux);
end proc;

> ff(3,f);

The option remember updates the remember table

>    fib := proc(n) option remember;
    if n<2 then n
           else fib(n-1)+fib(n-2) end if;
end proc;

> st:=time():
fib(100);
time()-st;

Accessing the remember table

> op(4,eval(fib));

> op(4,eval(cos));

> cos(Pi/12):=1/4*sqrt(6)+1/4*sqrt(2);

> op(4,eval(cos));

> cos(Pi/12);

Matrices (linalg)

> with(linalg);

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

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

> det(m);

> im:=inverse(m);

> id3:=evalm(m&*im);

> eigenvalues(m);

> map(evalf,[%]);

> charpoly(m,lambda);

> eig:=[eigenvectors(m)];

> eiv:=eig[2][1];

> eigm2:=op(eig[2][3]);

> evalm((m-eiv*id3)&*eigm2);

> nullspace(m);

Gaussian Elimination

> m:=matrix([[1,2,3,-2],[2,3,1,4],[3,2,1,7]]);

> m:=addrow(addrow(m,1,2,-2),1,3,-3);

> m:=addrow(m,2,3,-4);

> submatrix(m,1..2,1..3);

> minor(m,1,2);

Matrices (LinearAlgebra)

> with(LinearAlgebra);

Warning, the assigned name GramSchmidt now has a global binding

> M:=Matrix(3,[[1,2,3],[-1,-3,-5],[-5,7,9]]);

> sort(CharacteristicPolynomial(M,lambda));

> Eigenvalues(M);

> f:= (i,j) -> x^i+y^j;

> Matrix(3,f);

Differentiation, Integration (diff, int)

> restart;

> diff(x^n,x);

> diff(1/sqrt(x+1),x);

> diff(x^x,x);

> diff(exp(x*ln(x)),x);

> int(x^n,x);

> int(sin(x)^5,x);

> int(exp(-x^2),x);

> int(sqrt(1-x^2),x=-1..1);

> int(1/(1+x^3),x=0..infinity);

> int(1/(1+x^6),x=0..infinity);

> int(1/(1+x^7),x=0..infinity);

> int(1/(1+x^n),x=0..infinity);

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

Need to know the sign of --> n

Will now try indefinite integration and then take limits.

> int(cos(x)/(1+x^2),x=0..infinity);

> convert(%,exp);

> simplify(%);

> int(log(x)/(1+x)^3,x=0..infinity);

Number Theory (numtheory)

> with(numtheory);

Warning, the protected name order has been redefined and unprotected

$[GIgcd, bigomega, cfrac, cfracpol, cyclotomic, divisors, factorEQ, factorset, fermat, imagunit, index, integral_basis, invcfrac, invphi, issqrfree, jacobi, kronecker, lambda, legendre, mcombine, mersen...$
$[GIgcd, bigomega, cfrac, cfracpol, cyclotomic, divisors, factorEQ, factorset, fermat, imagunit, index, integral_basis, invcfrac, invphi, issqrfree, jacobi, kronecker, lambda, legendre, mcombine, mersen...$
$[GIgcd, bigomega, cfrac, cfracpol, cyclotomic, divisors, factorEQ, factorset, fermat, imagunit, index, integral_basis, invcfrac, invphi, issqrfree, jacobi, kronecker, lambda, legendre, mcombine, mersen...$

> ifactor(2^(2^6)+1);

> divisors(111111);

> n:=2352356;

> phi(n);

> pr:=convert(factorset(n),list);

> n*product('(1-1/pr[i])','i'=1..nops(pr));

> cyclotomic(3,x);

> cyclotomic(6,x);

> cyclotomic(7,x);

> cyclotomic(8,x);

> cyclotomic(100,x);

>    for n from 1 to 500 do
   if {coeffs(cyclotomic(n,x))} minus {1,-1} <> {}
      then lprint(n); fi;
od;

> s:={};
st:=time():
for n from 1 to 2000 do
s:=s union {coeffs(cyclotomic(n,x))}: od:
time()-st;
s;

Count the number of primes less than a given limit

> pi(100);

resultants

> p:=x^6+y^2-1;

> q:=x^3+y^3-5;

> factor(resultant(p,q,x));

> solve({p,q},{x,y});

${y = RootOf(24+_Z^2-10*_Z^3+_Z^6), x = RootOf(_Z^3-5+RootOf(24+_Z^2-10*_Z^3+_Z^6)^3)}$

> with(plots):

Warning, the name changecoords has been redefined

> p:=x^2+y^2-1;q:=x-y-1;

> a:=implicitplot(p,x=-1..1,y=-1..1):

> b:=implicitplot(q,x=-2..2,y=-2..1,color=green):

> display({a,b});

other packages

> with(orthopoly);

> seq(T(n,x),n=1..6);

> map(factor,[%]);

> with(stats);with(stats[statplots]):

> mark:=rand(1..100);

> marks:=[seq(mark(),i=1..30)];

> histogram(marks);

>