Dumas Gourdon, Maple code, livre 2, chap. 2

Maple
Son bon usage en mathématiques
Philippe Dumas, Xavier Gourdon
Springer-Verlag, 1997
Code Maple du Chapitre 2, Livre second.
# section 1, page 141.


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

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

B:=matrix(4,3,proc(i,j) i+j end);

C:=array(1..3,1..3,antisymmetric):
C[3,2]:=p:C[1,3]:=q:C[2,1]:=r:
print(C);

evalm(A&*B&*C); 
linalg[multiply](A,B,C);

coeffA:=rand(-10..10):
n:=5:
A:=matrix(n,n,[seq(coeffA(),i=1..n^2)]);

sol:=NULL:
for k to n^2 while sol=NULL do
  Equ:=evalm(add(lambda[i]*A^i,i=0..k-1)+A^k);
  Inc:={seq(lambda[i],i=0..k-1)};
  sol:=solve({seq(seq(Equ[i,j],j=1..n),i=1..n)},Inc)
od:
nu:=k-1:
mu:=subs(sol,add(lambda[i]*x^i,i=0..nu-1)+x^nu);

with(linalg):
P:= companion(x^3+x^2,x):
Q:=companion(x^3,x):
R:=companion(x^3-x^2,x):
Z := matrix(3,3,0):
A:=blockmatrix(3,3,Z,-R,Q,R,Z,-P,-Q,P,Z);

n:=coldim(A):
M:=array(1..n,1..n,identity):
d[0]:=0:
for k to n do 
  M:=evalm(A&*M);
  N[k]:=nullspace(M);
  d[k]:=nops(N[k]) 
od:
seq(d[k],k=0..n);

seq(N[k],k=0..9);

N[1];
N[2];

if d[1]>0 then 
  K:=augment(op(N[1]));
  r:=d[1];
  for k from 2 to n do
    if d[k]>d[k-1] then 
      for v in N[k] do
        tmpK:=augment(K,v);
        if rank(tmpK)>r then
          K:=eval(tmpK);
          r:=r+1;
        fi;
      od
    fi
  od
fi:
eval(K);

AK[1]:=K:
for k from 2 to n do AK[k]:=evalm(A &* AK[k-1]) od:
seq(eval(AK[k]),k=1..n);

S:=T;

K:=tmpK;

        if rank(tmpK)>r then
          K:=tmpK;
          r:=r+1
        fi;

f:=(x-2)*(x+1)^2:
rec:=add(coeff(f,x,k)*u(n+k),k=0..degree(f,x));

sol:=rsolve({rec,seq(u(k)=a[k],k=0..degree(f,x))},{u(n)});

g:=q^n*n^nu:
matrix(3,3,[seq(subs(n=nn,[seq(subs(i,g),
            i=[{q=2,nu=0},{q=-1,nu=0},{q=-1,nu=1}])]),nn=0..2)]);
linalg[det](");

sol:=rsolve({rec=n},{u(n)});

remove(has,subs(sol,u(n)),u);

# section 3, page 153.

# probleme page 155.
triangularmultiply:=proc(A::matrix,B::matrix)
  local M,n,i,j,k;
  n:=linalg[coldim](A);
  M:=array(1..n,1..n,sparse);
  for i to n do 
    for j from i to n do 
      for k from i to j do M[i,j]:=M[i,j]+A[i,k]*B[k,j] od
    od
  od;
  eval(M)
end: # triangularmultiply

n:=10:
H:=matrix(n,n,proc(i,j) 1/(i+j-1) end);
# fin probleme.

# probleme page 157.
a:=(x-2)^3:
b:=(x+1)^2:
p:=a*b;
q:=collect(p,x):
d:=degree(q,x):
C:=linalg[companion](q,x);

haphazard:=rand(-10..10):
P:=matrix(d,d,[seq(haphazard(),k=1..d^2)]):
if linalg[det](P)=0 then ERROR(`P is singular. Try again!`) fi;
L:=evalm(P^(-1) &* C &* P);

evalm(subs(x=L,p));

A:=evalm(subs(x=L,a)):
B:=evalm(subs(x=L,b)):

polytodiffop:=proc(p,x,t)
  local q,d,opalg,k,z;
  global _z,_opalg;
  if type(p,polynom(ratpoly,x)) then
    q:=collect(p,x);
    d:=degree(q,x);
    opalg :=add(coeff(q,x,k)*'diff'(z,[t$k]),k=0..d);
    subs(_opalg=opalg,z=_z,proc(_z::algebraic) _opalg end)
  else ERROR(`polytodiffop expects its first argument,p, to be
              a polynom in`,x,`with coefficients of type 
              ratpoly, but received`,p)
  fi
end: # polytodiffop

polytodiffop(x^2+1,x,t)(z(t));

dsolve(polytodiffop(x^2+1,x,t)(z(t)),z(t));

a:=(x-2)^3:
b:=(x+1)^2:
A:=polytodiffop(a,x,t):
B:=polytodiffop(b,x,t):
A(y(t)),B(y(t));

A(B(y(t))); 
eval(subs(y(t)=A(y(t)),B(y(t))));
p:=a*b: P:=polytodiffop(p,x,t): P(y(t));

# section 4, page 160.

# probleme page 160.
u:=binomial(n,k)*binomial(k,m):
imax:=1: jmax:=1:
equ[1]:=add(add(a[i,j]*subs(n=n-i,k=k-j,u),j=0..jmax),i=0..imax):

equ[2]:=expand(equ[1]/u);
equ[3]:=collect(numer(equ[2]),k);

sys:={seq(coeff(equ[3],k,i),i=0..degree(equ[3],k))}:
SOL:=solve(sys,{seq(seq(a[i,j],j=0..jmax),i=0..imax)});

rec[1]:=subs(SOL,add(add(a[i,j],j=0..jmax)*s(n-i),i=0..imax)):
rec[2]:=numer(rec[1]):
sol:=rsolve({rec[2]},s(n));

subs(n=m,sol):
solve("=1,s(0)):
sol:=subs(s(0)=",sol);
convert(",binomial):
combine(",power);
# fin probleme.

# section 5, page 164.

# exercice 2.1, page 164.
n:=3:
A:=array(1..n,1..n,antisymmetric):
sol:=NULL:
for k to n^2 while sol=NULL do
  Equ:=evalm(add(lambda[i]*A^i,i=0..k-1)+A^k);
  Inc:={seq(lambda[i],i=0..k-1)};
  sol:=solve({seq(seq(Equ[i,j],j=1..n),i=1..n)},Inc)
od:
nu:=k-1:
mu:=subs(sol,add(lambda[i]*x^i,i=0..nu-1)+x^nu):
collect(mu,x);
# fin exercice 2.1.

# exercice 2.2, page 165.
f:=(x-lambda)^3:
rec:=add(coeff(f,x,k)*u(n+k),k=0..degree(f,x));

sol:=rsolve({rec,seq(u(k)=a[k],k=0..degree(f,x))},{u(n)}):
collect(sol,{n},factor);

f:=(x^2+1)*(x+1)^2:
rec:=add(coeff(f,x,k)*u(n+k),k=0..degree(f,x));

sol:=rsolve({rec,seq(u(k)=a[k],k=0..degree(f,x))},{u(n)}):
collect(evalc(sol),{cos,sin,n},factor);
# fin exercice 2.2.

# exercice 2.3, page 166.
decompose:=proc(P,basis)
  local Q,n,d,Equ,coord,i,j,k,sys,inc,sol;
  Q:=evalm(P);
  n:=linalg[coldim](Q);
  d:=nops(basis);
  Equ:=evalm(Q-add(coord[k]*basis[k],k=1..d));
  sys:={seq(seq(Equ[i,j],j=1..n),i=1..n)};  
  inc:={seq(coord[k],k=1..d)};
  sol:=solve(sys,inc);
  if sol=NULL then FAIL
  else subs(sol,add(coord[k]*basis[k],k=1..d))
  fi
end: # decompose

M:=matrix(3,3,[x+y+z,-x-y,-x-z,-x-z,x+y+z,-x-y,-x-y,-x-z,x+y+z]):
B[1]:=subs(x=1,y=0,z=0,eval(M)):
B[2]:=subs(x=0,y=1,z=0,eval(M)):
B[3]:=subs(x=0,y=0,z=1,eval(M)):
Basis:=[B[1],B[2],B[3]]:
decompose(M,Basis);

Mu:=subs(x=xi,y=eta,z=zeta,eval(M)):
P:=evalm(M &* Mu):
decompose(P,Basis);

decompose(M&*Mu-Mu&*M,Basis);

Id:=array(1 .. 3,1 .. 3,identity):
decompose(Id,Basis);

linalg[det](M);

P:=evalm(M^(-1)):
decompose(P,Basis);

C:=matrix(3,3,[0,0,1,1,0,0,0,1,0]);

newBasis:=[Id,C,C^2]:
map(decompose,newBasis,Basis);

map(decompose,Basis,newBasis);

B[0]:=array(1..2,1..2,identity):
B[1]:=matrix(2,2,[0,1,-1,0]):
B[2]:=matrix(2,2,[0,I,I,0]):
B[3]:=matrix(2,2,[I,0,0,-I]):
Basis:=[seq(B[i],i=0..3)]:
H:=alpha*B[0]+beta*B[1]+gamma*B[2]+delta*B[3]:
HH:=a*B[0]+b*B[1]+c*B[2]+d*B[3]:
decompose(H&*HH,Basis);
# fin exercice 2.3.

# exercice 2.4, page 168.
linalg[sylvester](a*x^2+b*x+c,`a'`*x^2+`b'`*x+`c'`,x):
linalg[det](");

f:=x^3+p*x+q;
linalg[sylvester](f,diff(f,x),x):
linalg[det](");

f:= 7*x^2-7*x*y+7*y^2+146*x-149*y-189:
g:=x^2-2*y^2-1:
linalg[sylvester](f,g,y);

S:=collect(linalg[det]("),x);

sol:=[solve(S,x)]:
sol:=evalc(sol):
sol:=remove(has,sol,I);

for s in sol do
  s,op(solve(subs(s,{f,g}),y)) 
od;
# fin exercice 2.4.

# exercice 2.6, page 170.
H:=matrix(2,2,[2,1,1,2]):
h:=(2*x+1)/(x+2):
HH:=H:
hh:=h:
for k to 10 do 
  HH:=evalm(HH&*H); 
  hh:=normal(subs(x=h,hh)) 
od;

H:=matrix(2,2,[a,b,c,d]):
h:=(a*x+b)/(c*x+d):
`H'`:=matrix(2,2,[`a'`,`b'`,`c'`,`d'`]);
`h'`:=(`a'`*x+`b'`)/(`c'`*x+`d'`);
evalm(H &* `H'`);

collect(normal(subs(x=`h'`,h)),x);

EV:=[linalg[eigenvects](H)]:
lambda:=EV[1][1]; 
mu:=EV[2][1];

sol:=[solve(h=x,x)];

if radnormal(sol[1]-op([1,3,1],EV[1])/op([1,3,1],EV[2]))=0
  then alpha:=sol[1]: beta:=sol[2]: 
  else alpha:=sol[2]; beta:=sol[1] 
fi:

k:=normal(((h-alpha)/(h-beta))/((x-alpha)/(x-beta)),expanded):
radnormal(k-lambda/mu);
# fin exercice 2.6.