Maple
Son bon usage en mathématiques
Philippe Dumas, Xavier Gourdon
Springer-Verlag, 1997
Code Maple du Chapitre A, Annexe.
# section 1, page 371.
# probleme page 371.
decimalexpansion:=proc(r::rational,l::nonnegint)
local p,q,pp,k,de;
p:=numer(r); #1
q:=denom(r);
pp:=irem(p,q);
if pp<0 then pp:=pp+q fi;
for k to l do de[k]:=iquo(10*pp,q,'pp') od; #2
[seq(de[k],k=1..l)]
end: # decimalexpansion
nicedecimalexpansion:=proc(r::rational)
local p,q,pp,history,k,j,de,i;
p:=numer(r);
q:=denom(r);
pp:=irem(p,q);
if pp<0 then pp:=pp+q fi;
history[pp]:=0; #1
for k do
de[k]:=iquo(10*pp,q,'pp');
if assigned(history[pp]) then #2
RETURN([seq(de[i],i=1..history[pp])],
[seq(de[i],i=history[pp]+1..k)])
else history[pp]:=k; #1
fi
od
end: # nicedecimalexpansion
for i to 6 do
NDE:=nicedecimalexpansion(xx[i]);
patternsize[i]:=nops(NDE[2])
od:
seq(patternsize[i],i=1..6);
1+4/10+sum(142857/10^(6*k+1),k=1..infinity);
pleasantdecimalexpansion:=proc(r::rational)
local p,q,pp,qq,qqq,alpha,beta,mu,t,k,de;
p:=numer(r);
q:=denom(r);
pp:=irem(p,q);
if pp<0 then pp:=pp+q fi;
qq:=q;
for alpha from 0 while irem(qq,2,'qqq')=0 do qq:=qqq od; #1
for beta from 0 while irem(qq,5,'qqq')=0 do qq:=qqq od;
mu:=max(alpha,beta);
if qq=1 then t:=1 else t:=nu(10,qq) fi; #2
for k to mu+t do #3
de[k]:=iquo(10*pp,q,'pp')
od;
[seq(de[k],k=1..mu)],[seq(de[k+mu],k=1..t)]
end: # pleasantdecimalexpansion
# fin probleme.
# probleme page 373.
carmichael:=proc(N::integer)
local NN,f,p;
if N<2 then RETURN(false) fi;
NN:=op(2,readlib(ifactors)(N)); #1
if nops(NN)=1 then RETURN(false) fi; #2
for f in NN do #3
if op(2,f)>1 then RETURN(false) fi
od;
for f in NN do #4
p:=op(1,f);
if irem(N-1,p-1)<>0 then RETURN(false) fi
od;
true;
end: # carmichael
b:=2:
p[0]:=1: B:=1:
for i to b do
p[i]:=nextprime(p[i-1]);
p2[i]:=p[i]^2;
B:=p2[i]*B
od:
S:={seq(j,j=1..B)} minus
{seq(seq(p2[i]*k,k=0..iquo(B,p2[i])),i=1..b)}:
evalf(nops(S)/B);
BB:=10^5:
counter:=0:
for k to iquo(BB,B) do
for s in S do
n:=B*k+s;
if carmichael(n) then
counter:=counter+1;
T[counter]:=n
fi
od
od:
[seq(T[i],i=1..counter)];
threefactors:=proc(p::integer)
local q,r,x,y,aux,counter,T,i;
if not isprime(p) then ERROR(`threefactors expects its
argument to be a prime number`) fi;
counter:=0;
for x from 2 to p-1 do
for y from iquo(p^2,x)+1 to iquo(p*(p+1),x)+1 do
q:=1+(p-1)*(p+x)/(x*y-p^2);
if (type(q,integer) and q>p and isprime(q)) then
r:=(1+(q-1)*y)/p;
if type(r,integer) then
if r>p and isprime(r) and irem(q*r-1,p-1)=0 then
counter:=counter+1;
T[counter]:={p,q,r}
fi
fi
fi
od
od;
{seq(T[i],i=1..counter)}
end: # threefactors
threefactors(ithprime(10));
[seq(op(select(proc(s) evalb(convert(s,`*`)<10^6) end,
threefactors(ithprime(i)))),i=1..25)];
map(proc(s) convert(s,`*`) end,");
largecarmichael:=proc(primesetset)
local primeset,p,C,L,D,F,i,T,counter,opprimeset,TT;
global globalB;
for primeset in primesetset do
C:=mul(p,p=primeset);
L:=ilcm(seq(p-1,p=primeset));
if irem(C-1,L,'D')<>0 then
ERROR(`largecarmichael expects its argument to be
a set of primes whose product is
a Carmichael number`)
fi;
counter:=0;
for F from 2 to D while F < globalB do
if (irem(D,F)=0 and isprime(F*L+1)) then
counter:=counter+1;
T[counter]:=F*L+1
fi
od;
opprimeset:=op(primeset);
TT[primeset]:=seq({opprimeset,T[i]},i=1..counter)
od;
{seq(TT[primeset],primeset=primesetset)}
end: # largecarmichael
globalB:=1000:
SS[0]:=threefactors(ithprime(100)):
for k to 5 do SS[k]:=largecarmichael(SS[k-1]) od:
map(convert,",`*`):
max(op(map(length,")));
# fin probleme.
# section 2, page 377.
# probleme page 377.
Lambda:=array(1..2,[1,1]):
A:=array(1..2,1..2,[[0,1],[1,1]]):
Gamma:=array(1..2,[1,0]):
n:=10^5:
start:=time():
L:=convert(n,base,2):
AA:=A:
M:=&*():
for i to nops(L) do
if L[i]=1 then M:=evalm(M&*AA) fi;
AA:=evalm(AA&*AA)
od:
evalm(Lambda &* M &* Gamma):
T[n]:=time()-start:
seq(evalf(T[2^k-1]/log(2^k-1),3),k=2..10);
# fin probleme.
# probleme page 379.
with(linalg):
echelon := proc(A)
local M,p,q,i,j,line;
M:=A;
p:=rowdim(M);
q:=coldim(M);
line:=1; #1
for j from 1 to q do #2
for i from line to p while M[i,j]=0 do od; #3
if i <= p then #4
M:=swaprow(M,i,line); #5
M:=mulrow(M,line,1/M[line,j]);
M:=pivot(M,line,j); #6
line:=line+1;
fi
od;
eval(M)
end: # echelon
echelon := proc(A::matrix(rational))
local M,p,q,i,j,line,tmp,jj,c;
M:=copy(A);
p:=linalg[rowdim](M);
q:=linalg[coldim](M);
line:=1; #1
for j from 1 to q do #2
for i from line to p while M[i,j]=0 do od; #3
if i<=p then #4
for jj from j to q do #5
tmp:=M[line,jj];
M[line,jj]:=M[i,jj];
M[i,jj]:=tmp
od;
c:=M[line,j];
for jj from j to q do
M[line,jj]:=M[line,jj]/c
od;
for i from 1 to p do #6
if i <> line then
c:=M[i,j];
for jj from j to q do
M[i,jj]:=M[i,jj]-M[line,jj]*c
od
fi
od;
line:=line+1;
fi
od;
eval(M)
end: # echelon
n:=10:
H:=matrix(n,n,proc(i,j) 1/(i+j-1) end);
In:=array(1..n,1..n,identity);
M:=linalg[augment](H,In);
N:=echelon(M);
K:=linalg[submatrix](N,1..n,n+1..2*n);
evalm(H &* K);
n:=10:
A:=matrix(n,n,proc(i,j) i+j-1 end);
# fin probleme.
# probleme page 381.
gcdex(a,b,x,'u','v'):
pi[A]:=evalm(subs(x=L,collect(v*b,x))):
pi[B]:=evalm(subs(x=L,collect(u*a,x))):
linalg[iszero](evalm(pi[A]^2-pi[A])),
linalg[iszero](evalm(pi[B]^2-pi[B]));
linalg[iszero](evalm(pi[A]&*pi[B])),
linalg[iszero](evalm(pi[B]&*pi[A])),
linalg[iszero](evalm(pi[A]+pi[B]-&*()));
A:=evalm(subs(x=L,a)):
B:=evalm(subs(x=L,b)):
linalg[iszero](evalm(A &* pi[A])),
linalg[iszero](evalm(B &* pi[B]));
linalg[nullspace](A);
map(proc(z) evalm(pi[A]&*z-z) end,");
linalg[nullspace](B);
map(proc(z) evalm(pi[B]&*z-z) end,");
gcdex(a,b,x,'u','v'):
pi[A]:=polytodiffop(v*b,x,t):
pi[B]:=polytodiffop(u*a,x,t):
sol:=dsolve(P(z(t)),z(t));
eval(subs(sol,pi[A](z(t))));
eval(subs(sol,pi[B](z(t))));
a:=t^2*x^2+1:
b:=t*x-1:
p:=a*b;
A:=polytodiffop(a,x,t):
B:=polytodiffop(b,x,t):
P:=polytodiffop(p,x,t):
dsolve(P(z(t)),z(t));
dsolve(A(z(t)),z(t));
dsolve(B(z(t)),z(t));
# fin probleme.
# section 3, page 384.
# probleme page 384.
d:=6:
p:=3:q:=5:
A:=matrix(d,d):
R:=array(1..d,1..d,sparse):
for i to d do R[i,i]:=1 od:
R[p,p]:=cos(theta):
R[q,q]:=cos(theta):
R[p,q]:=sin(theta):
R[q,p]:=-sin(theta):
B:=evalm(transpose(R) &* A &* R):
B:=map(collect,map(combine,B,trig),{sin,cos});
B[p,q];
rotation:=proc(A,p,q)
local d,kappa,t,c,s,B,i,j;
d:=linalg[coldim](A); #1
kappa:=(A[q,q]-A[p,p])/2/A[p,q];
if kappa>=0 then t:=evalf(-kappa+sqrt(1+kappa^2))
else t:=evalf(-kappa-sqrt(1+kappa^2))
fi;
c:=evalf(1/sqrt(1+t^2));
s:=c*t;
B:=array(1..d,1..d,symmetric); #2
for i to d do for j to i do B[i,j]:=A[i,j] od;
for j to d do
B[p,j]:=c*A[p,j]-s*A[q,j];
B[q,j]:=s*A[p,j]+c*A[q,j];
od
od;
B[p,q]:=0;
B[p,p]:=A[p,p]-t*A[p,q];
B[q,q]:=A[q,q]+t*A[p,q];
eval(B)
end: # rotation
off:=proc(A)
local d,i,j;
d:=linalg[coldim](A);
add(add(evalf(A[i,((i+j-1) mod d)+1]^2),j=1..d-1),i=1..d)
end: # off
nicematrixplot:=proc(B,p,q,h)
plots[display](plots[matrixplot](map(abs,B),heights=histogram,
axes=frame,gap=0.25,style=patch),
plottools[sphere]([p+1/2,q+1/2,h],1/4,color=red),
scaling=constrained);
end: # nicematrixplot
visualjacobithreshold:=proc(A::matrix(realcons,square),
epsilon::realcons,maxsweepnb::posint)
local d,B,threshold,sweepnb,p,q,iternb,counter,pic;
d:=linalg[coldim](A); #1
if not (evalb(linalg['indexfunc'](A)='symmetric')
or convert([seq(seq(evalb(A[i,j]=A[j,i]),
j=i+1..d),i=1..d)],`and`))
then
ERROR(`expected a symmetric matrix, but received`,eval(A))
fi;
B:=map(evalf,A); #2
counter:=0:
pic[counter]:=nicematrixplot(B,-1,-1,0);
threshold:=evalf(sqrt(off(B,d)/d/(d-1)));
Digits:=Digits+5;
threshold:=evalf(sqrt(off(B,d)/d/(d-1)));
iternb:=0;
for sweepnb to maxsweepnb while threshold>=epsilon do #3
for p to d-1 do
for q from p+1 to d do
counter:=counter+1;
pic[counter]:=nicematrixplot(B,p,q,abs(B[p,q])+1);
if abs(B[p,q])>=threshold then
iternb:=iternb+1;
counter:=counter+1;
B:=rotation(B,p,q);
pic[counter]:=nicematrixplot(B,p,q,abs(B[p,q])+1);
fi
od
od;
threshold:=evalf(sqrt(off(B,d)/d/(d-1)));
od;
if threshold>=epsilon then #4
userinfo(2,'visualjacobithreshold',`process stopped
after`,sweepnb-1,`sweepings and`,iternb,`iterations`)
else userinfo(2,'visualjacobithreshold',`result obtained
after`,sweepnb-1,`sweepings and`,iternb,`iterations`)
fi;
plots[display]([seq(pic[i],i=0..counter)],
insequence=true,orientation=[-20,80]);
end: # visualjacobithreshold
d:=5:
A:=matrix(d,d,proc(i,j) d/(i+j-1) end):
infolevel[visualjacobithreshold]:=2:
visualjacobithreshold(A,10^(-1),10);
# fin probleme.
# probleme page 388.
gauss1:=proc(Q,x,Q1)
local QQ,a,b;
QQ:=collect(Q,x);
a:=coeff(QQ,x,2);
b:=coeff(QQ,x,1);
Q1:=expand(QQ-a*(x+b/2/a)^2);
a*(x+b/2/a)^2
end: # gauss1
gauss2:=proc(Q,x,y,Q1)
local QQ,a,b,c,u,v;
QQ:=collect(Q,{x,y});
a:=coeff(coeff(QQ,x,1),y,1);
b:=coeff(coeff(QQ,x,1),y,0);
c:=coeff(coeff(QQ,y,1),x,0);
u:=x+c/a;
v:=y+b/a;
Q1:=expand(QQ-a/4*(u+v)^2+a/4*(u-v)^2);
a/4*(u+v)^2-a/4*(u-v)^2
end: # gauss2
gauss:=proc(P)
local Q,ind,decomposition,path,n,counter,i,x,QQ,a,b,
result,j,y,c,u,v,k;
Q:=op(1,P);
if Q=0 then RETURN(P) fi;
ind:=op(2,P);
n:=nops(ind);
decomposition:=op(3,P);
path:=op(4,P);
counter:=0;
for i to n do
x:=ind[i];
QQ:=collect(Q,x,normal); #
if degree(QQ,x)=2 then
counter:=counter+1;
a:=coeff(QQ,x,2);
b:=coeff(QQ,x,1);
result[counter]:=[normal(QQ-a*(x+b/2/a)^2), #
[seq(ind[j],j=1..i-1),seq(ind[j],j=i+1..n)],
decomposition+a*(x+b/2/a)^2,[op(path),x^2]]
fi
od;
for i to n-1 do
for j from i+1 to n do
x:=ind[i];
y:=ind[j];
QQ:=collect(Q,{x,y},distributed,normal); #
if has(QQ,x*y) and degree(QQ,x)=1
and degree(QQ,y)=1 then
counter:=counter+1;
a:=coeff(coeff(QQ,x,1),y,1);
b:=coeff(coeff(QQ,x,1),y,0);
c:=coeff(coeff(QQ,y,1),x,0);
u:=x+c/a;
v:=y+b/a;
result[counter]:=[normal(QQ-a/4*(u+v)^2+a/4*(u-v)^2), #
[seq(ind[k],k=1..i-1),seq(ind[k],k=i+1..j-1),
seq(ind[k],k=j+1..n)],
decomposition+a/4*(u+v)^2-a/4*(u-v)^2,[op(path),x*y]]
fi
od
od;
seq(result[i],i=1..counter)
end: # gauss
Q:=x^2-y^2+y*z+z*t;
ind:=[x,y,z,t]:
P:=[Q,ind,0,[]]:
n:=nops(ind):
PP:=[P]:
for i to n while PP<>[] do
AAA:=map(gauss,PP);
QQ[i]:=select(proc(expr) op(1,expr)=0 end,AAA);
PP:=remove(proc(expr) op(1,expr)=0 end,AAA);
od:
map(proc(expr) [op(3,expr),op(4,expr)] end,
[seq(op(QQ[i]),i=1..n)]);
Q:=(1-cos(theta))*x[1]^2+2*sin(theta)*x[1]*x[2]+
(1+cos(theta))*x[2]^2-2*cos(theta)*x[1]*x[3]+
2*sin(theta)*x[2]*x[3]+sin(theta)*x[3]^2:
gauss1(Q,x[1],'Q1'):
collect(Q1,x[2],normal);
gaussall:= proc(P)
local PP,QQ,n,AAA,i;
n := nops(op(2,P));
PP := [P];
for i to n while PP <> [] do
AAA := map(gauss,PP);
QQ[i] := select(proc(expr) op(1,expr)=0 end, AAA);
PP := remove(proc(expr) op(1,expr)=0 end, AAA)
od;
map(proc(expr) [op(3,expr),op(4,expr)] end,
[seq(op(QQ[j]), j = 1 .. i - 1)])
end: # gaussall
Q:=(cos(theta)^2+sin(theta)^2-1)*x*y+y*z-z*x:
ind:=[x, y, z]:
gaussall([Q,ind,0,[]]);
Normalizer:=proc(expr) combine(normal(expr),trig) end:
gaussall([Q,ind,0,[]]);
# fin probleme.
# section 4, page 392.
# probleme page 392.
bernstein:=proc(n::posint,f,x::name)
local i;
add(binomial(n,i)*x^i*(1-x)^(n-i)*subs(x=i/n,f),i = 0 .. n)
end: # bernstein
a:=abs(x-1/2):
for i in {5,10,20} do Ba[i]:=bernstein(i,a,x) od:
plot([a,Ba[5],Ba[10],Ba[20]],x=0..1,
color=[black,red,green,blue]);
f:=piecewise(t<1/2,-4*(t-1/4),4*(t-3/4)):
g:=piecewise(t<1/4,4*t,t<3/4,-4*(t-1/2),4*(t-1)):
Bf:=bernstein(20,f,t):
Bg:=bernstein(20,g,t):
plot([[f,g,t=0..1],[Bf,Bg,t=0..1]],thickness=[2,1],
color=black,scaling=constrained);
f:=piecewise(x<1/2,0,x-1/2):
for k to 10 do B[k]:=evalf(subs(x=0.5,bernstein(10*k,f,x))) od:
seq(evalf(ln(B[k])/ln(10*k)),k=1..10);
seq(evalf(ln(B[2*k]/B[k])/ln(2)),k=1..5);
# fin probleme.
# probleme page 395.
beta:=2*log(2)-log(1+alpha^2)-2*alpha*arctan(1/alpha):
normal(diff(beta,alpha)):
alpha0:=fsolve(beta=0,alpha=0..infinity);
runge:=proc(alpha,m)
local abscissae,k,f,q;
f:=1/(x^2+alpha^2);
abscissae:=[seq((2*k+1)/(2*m),k=-m..m-1)];
q:=interp(abscissae,[seq(subs(x=k,f),k=abscissae)],x);
plot({f,q},x=-1..1);
end: # runge
runge(2/5,10);
# fin probleme.
# section 5, page 397.
# probleme page 397.
addI:=proc(I,J)
op(1,I)+op(1,J)..op(2,I)+op(2,J);
end: # addI
multiplyI:=proc(I,J)
local bounds;
bounds := op(1,I)*op(1,J), op(1,I)*op(2,J),
op(2,I)*op(1,J), op(2,I)*op(2,J);
min(bounds)..max(bounds);
end: # multiplyI
powerI:=proc(I,n)
local an,bn;
an:=op(1,I)^n;
bn:=op(2,I)^n;
if type(n,odd) then an..bn
else
if op(1,I)*op(2,I)<=0 then 0..max(an,bn)
else min(an,bn)..max(an,bn)
fi
fi
end: # powerI
imageI:=proc(f,x,I)
local k,image;
if op(0,f)=`*` then
image := imageI(op(1,f),x,I);
for k from 2 to nops(f) do
image := multiplyI(image, imageI(op(k,f),x,I));
od;
elif op(0,f)=`+` then
image := imageI(op(1,f),x,I);
for k from 2 to nops(f) do
image := addI(image, imageI(op(k,f),x,I));
od;
elif op(0,f)=`^` then
image := powerI(imageI(op(1,f),x,I),op(2,f));
elif f=x then
image := I;
elif type(f,rational) then
image := f..f;
else ERROR(`input expression is not polynomial in `,x)
fi;
image;
end: # imageI
for i to 4 do p:=orthopoly[T](i,x):
imageI(p,x,0..1);
imageI(convert(p,horner,x),x,0..1)
od;
easysolveI:=proc(f::polynom(rational),x::name,I::range,
epsilon::positive)
easyrsolveI(convert(f,horner),x,I,epsilon);
end: # easysolveI
easyrsolveI:=proc(f,x,I,epsilon)
local fI,middle;
fI:=imageI(f,x,I);
if (op(1,fI)>0 or op(2,fI)<0) then RETURN([]); fi;
if op(2,I)-op(1,I)0 or op(2,fI)<0) then RETURN([]) fi;
if op(2,I)-op(1,I)0 or op(2,DfI)<0) then #1
RETURN([newtonI(f,df,x,I,epsilon)])
fi;
middle:= (op(1,I)+op(2,I))/2;
[op(rsolveI(f,df,x,op(1,I)..middle,epsilon)),
op(rsolveI(f,df,x,middle..op(2,I),epsilon))]
end: # rsolveI
newtonI:=proc(f,df,x,I,epsilon)
local middle, fmiddle, DfI, bounds, boundmin, boundmax;
if op(2,I)-op(1,I)op(2,I)) then RETURN(NULL) fi;
boundmin:=max(boundmin,op(1,I));
boundmax:=min(boundmax,op(2,I));
newtonI(f,df,x,boundmin..boundmax,epsilon)
end: # newtonI
# fin probleme.
# probleme page 401.
crack:=proc(g,x)
local g1,g1da,`1/x/ln(x) part`,`non 1/x/ln(x) part`,
`1/x part`,`non (1/x/ln(x) or 1/x) part`;
g1:=normal(diff(g,x)/g); #1
g1da:=expand(asympt(g1,x));
if type(g1da,`+`) then `1/x/ln(x) part`:=select(has,g1da,ln(x))
elif has(g1da,ln(x)) then `1/x/ln(x) part`:=g1da #2
else `1/x/ln(x) part`:=0
fi;
`non 1/x/ln(x) part`:=g1da-`1/x/ln(x) part`; #3
if `non 1/x/ln(x) part`=0 then `1/x part`:=0 #4
elif type(`non 1/x/ln(x) part`,`+`) then
`1/x part`:=select(auxiliary1,`non 1/x/ln(x) part`,x)
elif limit(1/(`non 1/x/ln(x) part`*x),x=infinity)<>0 then
`1/x part`:= `non 1/x/ln(x) part`
else `1/x part`:=0
fi;
`non (1/x/ln(x) or 1/x) part`:= #5
g1da-`1/x/ln(x) part`-`1/x part`;
g1,`non (1/x/ln(x) or 1/x) part`,`1/x part`,`1/x/ln(x) part`
end: # crack
auxiliary1:=proc(z,x)
evalb(limit(1/(z*x),x=infinity)<>0)
end: # auxiliary1
integrationbyparts:=proc(g,x,g1,`non (1/x/ln(x) or 1/x) part`,
`1/x part`,`1/x/ln(x) part`)
local dominantterm,h,c,alpha,beta;
if `non (1/x/ln(x) or 1/x) part`<>0 then #1
h:=1/g1;
[g*h,expand(asympt(-g*diff(h,x),x))]
elif `1/x part`<>0 then
alpha:=normal(x*`1/x part`); #2
if `1/x/ln(x) part`<>0 then
beta:=normal(x*ln(x)*`1/x/ln(x) part`)
else
beta:=0
fi;
if alpha<>-1 then #3
[x*g/(alpha+1),expand(asympt(-beta/(alpha+1)*g/ln(x),x))]
else
[int(g,x),0]
fi
else #4
[x*g,expand(asympt(-x*diff(g,x),x))]
fi;
end: # integrationbyparts
intasympt:=proc(f,x,omega)
local F,G,k,g,prec,Prec,result,thingummy,todo,newPrec;
if nargs>2 then #1
Order:=omega
fi;
F:=expand(asympt(f,x)); #2
if type(F,`+`) then G:=[op(F)] else G:=[F] fi; #3
result:=0; #4
for k to Order while not has(G,O) do #5
for g in G do
thingummy:=integrationbyparts(g,x,crack(g,x));
result:=result+thingummy[1];
todo[g]:=thingummy[2]
od;
G:=[seq(flatten(todo[g]),g=G)];
od;
g:=select(has,G,O);
if g<>[] then #6
g:=op(g);
prec:=eval(subs(O=Identity,g));
Prec:=integrationbyparts(prec,x,crack(prec,x))[1]; #7
G:=remove(has,G,O);
while G<>[] do
for g in G do
if limit(prec/g,x=infinity)=0 then #8
thingummy:=integrationbyparts(g,x,crack(g,x));
result:=result+thingummy[1];
todo[g]:=thingummy[2]
else todo[g]:=0
fi
od:
G:=[seq(flatten(todo[g]),g=G)];
od;
else
if nops(G)>0 then #9
Prec:=integrationbyparts(G[1],x,crack(G[1],x))[1];
for k from 2 to nops(G) do
newPrec:=integrationbyparts(G[k],x,crack(G[k],x))[1];
if limit(Prec/newPrec,x=infinity)=0 then Prec:=newPrec fi
od
else #10
Prec:=0
fi
fi;
Prec:=expand(asympt(Prec,x)); #11
if type(Prec,`+`) then Prec:=remove(auxiliary2,Prec,Prec,x) fi;
Prec:=eval(subs(O=Identity,Prec));
result+O(1)*Prec #12
end: # intasympt
flatten:=proc(term)
if term=0 then NULL
elif type(term,`+`) then op(term)
else term fi
end: # flatten
Identity:=proc(z) z end:
auxiliary2:=proc(z,y,x)
evalb(limit(z/y,x=infinity)=0)
end: # auxiliary2
# fin probleme.
# section 6, page 405.
# probleme page 405.
kummer:=proc(F::ratpoly,p::posint,n::name)
local FF,N,D,d,u,i,G,alpha,j,delta,term;
FF:=normal(F,expanded);
N:=numer(FF); #1
D:=denom(FF);
d:=degree(N,n)-degree(D,n);
if d>-2 then
ERROR(`kummer expects its first argument to be the
term of a convergent series, but received`,F)
fi;
u:=0; #2
for i from -d to p-d-1 do
G:=alpha/mul(n+j,j=0..i-1);
delta:=normal(FF-G,expanded);
N:=numer(delta);
term:= subs(alpha=solve(lcoeff(N,n),alpha),G);
u:=u+term;
FF:=FF-term
od;
FF:=factor(FF); #3
u+FF
end: # kummer
F:=1/n-1/(n+x):
kummer(F,3,n):
map(factor,");
# fin probleme.
# probleme page 406.
n:=4:
z[0]:=1:z[1]:=1+I:z[2]:=-1+I:z[3]:=-I:z[4]:=1:
g:=piecewise(seq(op([t<=j/4,4*(z[j-1]*(j/4-t)+z[j]*(t-(j-1)/4))])
,j=1..4)):
z[-1]:=z[n-1]:
for j from -1 to n do t[j]:=j*2*Pi/n od:
for j from -1 to n-1 do v[j]:=(z[j+1]-z[j])/(t[j+1]-t[j]) od:
pic[-1]:=plots[complexplot](g,t=0..1):
c0:=add((z[j]+z[j+1])/2*(t[j+1]-t[j]),j=0..n-1)/2/Pi:
ck:=add((v[j-1]-v[j])*exp(-I*k*t[j]),j=0..n-1)/2/Pi/k^2:
nn:=5:
s[nn]:=add(subs(k=kk,ck)*exp(I*kk*t),kk=-nn..-1)+c0
+add(subs(k=kk,ck)*exp(I*kk*t),kk=1..nn):
pic[nn]:=plots[complexplot](s[nn],t=0..2*Pi):
plots[display]({pic[-1],pic[nn]},scaling=constrained);
zeta:=exp(2*I*Pi/n):
simplify(evalc(-(1-zeta)^2/zeta),trig);
n:=5:
ll:=4:
s:=(sin(Pi/n)/(Pi/n))^2
*add(subs(k=1+l*n,exp(I*k*t)/k^2),l=-ll..ll);
plots[complexplot](s,t=0..2*Pi,scaling=constrained);
n:=11:
ll:=4:
for r to 5 do
s:=add(subs(k=r+l*n,exp(I*k*t)/k^2),l=-ll..ll);
plots[complexplot](s,t=0..2*Pi,scaling=constrained);
od;
# fin probleme.
# section 7, page 409.
# probleme page 409.
VX:=X+Y+1:
VY:=X+X*Y:
criticallist:=[solve({VX,VY},{X,Y})];
Jac:=linalg[jacobian]([VX,VY],[X,Y]);
for i to nops(criticallist) do
J[i]:=subs(criticallist[i],eval(Jac));
V[i]:=linalg[eigenvects](J[i])
od:
seq({V[i]},i=1..nops(criticallist));
VX:=-Y-X*sqrt(X^2+Y^2):
VY:=X-Y*sqrt(X^2+Y^2):
criticallist:=[solve({VX,VY},{X,Y})]:
Jac:=linalg[jacobian]([VX,VY],[X,Y]);
sys:=subs(X=x(t),Y=y(t),{'diff'(X,t)=VX,'diff'(Y,t)=VY}):
vars:={x(t),y(t)}:
n:=6:
for i to nops(criticallist) do
inits:=subs(criticallist[i],[seq([x(0)=X+1/10*cos(2*k*Pi/n),
y(0)=Y+1/10*sin(2*k*Pi/n)],k=0..n-1)]);
trange:=-2..2;
pict[i]:=DEtools[phaseportrait](sys,vars,trange,inits,
color=black,linecolor=black,
axes=boxed,scaling=constrained)
od:
for i to nops(criticallist) do pict[i] od;
assume(r>0):
subs(X^2+Y^2=r^2,X=r*cos(theta),Y=r*sin(theta),eval(Jac)):
Jacpol:=map(normal,");
J[1]:=subs(r=0,eval(Jacpol));
linalg[eigenvects](J[1]);
# fin probleme.
# probleme page 413.
algtodiffeq:=proc(P,expr::function)
local x,y,Q,factorizedQ,d,A,B,U,V,G,elim,degy,T,i,k,sol,K;
y:=op(0,expr);
x:=op(1,expr);
if nops(expr)>1 or not type(x,name) then
ERROR(`algtodiffeq expects its second argument
to be a univariate function, but received`,expr)
fi;
Q:=subs(y(x)=Y,y=Y,x=X,P); #2
if not type(Q,polynom(ratpoly(rational,X),Y)) then
ERROR(`algtodiffeq expects its first argument to be
a polynomial with respect to`,expr,`whose coefficients are
rational fractions with respect to`,x,`but received`,P)
fi;
factorizedQ:=factor(Q);
if type(factorizedQ,`*`) and
nops(select(has,{op(factorizedQ)},Y))>1 then
ERROR(`algtodiffeq expects its first argument to be an
irreducible polynomial, but received`,
subs(Y=y(x),X=x,factor(Q)))
fi;
d:=degree(Q,Y); #3
A:=-diff(Q,X); B:=diff(Q,Y);
G:=gcdex(B,Q,Y,'U','V'); #4
T[1]:=collect(rem(A*U,Q,Y)/G,Y,normal);
if d=2 then #5
sol:=diff(y(x),x)*subs(X=x,Y=y(x),denom(T[1]))
-subs(X=x,Y=y(x),numer(T[1]));
RETURN(collect(sol,[diff(y(x),x),y(x)]))
fi;
elim:=matrix(d-2,1,[seq(coeff(T[1],Y,i),i=2..d-1)]); #6
K:=linalg[kernel](elim);
for k from 2 while nops(K)=0 do
T[k]:=collect(rem(T[1]*diff(T[k-1],Y)+diff(T[k-1],X),Q,Y),
Y,normal);
elim:=linalg[concat](elim,
vector([seq(coeff(T[k],Y,i),i=2..d-1)]));
K:=linalg[kernel](elim);
od;
K:=op(1,K); #7
sol:=subs(X=x,Y=y(x),add(K[i]*(diff(y(x),x$i)-T[i]),i=1..k-1));
collect(numer(normal(sol)),
[seq(diff(y(x),x$(d-i)),i=1..d-1),y(x)]);
end: # algtodiffeq
algtodiffeq((x^2-1)*y(x)^2-1,y(x));
dsolve(",y(x));
algtodiffeq(y-1-x*y^3,y(x));
with(share):
readshare(gfun,analysis):
with(gfun):
?gfun[algeqtodiffeq]
algeqtodiffeq(y-1-x*y^3,y(x));
# fin probleme.
# probleme page 415.
lindeq:=alpha(x)*diff(w(x),x,x)
+beta(x)*diff(w(x),x)+gamma(x)*w(x):#
x0:=x0: #
w0:=w0: #
w1:=w1: #
lininitcond:=w(x0)=w0,D(w)(x0)=w1: #
if type(lindeq,`=`) then lindeq:=op(1,lindeq)-op(2,lindeq) fi;
eval(subs(w(x)=exp(-int(y(t),t=x0..x)),lindeq)):
ricdeq:=select(has,factor("),diff);
ricinitcond:=y(x0)=subs(lininitcond,-D(w)(x0)/w(x0));
ricdeq:=diff(y(x),x)=a(x)*y(x)^2+b(x)*y(x)+c(x); #
x0:=x0; #
y0:=y0; #
ricinitcond:=y(x0)=y0; #
if type(ricdeq,`=`) then ricdeq:=op(1,ricdeq)-op(2,ricdeq) fi;
redricdeq:=subs(diff(y(x),x)=Y1,y(x)=Y,ricdeq);
redricdeq:=collect(redricdeq,{Y,Y1});
aa:=-coeff(redricdeq,Y,2)/coeff(redricdeq,Y1,1);
eval(subs(y(x)=-diff(w(x),x)/w(x)/aa,ricdeq));
lindeq:=numer(normal(w(x)*"))
lininitcond:=w(x0)=1,subs(ricinitcond,D(w)(x0)=-y(x0));
dsolve({diff(y(x),x)=1/x^4-y(x)^2,y(1)=eta},y(x)):
Y:=collect(normal(convert(subs(",y(x)),exp)),exp):
YY:=subs(exp(1/x)=E,Y);
ricdeq:={diff(y(x),x)=y(x)^2+x,y(0)=1}:
lindeq:=rictolin(ricdeq);
dsolve(lindeq,w(x));
W:=subs(",w(x));
WS:=map(simplify,series(W,x,20),GAMMA);
with(share):
readshare(gfun,analysis):
with(gfun):
seriestorec(WS,u(n));
WP:=convert(WS,polynom):
plot(WP,x=-1..5);
x_max:=fsolve(W,x=0.9..1);
dsolve(ricdeq,y(x));
Y:=subs(",y(x)):
YS:=map(simplify,series(Y,x,20),GAMMA):
YP:=convert(YS,polynom):
eps:=10^(-2):
plot(YP,x=-2..x_max-eps,view=[-2..1,-10..10]);
# fin probleme.
# section 8, page 420.
# probleme page 420.
exclusiontest:=proc(eqn,x,y,r)
local p;
p:=collect(eqn,{x,y},abs); #1
2*coeff(coeff(p,x,0),y,0)-subs(x=r,y=r,p);
end: # exclusiontest
dichotomy:=proc(eqn,x,y,Square,eps)
local p,d1,d2,d3;
p:=expand(subs(x=x+Square[1],y=y+Square[2],eqn));
if exclusiontest(p,x,y,Square[3])>0 then RETURN() fi; #2
if Square[3]<=eps then RETURN(Square) fi;
d3:=Square[3]*.5; d1:=Square[1]-d3; d2:=Square[2]-d3;
dichotomy(eqn,x,y,[d1,d2,d3],eps),
dichotomy(eqn,x,y,[d1+Square[3],d2,d3],eps),
dichotomy(eqn,x,y,[d1+Square[3],d2+Square[3],d3],eps),
dichotomy(eqn,x,y,[d1,d2+Square[3],d3],eps);
end: # dichotomy
do_plot:=proc(l)
PLOT(POLYGONS(op(map(squaretransform,l)),STYLE(LINE)))
end: # do_plot
squaretransform:=proc(l)
[ [l[1]-l[3],l[2]-l[3]], [l[1]+l[3],l[2]-l[3]],
[l[1]+l[3],l[2]+l[3]], [l[1]-l[3],l[2]+l[3]] ];
end: # squaretransform
exclusionplot:=proc(eqn::polynom,x,y,R::{numeric,positive},
eps::{numeric,positive})
do_plot([dichotomy(eqn,x,y,[0,0,R],eps)]);
end: # exclusionplot
exclusionplot(p2,x,y,2.5,0.01);
exclusionplot(p4,x,y,5,0.05);
# fin probleme.
# probleme page 421.
JP:=linalg[jacobian](P,[u,v]);
tu:=linalg[col](JP,1):
tv:=linalg[col](JP,2):
E:=normalizer(linalg[dotprod](tu,tu)):
F:=normalizer(linalg[dotprod](tu,tv)):
G:=normalizer(linalg[dotprod](tv,tv)):
metric:=matrix(2,2,[E,F,F,G]);
P:=[R*sin(u)*cos(v),R*sin(u)*sin(v),R*cos(u)]:
domain:=u=0..Pi,v=-Pi..Pi:
normalizer:=readlib(`combine/trig`):
illustration:=R=1: #a
JP:=array(1..3,1..2,[seq([seq(diff(P[i],j),j=[u,v])],i=1..3)]):
tu:=linalg[col](JP,1):
tv:=linalg[col](JP,2):
E:=normalizer(linalg[dotprod](tu,tu)):
F:=normalizer(linalg[dotprod](tu,tv)):
G:=normalizer(linalg[dotprod](tv,tv)):
metric:=matrix(2,2,[E,F,F,G]); #z
arc:=[Pi/2*(2-t),t*Pi/2]:
arcdomain:=t=0..1:
Jarc:=vector(2,[seq(diff(arc[i],t),i=1..2)]):
dlength2:=evalm(transpose(Jarc)
&*subs(u=arc[1],v=arc[2],eval(metric))&*Jarc);
dlength:=normalizer(dlength2^(1/2));
evalf(subs(illustration,Int(dlength,arcdomain)));
P:=[(a+r*cos(v))*cos(u),(a+r*cos(v))*sin(u),r*sin(v)];
domain:=u=-Pi..Pi,v=-Pi..Pi;
normalizer:=readlib(`combine/trig`): #a
#z
darea2:=normalizer(linalg[det](metric)):
area:=student[Doubleint](sqrt(darea2),domain);
value(area);
normalvector:=map(normalizer,linalg[crossprod](tu,tv)):
darea2:=map(normalizer,
linalg[dotprod](normalvector,normalvector));
darea:=simplify(factor(expand(darea2))^(1/2),power,symbolic);
N:=evalm(map(expand,normalvector)/darea):
N:=map(normal,N):
N:=map(readlib(`simplify/trig`),N);
JN:=linalg[jacobian](N,[u,v]);
Ntu:=linalg[col](JN,1):
Ntv:=linalg[col](JN,2):
NE:=normalizer(linalg[dotprod](Ntu,Ntu)):
NF:=normalizer(linalg[dotprod](Ntu,Ntv)):
NG:=normalizer(linalg[dotprod](Ntv,Ntv)):
Nmetric:=matrix(2,2,[NE,NF,NF,NG]);
Ndarea2:=normalizer(linalg[det](Nmetric)):
Gcurvature2:=Ndarea2/darea2:
factor(expand(Gcurvature2)):
Gcurvature:=simplify("^(1/2),power,symbolic);
PP:=subs(illustration,P):
surface:=plot3d(PP,domain,scaling=constrained,
color=blue,style=wireframe):
NN:=subs(illustration,N):
Nsphere:=plot3d(NN,domain,scaling=constrained,
color=blue,style=wireframe):
eps:=0.1:
u0:=0:v0:=3*Pi/8:
littledomain:=u=u0-eps..u0+eps,v=v0-eps..v0+eps:
pieceofsurface:=plot3d(PP,littledomain,scaling=constrained,
color=red,style=patch):
plots[display]({surface,pieceofsurface},orientation=[45,45]);
pieceofNsphere:=plot3d(NN,littledomain,scaling=constrained,
color=red,style=patch):
plots[display]({Nsphere,pieceofNsphere},orientation=[45,45]);
ddarea:=subs(illustration,sqrt(darea2));
Nddarea:=subs(illustration,sqrt(Ndarea2));
littlearea:=value(student[Doubleint](ddarea,littledomain));
Nlittlearea:=value(student[Doubleint](Nddarea,littledomain));
Nlittlearea/littlearea,
evalf(subs(illustration,subs(u=u0,v=v0,Gcurvature)));
# fin probleme.