-
Structure des objets Maple.
Nous décrivons, aux pages 10-11, les
objets Maple comme des arbres. Cette description, suffisante en
premier approche, est fausse : les objets Maple sont des
DAG (Directed Acyclic Graph). La raison d'être de
cette structure est fort bien expliquée pages 68-69 de
Calcul formel : Mode d'emploi - Exemples en Maple,
Bruno Salvy, Claude Gomez, Paul Zimmermann, Masson ,1995. On en voit
aussi l'intérêt dans la section 2.2 de
Computing the Dimension of a Projective Variety: the
Projective Noether Maple Package,
(RR
3224), 1997, Marc Giusti, Klemens Hägele, Grégoire
Lecerf, Joël Marchand, Bruno Salvy
(référence Bibtex,
version
HTML). La procédure de substitution proposée
s'appuie sur une option remember; si un sous-arbre apparaît deux
fois dans l'expression considérée la descente dans le
sous-arbre ne sera effectuée qu'une fois; on utilise ainsi la
structure de DAG, puisqu'une sous-expression n'est traitée qu'une
fois et non pas autant de fois qu'elle apparaît dans l'expression.
-
Représentation arborescente.
Pour illustrer la structure des objets Maple, nous avons
dessiné plusieurs arbres (pages 11, 32,
51). Ceci a été fait en employant une
procédure dont nous fournissons le
code
et un fichier
test.
Cette procédure permet de se familiariser avec la structure
arborescente des expressions. Son fonctionnement est relativement
satisfaisant mais elle ne peut évidemment pas traiter de
grosses expressions. De plus le fait que le logiciel remplisse
toujours les images au mieux pose problème. Nous
espérons cependant que ce gadget vous amusera
(
feuille de travail V.4
).
-
Promenade sans collision.
Nous demandons de produire les promenades sans collision, pages 27-28. Un référence plus
complète que celle qui est donnée dans la solution,
pages 57-58, est le livre de Madras et Slade
[MaSl93].
-
Démontage des expressions Maple.
À la page 33, exercice 1.26, nous
demandons d'écrire une
procédure qui fait passer d'une expression Maple à sa
version Lisp. La solution fournie page 61
peut être
améliorée en une
procédure qui permet de traiter
une classe un peu plus large d'exemples. Rappelons que args est,
à l'intérieur d'une procédure, la séquence
des arguments passés à la procédure.
-
Itération sur les pentagones.
Pages 80-81, nous considérons une
transformation sur les pentagones qui est prétexte à
écrire une procédure. Après quelques
expériences, une conjecture émerge, mais un peu de
réflexion permet de montrer à la page 100 que la
conjecture est fausse. Autrement dit, une centaine d'essais est
démentie par une seule phrase. Pour comprendre le
phénomène un peu de
théorie
est nécessaire.
-
Préfixes des puissances de 2.
À la page 83, nous nous
intéressons aux préfixes de
l'écriture décimale des puissances de 2. Dans le livre
de G. Rauzy, Propriétés statistiques de suites
arithmétiques (
référence Bibtex), pages 13 et 14, on trouve une
explication du fait que parmi les N premières puissances de 2,
écrites en base 10, le chiffre 1 est celui qui apparaît
le plus fréquemment en tête (environ 0.30 N fois) alors
que le chiffre 9 est le plus rare à figurer en tête
(environ 0.04 N fois). Une petite
expérience le montre
aisément. Le raisonnement se généralise
à tout entier naturel n; si la longueur de son écriture
en base dix est L+1, il suffit de le considérer comme un chiffre
dans l'écriture en base 10^L pour conclure que parmi les N
premières puissances de 2 l'entier n apparaît environ N
(ln(n+1)-ln(n))/ ln(10^L) fois.
-
Préfixes des puissances de 2.
La phrase de la page 83 "On notera qu'il est possible
d'être plus efficace en inversant le problème." est un
peu obscure (volontairement). Elle a le sens suivant : au lieu de
chercher chaque motif dans les préfixes, on calcule tous
les préfixes et on garde en mémoire la première
occurrence de chaque motif. C'est
ainsi qu'a été
fabriqué le dessin de couverture.
-
Emploi d'hypothèses. L'exercice 2.17 de
la page 95 fait utiliser
assume. Voici un autre exemple d'emploi de cette
procédure.
> int(x^alpha,x=0..1);
(alpha + 1)
x 1
lim - ------------ + ---------
x -> 0+ alpha + 1 alpha + 1
> assume(alpha>-1):
> int(x^alpha,x=0..1);
1
----------
alpha~ + 1
Sur des cas simples, cette procédure est agréable, mais
dans un calcul compliqué elle devient une
boîte noire sur laquelle l'utilisateur n'a aucun
contrôle. De plus les hypothèses faites ne sont pas
prises en compte dans tout le logiciel. Nous préférons donc
éviter son emploi.
-
Développement décimal d'un
rationnel.
La procédure nicedecimalexpansion fournie page 371, en
réponse à la question 1.b de la page
125 a un coût en O(t) comme demandé. La version naïve, qui comporte
deux boucles imbriquées, a par contre un coût en O(t^2).
-
Calcul de la fonction pi.
Nous avons signalé à la page
137 l'existence d'un
algorithme de calcul de la fonction pi dont la complexité
en temps est (grosso modo) en x^(1/3). Cette méthode est due
à Lagarias, Miller et
Odlyzko
( référence
Bibtex). On en trouve une
présentation dans les résumés du Séminaire Algorithmes,
année 1995-1996, section Symbolic Computation,
Le calcul de grandes
valeurs de la fonction $\pi(x)$. Nous proposons une
implémentation de cet algorithme en Maple. Précisons
qu'elle n'est qu'un complément de l'article de Lagarias, Miller
et Odlyzko. Il est nécessaire de lire cet article pour
comprendre la séquence de calcul fournie. Cette
implémentation ne prétend absolument pas être
efficace. Son seul but est d'aider à la compréhension de
l'article.
On trouvera dans
[MoNi95,
problème 5]
un traitement de l'algorithme de Meissel et Lehmer qui est
intermédiaire entre la méthode élémentaire
que nous avons présentée et la méthode de
Lagarias, Miller et Odlyzko.
-
Syntaxe de rsolve.
On notera, pages 148-149, la variation
subtile qui apparaît dans l'écriture de la solution,
suivant que l'on entoure ou que l'on n'entoure pas la suite inconnue
d'accolades. Nous laissons le lecteur effectuer la comparaison par
lui-même.
-
Technique de Soeur Céline.
Page 160, nous présentons la technique
de Soeur Céline
dont le but est l'évaluation de sommes. Cette technique est le
point de départ d'un sujet en plein développement et qui
est fortement illustré dans les résumés du Séminaire Algorithmes,
section Symbolic Computation. On y trouve pour l'année
1992-1993 Symbolic
Computation with P-finite Sequences,
par Marko Petkovsek; pour l'année
1994-1995 Holonomic Systems
and Automatic Proofs of Identities, par Frédéric
Chyzak; pour l'année
1995-1996 Creative
Telescoping and Applications, par
Frédéric Chyzak; pour l'année
1996-1997 New Algorithms
for Definite Summation and Integration, par Frédéric
Chyzak. On l'aura
compris, Frédéric
Chyzak est le spécialiste du Projet Algorithme pour ce qui
concerne les algorithmes de preuve d'identités.
D'autre part Doron Zeilbeger
propose, entre autres, une présentation du livre A=B qu'il a
écrit avec Marko Petkovsek et Herb Wilf sur les algorithmes de
sommation et des programmes Maple. Vous trouverez aussi dans sa page
un portrait de la Soeur Céline.
-
Homographies.
Page 171, nous introduisons les homographies
de la sphère de Riemann et nous faisons allusion à une
définition plus élégante que celle que nous
donnons. Il s'agit de la notion d'espace projectif. Dans un espace
projectif, les droites parallèles se rencontrent à
l'infini, ce qui simplifie bien la vie du géomètre. On
trouvera une présentation des espaces projectifs dans le
plaisant livre de Pierre Samuel [Samuel86] ou
dans le livre, plus abordable par le débutant, de Daniel Lehmann et
Rudolphe Bkouche [LeBk88].
-
Groupe de transformations orthogonales
tridimensionnelles. Page 180, dans
l'exercice 3.3, on étend value et on introduit
value/&*. Tout cela fonctionne très bien, mais il
aurait été plus normal de ne pas définir cette
procédure. Cela a été rendu nécessaire par
le fait que value ne traite pas les arguments de f si
value/f n'est pas défini. Il aurait cependant
été plus clair de définir value/&* par
`value/&*`:=proc(expr) map(value,expr) end:
et de traiter l'exemple qui suit par
evalm(value(Rotation3d([0,0,1],Pi/3)&*Reflection3d([1,0,0])));
Dans la version V.5, value descend récursivement dans
les arguments de f même si value/f n'est pas
défini et il devient inutile de définir une
procédure value/&*.
-
Groupe de transformations orthogonales
tridimensionnelles. Page 181,
on fait étudier un groupe de transformations orthogonales
défini par un système de générateurs.
On obtient un autre exemple plaisant avec les instructions suivantes,
à employer pour modifier la correction des pages 189-190.
R:=rotation3d([0,0,1],Pi/2):
S:=rotation3d([1,0,0],Pi/2):
G:=evaln(G):
G[1]:=eval(R):
G[2]:=eval(S):
globalcounter:=2:
todocounter[1]:=2:
On engendre ainsi le groupe des rotations du cube, qui comporte
vingt-quatre éléments.
-
Méthode de Newton et séries
formelles.
Pages 201-202, exercice 4.7, nous
introduisons l'anneau des
séries formelles et nous montrons comment calculer un inverse
par la méthode de Newton. La technique employée a une
portée plus générale, comme on pourra le
constater en consultant dans les résumés du Séminaire Algorithmes,
section Symbolic Computation, année 1995-1996,
l'exposé de Bruno Salvy Linear Recurrences, Linear
Differential Equations and Fast Computation.
-
Solutions asymptotiques.
À la page 209, nous proposons un
problème qui porte sur
la résolution approchée (au sens du calcul asymptotique)
d'équations. On trouve dans les résumés du Séminaire Algorithmes,
année 1996-1997, section Symbolic Computation, un
exposé de Bruno Salvy, Asymptotique des fonctions
implicites et calcul formel, qui donne une vision plus large de la
question.
-
Solutions asymptotiques.
Pages 210-211, nous étudions une
équation différentielle avec retard qui fait
considérer l'équation s=exp(-s). Le calcul,
fortement guidé par l'énoncé est le suivant.
solve(s=exp(-s),s);
for n from -2 to 2 do LambertW(n,1.) od;n:='n':
On effectue un petit dessin.
nmax:=10:
plots[complexplot]([seq(LambertW(n,1.),n=-nmax..nmax)],
style=point,symbol=circle);
On pose s=-x+I*t et t=2*n*Pi-Pi/2-u.
L'étude du module permet de lier x, t et
z=exp(-x).
t:=sqrt(1/z^2-x^2);
En séparant parties réelle et imaginaire on a deux
équations.
equ[1]:=x*z=sin(u);
equ[2]:=2*n*Pi-Pi/2-u=1/z*cos(u);
On verrouille le produit x*z et on determine u en
fonction de x*z par réversion de série.
Equ[1]:=xz=series(sin(u),u);
U:=solve(Equ[1],u);
On pose nu=1/n pour se ramener en 0 et on explicite
l'autre équation.
Equ[2]:=nu=series(1/subs(u=U,xz=x*z,solve(equ[2],n)),z);
On résout en z. On a ainsi un développement de
z par rapport à nu avec des coefficients en
x.
Z:=solve(Equ[2],z);
On s'intéresse à x qui est -ln(z). On
sort le -log(nu/2/Pi) pour aider series. On a un
développement de -ln(z) par rapport à
nu avec des coefficients en x.
-log(nu/2/Pi)-series(log(Z/(nu/2/Pi)),nu);
On résout en x. On a ainsi un développement de
x par rapport à nu avec des coefficients en
log(nu).
solve(x=",x);
On repasse en n=1/nu et on fait un peu de ménage.
X:=map(collect,asympt(subs(nu=1/n,"),n),ln(n));
Maintenant qu'on a x, on le substitue dans le
développement de z qui comportait des coefficients en
x et on repasse en n=1/nu.
Z:=subs(x=X,nu=1/n,Z);
On injecte le résultat dans l'expression de t en
fonction de z et x et on développe. Ensuite
on nettoie.
T:=asympt(subs(x=X,z=Z,t),n);
T:=map(collect,T,ln(n));
On vérifie que tout baigne dans l'huile par le
numérique.
for k from 100 to 110 do
evalf(eval(subs(O=0,n=k,-X+I*T)))-LambertW(k,1.)
od;
-
Champ d'application de la formule
d'Euler-Maclaurin. Page 218, nous
affirmons qu'un certain exemple ne rentre pas dans le champ
d'application de la formule d'Euler-Maclaurin, parce que le terme
général de la somme dépend de l'indice et de
la borne de la somme. Ceci est une stupidité. La formule
d'Euler-Maclaurin permet par exemple d'obtenir des
développements asymptotiques de sommes de Riemann. On trouve
une jolie présentation de ceci dans [SeFl96a] ou dans sa version
française [SeFl96b].
-
Arithmétique d'intervalles. Pages 237-239, nous présentons
l'arithmétique d'intervalles. La thèse d'Olivier
Didrit (juin 1997) fournit une introduction
pédagogique à ce domaine.
-
Intégration par parties et
développement asymptotique d'une primitive.
Pages 239-243, nous appliquons
l'intégration par parties à la recherche du
développement asymptotique d'une primitive. Nous sommes
amenés à décortiquer une expression de la forme
c_1 gamma_1 x^{gamma_1-1} + ... + c_r gamma_r x^{gamma_r-1} +
alpha/x + beta/(x\ln{x}). Il serait naturel d'employer series pour
cela, mais les exposants gamma_i peuvent être rationnels ce qui
empêche un fonctionnement adéquat de series. Nous
utilisons donc limit et select; pourtant nous évitons
d'habitude l'emploi de limit. Ici la classe d'expressions
utilisée est suffisamment bien délimitée pour que
cet emploi ne pose pas de problème.
-
Fonction d'Eisenstein.
Page 250, nous ne terminons pas l'étude
de la fonction d'Eisenstein en laissant dans l'ombre le cas 0 < x <
e^(-e). On trouve une étude propre et complète de la
question dans Analyse, Exercices avec solutions, 1, E. Ramis,
C. Deschamps, J. Odoux, Masson, 1991, à la page 110, exercice
3.4.12.
-
Écriture d'une somme
hypergéométrique. Page
261, nous écrivons une procédure
convert/Sum qui consiste essentiellement en une
substitution. Il serait possible de procéder
récursivement. Dans un souci d'efficacité, nous
procédons autrement. En effet les expressions Maple ne sont pas des
arbres mais des dags (directed acyclic graphs). Les
sous-expressions communes sont partagées. Ceci est bien
expliqué dans [GoSaZi95]. Si l'on
procède récursivement, une sous-expression va être
traitée autant de fois qu'elle apparaît dans
l'expression. Avec la procédure que nous écrivons, les
sous-expressions qui doivent subir la substitution sont d'abord
extraites de l'expression, à l'aide de indets, qui
emploie réellement la structure de dag. Pour chacune de ces
sous-expressions la substitution à opérer est
calculée et gardée dans une table. Enfin toutes les
substitutions sont effectuées. Le malheur est que subs
n'emploie pas la structure de dag mais seulement la structure
d'arbre. Cela est bien choquant et cette locution ne traduit que
faiblement notre pensée. En tout cas notre coeur est pur et
notre programmation sans tache.
-
Fonstions holonomes.
Page 262, nous employons le package
gfun pour déterminer une équation différentielle
vérifiée par une fonction donnée. Le calcul n'est pas
réellement satisfaisant puisque la procédure
employée devine l'équation différentielle, sans
garantie. Le package gfun comporte actuellement une
procédure holexprtodiffeq nettement plus
satisfaisante, car garantie, mais qui n'est pas intégrée
à la version V.4. Par contre la version de gfun qui est dans la
bibliothèque partagée des utilisateurs (share
library) pour la version V.5 comporte
holexprtodiffeq. Avec la version V.5, on peut donc employer
le package à l'aide des instructions suivantes (la
deuxième n'est pas indispensable elle fournit une
description du package).
> with(share);
> ?share,gfun
> with(gfun);
Comme indiqué dans le livre, page 432,
on peut trouver la dernière version de gfun dans la
bibliothèque du
Projet Algorithmes.
Une fois que gfun est disponible, on obtient de l'aide sur le
package et sur la procédure holexprtodiffeq
par les instructions suivantes.
> ?gfun[gfun]
> ?gfun[holexprtodiffeq]
À l'aide la version actuelle de gfun et en Maple V.4, le
calcul des pages 262-263 prend la tournure
suivante.
> f:=arcsin(sqrt(x))/sqrt(x-x^2):
> holexprtodiffeq(f,y(x));
/ 2 \
/d \ 2 |d |
2 y(x) + (6 x - 3) |-- y(x)| + (-2 x + 2 x ) |--- y(x)|
\dx / | 2 |
\dx /
La procédure a fourni une équation différentielle
linéaire d'ordre 2 satisfaite par la fonction. On peut noter
que cette équation n'est pas minimale, puisque dans le livre on
a obtenu une équation linéaire d'ordre 1. D'ailleurs
dsolve nous montre bien que l'espace des solutions est plus
gros que la droite engendrée par la fonction.
> dsolve(",y(x));
1/2
_C1 _C2 arcsin(x )
y(x) = ---------------- + ----------------
1/2 1/2 1/2 1/2
x (-1 + x) x (-1 + x)
Cependant l'équation fournie par holexprtodiffeq est
homogène alors que celle fournie par seriestodiffeq ne
l'est pas. Ceci explique la différence de comportement des deux
procédures.
Dans la
bibliothèque du Projet Algorithmes
on trouvera une présentation du package gfun et de la
théorie sous-jacente. Disons simplement que les fonctions
holonomes sont les solutions d'équations
différentielles linéaires à coefficients
polynomiaux. On montre que cet ensemble de fonctions possède
des propriétés de stabilité par la somme et le
produit, entre autres. De plus les fonctions algébriques sont
holonomes (problème de la page 320) et
la composition par une fonction algébrique préserve
l'holonomie. Le calcul précédent est
fondé sur la recherche d'équations
différentielles linéaires satisfaites par arcsin(x),
puis arcsin(sqrt(x)) d'une part, par 1/sqrt(x-x^2) d'autre part, par
leur produit enfin.
-
Phénomène de Gibbs.
L'encart de la figure 6.2, page 277, a
été obtenu par le code suivant, qui complète le
code déjà fourni.
window:=[0..0.1,0.8..1.2]:
picture[0]:=plot(signal,t=0..0.1,color=black,
thickness=3,view=window):
for k to 5 do
N:=100*k;
signalapprox:=a0/2+add(b*sin(n*t),n=1..N);
picture[k]:=plot(signalapprox,t=0..0.1,
color=black,view=window):
od:
plots[display]({seq(picture[k],k=0..5)});
-
Sommation. Page 288,
nous notons que la réponse suivante, fournie par
Maple V.4, n'est pas satisfaisante.
> sum(z^(2*n+1)/(2*n+1),n=0..infinity);
2 1/2
1 + (z )
z ln(-----------)
2 1/2
1 - (z )
1/2 -----------------
2 1/2
(z )
Par contre dans la version V.5, la réponse est bien ce que l'on
attend.
> sum(z^(2*n+1)/(2*n+1),n=0..infinity);
1 + z
1/2 ln(-----)
1 - z
-
Orbites et champ de vecteurs.
On peut se demander pourquoi nous employons l'option
arrows=NONE à la fin du paquet d'instructions qui
figure dans la page 306. En effet sans cette
option, on obtient exactement le même dessin. L'option sur les
flèches fonctionne de la manière suivante. Les solutions
du système sont numériquement calculées ;
puis un rectangle minimal contenant les orbites est
déterminé ; enfin les flèches sont
tracées suivant un maillage régulier du rectangle. Ici
nous avons ajouté l'option view et nous avons choisi
un rectangle beaucoup plus petit que le précédent. Du
coup ce petit rectangle ne contient pas de point du maillage pour les
flèches. Autrement dit, l'option view est mal prise en
compte dans la procédure DEtools/DEplot.
-
Connexion de développements locaux pour des
solutions d'équations différentielles. Dans
l'exercice 7.5, page 312, nous faisons
établir des développements asymptotiques localement
convergents en 0 et à l'infini pour les solutions d'une
équation différentielle linéaire du second
ordre. Cependant rien ne permet de raccorder ces
développements. Sur ce sujet difficile, on pourra consulter le
livre de Nico Temme [Temme96].
-
Méthode d'exclusion.
Page 353, nous proposons l'étude de la
méthode
d'exclusion. Elle est présentée dans un cadre plus
général par Jean-Claude Yakoubsohn dans les
résumés du
Séminaire Algorithmes, section Symbolic
Computation, année 1992-1993, The exclusion
algorithm.
-
Complexité de calcul d'un produit
matriciel. À la page 381, nous
déclarons que le calcul d'un inverse est du même
ordre de complexité que le calcul d'un produit par la
méthode usuelle. Cette toute dernière
précision tient au fait qu'il existe des méthodes plus
performantes de calcul d'un produit matriciel, comme la méthode
de Strassen [Knuth81, page 481], [Sedgewick90], qui a une complexité
en n^log_2(7).
-
Méthode de Kummer. Les estimations de
la deuxième question, page 406 ont
été obtenues de la manière suivante. On a
commencé par modifier la procédure kummer en
changeant à la fin u+FF en [u,FF]. Ensuite on
a utilisé le code que voici.
F:=1/(n^3+1):
epsilon:=10^(-20):
Digits:=25:
pmax:=8:
Kummer[0]:=[0,F]:
for p to pmax do
Kummer[p]:=kummer(F,p,n)
od:
for p from 0 to pmax do
term:=Kummer[p][2];
Numer:=numer(diff(term,n));
d:=degree(Numer,n);
if d>0 then
bound:=add(abs(coeff(Numer,n,d-i)),i=1..d);
bound:=ceil(bound/abs(lcoeff(Numer,n)))
else
bound:=1
fi;
if term<>0 then
asymptK:=convert(asympt(term,n),polynom);
for k from Order+1 while asymptK=0 do
asymptK:=convert(asympt(term,n,k),polynom);
od;
d:=degree(asymptK,n);
approximatebound:=
ceil(evalf((-(d+1)*epsilon/abs(coeff(asymptK,n,d)))^(1/(d+1))));
upperbound:=int(term,n=N..infinity);
for b from max(bound,approximatebound) while
evalf(abs(evalf(subs(N=b,upperbound)))-epsilon)>0
do od;
fi;
Bound[p]:=b;
Estimation[p][partialsum1]:=evalf(sum(Kummer[p][1],n=1..infinity));
Estimation[p][partialsum2]:=evalf(sum(term,n=1..b));
Estimation[p][reste]:=evalf(sum(term,n=b+1..infinity));
od:
seq(Bound[p],p=1..pmax);
seq(Estimation[p][partialsum1],p=1..pmax);
seq(Estimation[p][partialsum2],p=1..pmax);
seq(Estimation[p][reste],p=1..pmax);
seq(Estimation[p][partialsum1]+Estimation[p][partialsum2]+Estimation[p][reste],p=1..pmax);
Bien sûr, le calcul est incohérent puisque nous employons
sum dans les dernières instructions, ce qui signifie
que l'on dispose de formules sommatoires et que le calcul
numérique est inutile. Cependant il ne s'agissait que de tester
la méthode. On notera aussi qu'au delà de p=5
la borne remonte. Comme nous l'avons indiqué, cela tient
à une approximation trop grossière des racines des
polynômes.
Maple