Maple
Son bon usage en mathématiques

Philippe Dumas,
Xavier Gourdon

Springer-Verlag, 1997

Commentaires

mise à jour du 11 avril 2017

• 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.