Maple
Son bon usage en mathématiques
Philippe Dumas, Xavier Gourdon
Springer-Verlag, 1997
Code Maple du Chapitre 7, Livre second.
# section 1, page 297.
a:=Y/sqrt(1+Y^2):
b:=sin(X):
deq:=subs(Y=y(x),a)*diff(y(x),x)=subs(X=x,b);
deqsol:=dsolve(deq,y(x));
constantname:=op(indets(deqsol,name) minus indets(deq)):
Sol:=[solve(deqsol,y(x))];
kmin:=-5:
kmax:=5:
xmin:=-5:
xmax:=5:
for i to nops(Sol) do
picture[i]:=plot({seq(subs(constantname=k,Sol[i]),
k=kmin..kmax)},x=xmin..xmax,color=black);
od:
plots[display]({seq(picture[i],i=1..nops(Sol))});
deq:=x*diff(y(x),x)-y(x)=(x+y(x))^2;
DEtools[Dchangevar](y(x)=C*x,deq,x);
newdeq:=DEtools[Dchangevar](y(x)=t(x)*x,deq,x);
dsolve(newdeq,t(x)):
T:=solve(",t(x));
eval(subs(y(x)=t(x)*x,deq));
polardeq:=DEtools[Dchangevar]({x=rho(theta)*cos(theta),
y(x)=rho(theta)*sin(theta)},
deq,x,theta);
polarsol:=dsolve(polardeq,rho(theta)):
Rho:=subs(polarsol,rho(theta));
lambda:=1:
deq:=sin(x)*diff(y(x),x)-lambda*y(x)=x:
dsolve(deq,y(x)):
Y:=subs(",y(x)):
Y:=collect(Y,ln,factor);
mirror:=proc(expr::algebraic,eps::name)
local body,bodyseries,tc;
if type(expr,function) and op(0,expr)=ln
then body:=op(1,expr);
bodyseries:=series(body,eps);
tc:=op(1,bodyseries);
if tc<0 then ln(-body) else expr fi
elif has(expr,ln) then map(mirror,expr,eps)
else expr
fi
end: # mirror
assume(k,integer):
x0:=2*k*Pi:
rY:=eval(subs(x=x0+epsilon,Y)):
rYseries:=series(rY,epsilon,4);
lY:=mirror(eval(subs(x=x0-epsilon,Y)),epsilon):
lYseries:=series(lY,epsilon,4);
x0:=(2*k+1)*Pi:
rY:=mirror(eval(subs(x=x0+epsilon,Y)),epsilon):
rYseries:=series(rY,epsilon,6);
lY:=eval(subs(x=x0-epsilon,Y)):
lYseries:=series(lY,epsilon,6);
YY:=eval(subs(ln=proc(z) ln(abs(z)) end,Y)):
constantname:=op(indets(Y,name) minus indets(deq)):
nmin:=-2:
nmax:=3:
window:=[nmin*Pi..(nmax+1)*Pi,-20..20]:
cset:={seq(cc/2,cc=-3..3)}:
for n from nmin to nmax do
picture[n]:=plot({seq(subs(constantname=c,YY),c=cset)},
x=n*Pi..(n+1)*Pi,view=window)
od:
plots[display]({seq(picture[n],n=nmin..nmax)});
sys:={diff(x(t),t)=y(t)+x(t)-1,
diff(y(t),t)=x(t)+x(t)*y(t)};
vars:={x(t),y(t)}:
xmin:=-2: xmax:=3: ymin:=-2: ymax:=2:
xstep:=0.5: ystep:=0.5:
inits:=[seq(seq([x(0)=k*xstep,y(0)=l*ystep],
k=ceil(xmin/xstep)..floor(xmax/xstep)),
l=ceil(ymin/ystep)..floor(ymax/ystep))]:
trange:=-1..1:
DEtools[DEplot](sys,vars,trange,inits,stepsize=0.2,color=black,
linecolor=black,axes=boxed,view=[xmin..xmax,ymin..ymax],
arrows=NONE,scaling=constrained);
DEtools[dfieldplot](sys,vars,t=trange,x=xmin..xmax,y=ymin..ymax,
color=black, scaling=constrained);
A[I]:=array(1..2,1..2,[[lambda,0],[0,lambda]]); #1
A[II]:=array(1..2,1..2,[[lambda,1],[0,lambda]]);
A[III]:=array(1..2,1..2,[[lambda,0],[0,mu]]);
A[IV]:=array(1..2,1..2,[[lambda,0],[0,mu]]);
A[V]:=array(1..2,1..2,[[alpha,-beta],[beta,alpha]]);
A[VI]:=array(1..2,1..2,[[0,-beta],[beta,0]]);
cases:=[I,II,III,IV,V,VI];
for i in cases do
sys[i]:={diff(x(t),t)=A[i][1,1]*x(t)+A[i][1,2]*y(t),
diff(y(t),t)=A[i][2,1]*x(t)+A[i][2,2]*y(t)}
od:
vars:={x(t),y(t)};
init:={x(0)=x0,y(0)=y0};
for i in cases do #2
S[i]:=dsolve(sys[i] union init,vars) od;
n:=12;
trange:=-1..1;
nber[I]:=2; #3
SUBS[I,1]:={lambda=1}; SUBS[I,2]:={lambda=-1};
nber[II]:=2;
SUBS[II,1]:={lambda=-2}; SUBS[II,2]:={lambda=2};
nber[III]:=2;
SUBS[III,1]:={lambda=-2,mu=-1}; SUBS[III,2]:={lambda=1,mu=2};
nber[IV]:=1;
SUBS[IV,1]:={lambda=-1,mu=2};
nber[V]:=2;
SUBS[V,1]:={alpha=-1,beta=1}; SUBS[V,2]:={alpha=1,beta=1};
nber[VI]:=2;
SUBS[VI,1]:={beta=1}; SUBS[VI,2]:={beta=-1};
printlevel:=2: #4
for i in cases do for j to nber[i] do
X:=subs(S[i],SUBS[i,j],x(t));
Y:=subs(S[i],SUBS[i,j],y(t));
plot([seq(subs(x0=cos(2*k*Pi/n),y0=sin(2*k*Pi/n),
[X,Y,t=trange]),k=0..n-1)], color=black,thickness=3,
axes=boxed,title=cat(`case `,convert(i,string),` `,
convert(SUBS[i,j],string)));
od od;
Besseldeq:=x^2*diff(w(x),x,x)+x*diff(w(x),x)+(x^2-nu^2)*w(x):
dsolve(Besseldeq,w(x),type=series);
normal(subs(w(x)=x^alpha*(1+add(c(n)*x^n,n=1..Order)),
Besseldeq)/x^alpha,expanded):
series(",x):
map(collect,",c);
Bessel0deq:=subs(nu=0,Besseldeq);
dsolve(Bessel0deq,w(x),type=series);
# section 2, page 312.
# exercice 7.7, page 314.
J:=linalg[JordanBlock](lambda,5);
matseries(J,cos(alpha*x),x);
# fin exercice 7.7.
# section 4, page 320.
# probleme page 320.
P:=Y-1-X*Y^3:
A:=-diff(P,X):
B:=diff(P,Y):
gcdex(B,P,Y,'U','V'):
G:=collect(U*B+V*P,Y,normal):
normal(rem(A*U/G,P,Y));
algtodiffeq(1+y+x*y^2,y(x));
# fin probleme.
# probleme page 323.
dsolve({diff(diff(w(x),x),x)-diff(w(x),x)+w(x),
w(0)=1,D(w)(0)=-eta},w(x));
Y:=collect(normal(eval(subs(",-diff(w(x),x)/w(x)))),
{sin,cos},factor);
# fin probleme.
# section 5, page 325.
# exercice 7.2, page 326.
deq:=diff(y(x),x)=(3*x+4*y(x)+5)/(2*x+3*y(x)-6):
newdeq:=DEtools[Dchangevar]({x=t+a,y(x)=w(t)+b},deq,x,t);
rightterm:=op(remove(has,[op(newdeq)],diff));
equ1:=remove(has,numer(rightterm),{t,w(t)});
equ2:=remove(has,denom(rightterm),{t,w(t)});
sol:=solve({equ1,equ2},{a,b});
dsolve(subs(sol,newdeq),w(t));
# fin exercice 7.2.
# exercice 7.5, page 327.
deq0:=(9*x+4)*x^2*diff(w(x),x,x)+(2*x+1)*w(x):
deqi:=DEtools[Dchangevar]({x=1/t,w(x)=y(t)},deq0,x,t);
Order:=5:
dsolve(deq0,w(x),type=series);
dsolve(deqi,y(t),type=series);
# fin exercice 7.5.
# exercice 7.6, page 328.
A:=matrix(2,2,[1,1,0,0]):
B:=matrix(2,2,[0,0,0,1]):
expA:=linalg[exponential](A):
expB:=linalg[exponential](B):
evalm(expA &* expB),evalm(expB &* expA),linalg[exponential](A+B);
for i in cases do
sys[i]:={diff(x(t),t)=A[i][1,1]*x(t)+A[i][1,2]*y(t),
diff(y(t),t)=A[i][2,1]*x(t)+A[i][2,2]*y(t)};
S[i]:=dsolve(sys[i] union init,vars);
Phi[i]:=array(1..2,1..2);
X:=subs(S[i],x(t));
Y:=subs(S[i],y(t));
Phi[i][1,1]:=coeff(X,x0);
Phi[i][1,2]:=coeff(X,y0);
Phi[i][2,1]:=coeff(Y,x0);
Phi[i][2,2]:=coeff(Y,y0);
od:
seq(eval(Phi[i]),i=cases);
d:=10:
r:=10:
haphazard:=rand(-r..r):
A:=matrix(d,d,[seq(1.*haphazard(),i=1..d^2)]):
tr:=linalg[trace](A):
for i to d do A[i,i]:=A[i,i]-tr/d od:
Digits:=20:
E:=linalg[exponential](A):
linalg[det](E);
# fin exercice 7.6.
# exercice 7.7, page 330.
matseries:=proc(A,t,x)
local s,n,p;
s:=series(t,x);
p:=convert(convert(s,polynom),horner,x);
evalm(subs(x=A,p))
end:
readlib(`linalg/matfunc`)(J,cos(alpha*x),x);
Order:=20:
Digits:=5:
A:=matrix(4,4,0.1);
matseries(A,1/(1+x)^(1/2),x);
evalm("^(-2)-1);
# fin exercice 7.7.
# exercice 7.8, page 332.
A:=array(1..4,1..4,sparse):
s2:=2^(1/2):
s3:=3^(1/3):
s6:=6^(1/2):
A[1,1]:=s6/4-s2/4:
A[1,2]:=3*s6/20+3*s2/20:
A[1,4]:=s6/5+s2/5:
A[2,1]:=-3*s6/20-3*s2/20:
A[2,2]:=16/25+9*s6/100-9*s2/100:
A[2,4]:=3*s6/25-3*s2/25-12/25:
A[3,3]:=1:
A[4,1]:=-s6/5-s2/5:
A[4,2]:=3*s6/25-3*s2/25-12/25:
A[4,4]:=9/25+4*s6/25-4*s2/25:
evalm((A-1)&*((1-t)+t*A)^(-1)): #
map(normal,"): #
L:=map(int,",t=0..1): #
L:=map(radnormal,eval(L));
evalf(arctan(3^(1/2)*2^(1/2)+3^(1/2)-2^(1/2)-2)/Pi):
convert(",confrac,red);
red;
readlib(`linalg/matfunc`)(A,ln(x),x);
A:=matrix(4,4,[3,1,-5,-1,3,1,1,5,-3,-3,-3,3,-3,5,-1,1]):
A:=evalm(1/6*A):
evalf(Eigenvals(A));
evalm((A-1)&*((1-t)+t*A)^(-1)):
map(normal,"):
L:=map(int,",t=0..1);
Gamma:=linalg[exponential](L,t);
Proj:=array(1..3,1..4,sparse):
for i to 3 do Proj[i,i]:=1 od:
ps:=combinat[powerset](4) minus {{}}:
for ss in ps do
for e in ss do
ss1:=ss minus {e};
if ss1={} then spt:=[0,0,0]
else spt:=evalm((Proj&*add(linalg[col](Gamma,f),f=ss1)))
fi;
ept:=evalm(spt+Proj&*linalg[col](Gamma,e));
segment[ss,e]:=[convert(spt,list),convert(ept,list)];
od
od;
cube:={seq(seq(segment[ss,e],e=ss),ss=ps)}:
fnb:=10:
for n from 0 to fnb do
pic[n]:=plots[spacecurve](subs(t=n/fnb,cube),color=black)
od:
plots[display]([seq(pic[n],n=0..fnb)],scaling=constrained,
insequence=true);
# fin exercice 7.8.