Q1. 30 décimales sans accélération 

La somme étant alternée, l'erreur est majorée par le premier terme négligé.  

>	f:=sum((-1)^(i+1)/i,i=1..n);

 

Pour simplifier, on peut prendre la limite de la somme restante et on se demande quand cette valeur est inférieure à  

>	limit(f,n=infinity)/n^2-10^(-30);

 

>	solve(%,n);

 

>	evalf(%[2]);

 

Il faut donc de l'ordre de termes. C'est trop ! 

Q2. Les premiers termes de  

>	S:=proc(N) local i, k; add(add((-1)^i/i,i=1..k)*(-1)^k/k^2,k=1..N) end;

 

>	Digits:=30:

 

>	powmax:=10:

 

>	v:=[seq(evalf(S(2^m)),m=1..powmax)];

 

>	v[-1]-v[-2];

 

Apparemment, ce calcul donne environ 6 décimales. 

Q3. Accélération de la convergence 

>	estimrho:=proc(L) (L[-1]-L[-2])/(L[-2]-L[-3]) end;

 

>	estimrho(v);

 

>	euleracc:=proc(L,rho) local i; [seq((L[i+1]-rho*L[i])/(1-rho),i=1..nops(L)-1)] end;

 

>	euleracc(v,1/4);

 

>	estimrho(%);

 

>	T[0]:=v: for i while T[i-1]<>[] do T[i]:=euleracc(T[i-1],1/2^(i+1)) od;

 

Le nombre de décimales obtenues est estimé par   

>	T[9][-1]-T[8][-1];

 

qui suggère 19. 

Commentaire sur le choix de  

L'estimation empirique ci-dessus peut-être remplacée par une étude du comportement asymptotique du reste, qui est faite en bas de cette session. 

Q4. Les 100 premières valeurs 

>	L:=[seq(add(add((-1)^(i)/i,i=1..k)*(-1)^k/k^2,k=1..N),N=0..100)]:

 

Q5. Une récurrence 

>	rec:=gfun:-listtorec(L,u(n),['ogf']);

 

Q6. Une procédure 

>	p:=gfun:-rectoproc(rec[1],u(n),evalfun=evalf);

 

Q7. Vérification 

>	evalf(p(200)-S(200));

 

>	evalf(p(200)-S(200),40);

 

Q8. Calcul des 15 premiers termes 

>	powmax:=15:

 

>	v:=[seq(p(2^m),m=1..powmax)]:v;

 

Q9. Gain de temps 

>	for i from 10 to 16 do i,time(p(2^i)) od;

 

Le temps de calcul est à peu près linéaire en l'indice. Il faut donc a l'apporche directe un nombre de secondes de l'ordre de  

>	10^15/2^%[1]*%[2];

 

ou en années : 

>	%/365/24/3600;

 

Il est possible de gagner du temps en relachant la précision et en utilisant des flottants machines avec la commande evalhf: 

>	for i from 10 to 20 do i,time(evalhf(p(2^i))) od;

 

ce qui réduit le nombre d'années à attendre à : 

>	10^15/2^%[1]*%[2]/365/24/3600;

 

Le gain apporté par la récurrence est net : on passe d'une complexité quadratique à une complexité linéaire 

>	for i from 10 to 13 do i,time(evalhf(S(2^i)))/time(evalhf(p(2^i))) od;

 

Q10. Accélération de la convergence 

>	T[0]:=v: for i while T[i-1]<>[] do T[i]:=euleracc(T[i-1],1/2^(i+1)) od;

 

>	T[14][-1]-T[13][-1];

 

Ce résultat est trop proche de la précision de calcul. Il vaut donc mieux tout recommencer avec une précision supérieure : 

>	Digits:=50:

 

>	v:=[seq(p(2^m),m=1..powmax)]:

 

>	T[0]:=v: for i while T[i-1]<>[] do T[i]:=euleracc(T[i-1],1/2^(i+1)) od:

 

>	T[14][-1]-T[13][-1];

 

Q11. Identification 

>	res:=evalf(T[14][-1],35);

 

>

?identify

 

>	i:='i':sum((-1)^(k+1)/k^2*sum((-1)^(i+1)/i,i=1..infinity),k=1..infinity);

 

On prend les éléments de BasisSumConst de la page d'aide d'identify, auxquels on ajoute ce nombre : 

>	basis:=[%,1,sqrt(2),sqrt(3),Pi,ln(2),ln(3),Zeta(3),Zeta(5)];

 

>	exact:=identify(res,BasisSumConst=basis);

 

Vérification avec un peu plus de décimales 

>	evalf(exact,50)-T[14][-1];

 

Le comportement asymptotique du reste 

Pour justifier l'utilisation de rho=2^(-i-1) lors de l'accélération, on peut calculer un développement asymptotique du reste. On part du sommant : 

>

i:='i':

 

>	Order:=20:

 

>	asympt((-1)^(k+1)/k^2*sum((-1)^(i+1)/i,i=1..k),k);

 

>	map(simplify,%) assuming k::posint;

 

Le reste doit donc se comporter comme  

>	reste:=add((a[i]*(-1)^n+b[i])/n^(2+i),i=0..Order)+O(1/n^(Order+1));

 

pour des coefficients à déterminer. On construit d'abord un système (linéaire, triangulaire) d'équations pour ces coefficients : 

>	asympt(reste-subs(n=n-1,reste)-subs(k=n,%%),n):

 

>	map(simplify,%) assuming n::posint:

 

>	convert(%,polynom): # supprime le O()

 

>	subs((-1)^n=T, n=1/N,%):

 

>	sys:={coeffs(expand(%),[T,N])};

 

Il n'y a plus qu'à résoudre. 

>	solve(sys,indets(sys,name));

 

>	reste:=asympt(subs(%,reste),n,Order-1);

 

Dans les calculs, ce reste est toujours utilisé pour n pair : 

>	map(simplify,reste) assuming n::even;

 

On note le caractère divergent de ce développement : les coefficients ne tendent pas vers 0. 

>	evalf(%,10);

 

Ceci participe à l'explication du ralentissement de la convergence sur les dernières lignes: 

>	seq(evalf(abs(exact-T[i][-1])),i=1..14);

 

>