Maple
Son bon usage en mathématiques
Philippe Dumas, Xavier Gourdon
Springer-Verlag, 1997
Code Maple du Chapitre 3, Livre second.
# section 1, page 173.
reflection3d:=proc(normalvector::{list,vector})
local d,nv;
if type(normalvector,list) then
nv:=vector(normalvector)
else nv:=normalvector
fi;
d:=linalg[vectdim](nv);
if d=0 then
ERROR(`reflection3d expects its argument normalvector to be
of size>0 but received`,eval(normalvector))
fi;
evalm(array(1..d,1..d,identity)
-2/linalg[dotprod](nv,nv)*nv&*linalg[transpose](nv))
end: # reflection3d
N:=[alpha,beta,gamma]:
reflection3d(N);
rr:=rand(-10..10):
d:=4:
M:=matrix(4,4,rr);
for i to d do
E[i]:=array(1..d,sparse);
E[i][i]:=1
od:
MM:=M:
for k to d-1 do
A:=proj(linalg[col](MM,k),k);
S:=s(A,E[k]);
MM:=map(radnormal,evalm(S&*MM));
print(MM);
print(map(evalf,MM))
od:
s:=proc(a,e)
local nu,d,alpha,epsilon,sigma,tau,v;
nu:=linalg[dotprod](a,a);
if nu=0 then
d:=nops(convert(eval(e),list));
RETURN(array(1..d,1..d,identity))
fi;
alpha:=linalg[dotprod](e,a);
if evalf(alpha)<0 then epsilon:=-1 else epsilon:=1 fi;
sigma:=epsilon*linalg[dotprod](a,a)^(1/2);
tau:=1/sigma/(sigma+alpha);
v:=evalm(a+sigma*e);
map(radnormal,
evalm(&*()-2/linalg[dotprod](v,v)*v&*linalg[transpose](v)))
end: # s
proj:=proc(v,k)
local d,i;
d:=nops(convert(eval(v),list));
array(1..d,[seq(0,i=1..k-1),seq(v[i],i=k..d)])
end: # proj
# section 2, page 178.
# exercice 3.3, page 179.
R1:=rotation3d([1,1,2],Pi):
o3danalysis(-R1);
R2:=rotation3d([1,1,0],Pi/3):
o3danalysis(R2^6),o3danalysis(-R1^2);
o3danalysis(-R1&* R2^2);
`value/Identity`:=proc(expr)
subs(Identity=proc(d) array(1..d,1..d,identity) end,expr)
end: # value/Identity
`value/&*`:=proc(expr) evalm(map(value,expr)) end:
value(Rotation3d([0,0,1],Pi/3)&*Reflection3d([1,0,0]));
R:=rotation3d([0,0,1],2*Pi/3):
S:=rotation3d([1,0,0],Pi):
Sigma:=reflection3d([0,0,1]);
# fin exercice 3.3.
# section 3, page 181.
# probleme page 181.
A:=matrix(4,4,proc(i,j) i+j end):
rotation(A,1,2);
d:=5:
A:=linalg[hilbert](5);
# fin probleme.
# section 4, page 183.
# probleme page 183.
Q:=2*(x[1]+x[2])^2-(x[1]-x[3])^2+x[1]*x[4];
gauss1(Q,x[1],'Q1');
Q1;
Q:=x*y+y*z+z*x;
gauss2(Q,x,y,'Q1');
Q1;
# fin probleme.
# section 5, page 186.
# exercice 3.1, page 186.
rotation3d:=proc(axisvector::{list,vector},theta)
local av,d,norm2,n;
if type(axisvector,list) then #1
av:=vector(axisvector)
else av:=axisvector
fi;
d:=linalg[vectdim](av);
if d<>3 then ERROR(`rotation3d expects its argument
axisvector to be of size 3 but
received`,eval(axisvector))
fi;
norm2:=linalg[dotprod](av,av);
if norm2=0 then ERROR(`rotation3d expects its argument
axisvector to be nonzero but
received`,eval(axisvector))
fi;
n:=evalm(av/sqrt(norm2)); #2
evalm(cos(theta)*array(1..3,1..3,identity) #3
+(1-cos(theta))*n&*linalg[transpose](n)
+sin(theta)*array(1..3,1..3,antisymmetric,
[(3,2)=n[1],(1,3)=n[2],(2,1)=n[3]]))
end: # rotation3d
normal(linalg[trace](rotation3d([a,b,c],theta)));
# fin exercice 3.1.
# exercice 3.2, page 187.
e:=proc(n) cos(t)^n end:
scalarproduct:=proc(v,w)
option remember;
int(v*w,t=-Pi..Pi)/Pi
end: # scalarproduct
B:=10:
f(0):=e(0);
for n to B do
f(n):=e(n)-add(scalarproduct(f(k),e(n))
/scalarproduct(f(k),f(k))*f(k),k=1..n-1)
od:
for n to B do combine(f(n),trig) od;
# fin exercice 3.2.
# exercice 3.3, page 188.
o3danalysis:=proc(A)
local AA,p,q,I3,Delta,tr,v,n2,theta,asp,axisvector,toto;
AA:=evalm(A); #1
p:=linalg[rowdim](AA);
q:=linalg[coldim](AA);
if p<>q or p<>3 then
ERROR(`o3danalysis expects its argument A to be a square
matrix, but received`,eval(AA))
fi;
I3:=array(1..3,1..3,identity);
if not linalg[iszero](map(radnormal,
evalm(AA &* transpose(AA)-I3)))
then ERROR(`o3danalysis expects its argument A to be an
orthogonal matrix of size 3, but received`,eval(AA))
fi;
Delta:=linalg[det](AA); #2
if radnormal(Delta-1)=0 then #3
tr:=linalg[trace](AA);
if radnormal(tr-3)=0 then Identity(3)
else
if radnormal(tr+1)=0 then
v:=op(1,linalg[nullspace](AA-I3));
n2:=linalg[dotprod](v,v);
theta:=Pi
else
asp:=evalm((AA-transpose(AA))/2);
v:=[asp[3,2],asp[1,3],asp[2,1]];
n2:=linalg[dotprod](v,v); #a
theta:=arccos(radnormal((tr-1)/2,rationalized))
fi;
axisvector:=map(radnormal,evalm(v/sqrt(n2)),rationalized);
Rotation3d(eval(axisvector),theta) #b
fi
else #4
toto:=o3danalysis(-AA);
if toto=Identity(3) then -Identity(3)
elif op(2,toto)=Pi then Reflection3d(op(1,toto))
else axisvector:=evalm(-op(1,toto));
Reflection3d(eval(axisvector)) &*
Rotation3d(eval(axisvector),op(2,toto))
fi
fi
end: # o3danalysis
`value/Rotation3d`:=proc(expr)
subs(Rotation3d=rotation3d,expr)
end: # value/Rotation3d
`value/Reflection3d`:=proc(expr)
subs(Reflection3d=reflection3d,expr)
end: # value/Reflection3d
G[1]:=eval(R): #1
G[2]:=eval(S):
G[3]:=eval(Sigma):
globalcounter:=3: #2
todocounter[1]:=3:
for i while todocounter[i] > 0 do
todocounter[i+1]:=0; #3
for j from globalcounter-todocounter[i]+1 to globalcounter do#4
for k to 2*j-1 do #5
if k globalcounter then
globalcounter:=c;
G[globalcounter]:=eval(g1g2);
todocounter[i+1]:=todocounter[i+1]+1;
fi
od
od
od:
seq(todocounter[i],i=1..4);
seq(o3danalysis(G[i]),i=1..globalcounter);
# fin exercice 3.3.