Calcul formel : Mode d'emploi - Exemples en Maple
Claude Gomez, Bruno Salvy, Paul Zimmermann
Masson, 1995
Chapitre VIII, section 1.4, exercice 1, page 202.
Philippe.Dumas@inria.fr
http://algo.inria.fr/dumas/Maple/
|
|
Nous employons la seconde des deux formules. De plus nous considérons que le calcul du premier terme de l'expression est effectué par evalf. C'est ce que semble signifier l'énoncé.
Le terme général de la série a l'expression suivante.
> u:=1/(2*n+1)^2/binomial(2*n,n);
Nous calculons le quotient de deux termes consécutifs de façon à minimiser le nombre d'opérations.
> r:=convert(u/subs(n=n-1,u),GAMMA);
> simplify(r,GAMMA);
Nous voyons que le passage du terme de rang n-1 au terme de rang n demande sept multiplications ou divisions. Le calcul des n premiers termes et de la somme partielle d'indice n a donc un coût qui est un O(n). Si nous appliquions brutalement la formule, le coût serait un O(n^2).
L'énoncé demande de passer en paramètre le nombre de chiffres désiré. Il est plus dans l'esprit du logiciel d'employer la variable d'environnement Digits pour contrôler la précision du calcul. D'autre part, nous n'effectuons pas de calcul si le nombre de chiffres demandé est faible. Dans ce cas le résultat est tout simplement stocké en mémoire. C'est ce qui est fait pour les constantes connues du système. On pourra s'en convaincre par l'expérience suivante.
> interface(verboseproc=2):
readlib(`evalf/constant/Pi`);
readlib(`evalf/constant/bigPi`);
Si le nombre de chiffres demandé dépasse 50, nous augmentons la précision des calculs pour éviter les erreurs d'arrondis, puis nous sommons les termes de la série, jusqu'à ce qu'ils deviennent négligeables.
> catalan:=proc()
local t,s,n,oldDigits;
if Digits<=50 then
evalf(.915965594177219015054603514932384110774149374281672)
else
oldDigits:=Digits;
Digits:=floor(1.1*Digits);
t:=1;
s:=1;
for n while t>10^(-oldDigits-5) do
t:=evalf(t/2*n*(2*n-1)/((2*n+1)^2));
s:=s+t
od;
evalf(evalf(Pi/8*ln(2+3^(1/2))+3/8*s),oldDigits)
fi
end:
Comparons le résultat fourni avec celui que donne la procédure du système.
> Digits:=200:
> catalan();
> evalf(Catalan);
> Digits:=10:
Tout semble correct.
Le test d'arrêt de la boucle n'a pas été justifié. Il ne suffit pas que le terme d'indice N d'une série soit plus petit que 10^(-D) pour que le reste d'indice N vérifie la même inégalité. Par exemple le terme 1/N^2 est plus petit que 10^(-4) pour N > 100, mais le reste d'indice 101 vaut environ 0.009. Ici les termes décroissent tellement vite vers 0, que le reste d'indice N est plus petit que le terme d'indice N. Pour nous en convaincre nous regardons le comportement du terme général de la série.
> asympt(u,n,3);
> combine(%,power,symbolic);
Nous avons obtenu un équivalent du terme général de la série. Comme il s'agit de séries convergentes, les restes des deux séries sont donc équivalents. Pour cette nouvelle série, la suite des termes décroît à partir d'un certain rang. Le reste d'indice N est donc plus petit que l'intégrale de N à l'infini de l'équivalent obtenu (et plus grand que l'intégrale de N+1 à l'infini).
> int(sqrt(Pi)/4/4^n/n^(3/2),n=N..infinity);
> asympt(%,N,2);
Nous voyons que cette intégrale vaut environ 1/2/ln(2) fois le terme d'indice N, ce qui justifie la méthode employée.
> evalf(1/2/ln(2));