Chapitre 13  Algorithmes sur les suites et les séries

L’objectif ici est de traduire les algorithmes en l’écriture de programmes. On écrit ici des programmes permettant d’avoir les termes d’une suite ou d’une série et de trouver des valeurs approchées de leur limites.
Mais, pour édutier les suites et les séries, on peut aussi utiliser le tableur ce qui est souvent plus facile que d’écrire un programme.

13.1  Les suites

Soit u_n une suite de réels définie soit par u_n=f(n), soit par une relation de récurrence u_n=f(u_n−m,..,u_n−1) et la donnée de ses premiers termes. On veut ici, calculer les valeurs de u_n. Pour les fonctions qui suivent, il suffira de rajouter la fonction evalf dans le return pour avoir une valeur approchée de u_n : par exemple return evalf(uk).

13.1.1  Les suites u_n=f(n)

Pour avoir le n-ième terme u_n il suffit :
de définir la fonction f et de taper f(n).
On peut aussi mettre f comme paramètre et taper :

u(f,n):=f(n)

Ainsi u(sq,3) vaut 9 et u(sqrt,3) vaut sqrt(3).
On remarquera qu’il est souvent préférable de simplifier l’écriture de u(f,n) avec la commande normal : mettre plutôt normal dans la définition de f.
Par exemple on définit : f(x):=normal(x/sqrt(3)+sqrt(3))).
On tape :
u(f,3)
On obtient :
2*sqrt(3)
On peut aussi considérer qu’il n’y a qu’un paramètre l qui est la séquence f,n et définir u par :

u(l):=l[0](l[1])

Pour avoir la suite des termes u_k, pour k allant de k0 à n, on écrit :

utermes(f,k0,n):={
local k,lres;
lres:=NULL;
for (k:=k0;k<=n;k++){
     lres:=lres,f(k);
}
return lres;
}

On a choisit de mettre tous les termes cherchés dans une séquence.
On a : lres:=NULL; initialise la séquence à vide.
Par exemple, avec la fonction :
f(x):=normal(x/sqrt(3)+sqrt(3))).
On tape :
utermes(f,0,5)
On obtient :
sqrt(3),4*sqrt(3)/3,5*sqrt(3)/3,2*sqrt(3),7*sqrt(3)/3,8*sqrt(3)/3

13.1.2  La représentation des suites u_n=f(n)

On va représenter une suite u_n par des segments verticaux : le terme u_p sera représenté par le segment joignant le point (p,0) au point (p,f(p)).
Pour faire cette représentation, on définit la fonction f et on valide le programme suivant qui permet de représenter u_j1=f(j1)..u_j2=f(j2).

plotsuite(f,j1,j2):={
local j,P,L;
L:=NULL;
for (j:=j1;j<=j2;j++) {
P:=point(j+i*u(j),couleur=point_width_4+noir);
L:=L,segment(j,P,couleur=ligne_tiret+rouge),P;
}
return L;
};

Exemple
u(n):=1+(-1)^n/n
puis,
plotsuite(u,0,10)
puis,
plotsuite(u,20,30)

13.1.3  La représentation des suites récurrentes u₀=a,u_n=f(u_n−1)

Pour avoir les termes u₀=a,...u_p, on définit la fonction f puis on tape :

plotsuiterec1(f,a,p):={
local j,P,L;
L:=NULL;
a:=evalf(a);
for (j:=0;j<=p;j++) {
P:=point(j+i*a,couleur=point_width_4+noir);
L:=L,segment(j,P,couleur=ligne_tiret+rouge),P;
a:=f(a);
}
return L;
};

Exemple
f(x):=x^2-2
puis,
plotsuiterec1(f,0,10)
ou,
plotsuiterec1(sq-2,0,10)

13.1.4  La représentation des suites récurrentes [u₀,u₁,...u_s−1]=la, u_n=f(u_n−s...,u_n−1) si n>=s

Pour avoir les termes u₀=la[0],u₁=la[1],...u_p, on définit la fonction f puis on tape :

plotsuiterec(f,la,p):={
local j,P,L,s,a;
L:=NULL;
s:=size(la);
la:=evalf(la);
for (j:=0;j<s;j++) {
P:=point(j+i*la[j],couleur=point_width_4+noir);
L:=L,segment(j,P,couleur=ligne_tiret+rouge),P;
}
for (j:=s;j<=p;j++) {
a:=f(op(la));
P:=point(j+i*a,couleur=point_width_4+noir);
L:=L,segment(j,P,couleur=ligne_tiret+rouge),P;
la:=append(tail(la),a);
}
return L;
};

Exemple
f(x,y):=x+y
puis,
plotsuiterec(f,[0,1],6)

13.1.5  L’escargot des suites récurrentes u(0)=a, u(n)=f(u(n−1) si n>0

Cette session se trouve dans plottoile.xws.
On rappelle que la commande plotseq(f(x),a,p) permet de visualiser les p premiers termes de la suite récurrente u₀=a, u_n=f(u_n−1) si n>0 en visualisant "l’escargot". On se propose de réécrire cette commande de façon a bien mettre en evidence la construction des différents termes de la suite.
À la différence de plotseq la fonction plottoile a comme premier argument la fonction f et non l’expression f(x). On représente le premier terme u₀ par une croix noire sur l’axe des x d’abscisse le deuxième argument. Pour avoir u₁, on trace le segment vertical allant de la croix au graphe de f(x), puis le segment horizontal allant, du point du graphe de f(x) au graphe de la première bissectrice (pour reporter la valeur de f(u₀) sur l’axe des x), puis, un segment vertical en pointillés allant, du graphe de la première bissectrice à l’axe des x pour tracer une croix rouge.
Pour avoir les termes suivants, on trace ensuite un segment vertical allant du point du graphe de la première bissectrice, au graphe de f(x), puis le segment horizontal allant, du point du graphe de f(x) au graphe de la première bissectrice, puis, un segment vertical jusqu’à l’axe des x pour tracer une croix rouge etc... On s’arrête lorsque l’on a dessiné p croix rouges, p etant le troisième argument.
On tape :

plottoile(f,u,n):={
local j,v,L,P;
u:=evalf(u);
P:=point(u,couleur=point_width_4+noir);
if (n<=0 ) {return P;}
v:=f(u);
L:=segment(u,u+i*v,couleur=rouge),P;
L:=L,segment(u+i*v,v+i*v,couleur=rouge);
u:=v;
v:=f(u);
P:=point(u,couleur=point_width_4+rouge);
L:=L,segment(u,u+i*u,couleur=ligne_tiret+rouge),P;
for (j:=2;j<=n;j++) {
L:=L,segment(u+i*u,u+i*v,couleur=rouge);
L:=L,segment(u+i*v,v+i*v,couleur=rouge);
u:=v;
v:=f(u);
P:=point(u,couleur=point_width_4+rouge);
L:=L,segment(u,u+i*u,couleur=ligne_tiret+rouge),P;}
return plotfunc(f(x),x,couleur=vert),plotfunc(x,x,couleur=bleu),L;
};

Par exemple pour voir les premiers termes de :
u₀=2, si n ≥ 1, u_n=cos(u_n−1), on tape : plottoile(cos,2,5)
On peut aussi noter les indices des termes de la suite pour la croix représentant u_j en rajoutant legende(u,j,quadrant4).

plottoilegende(f,u,n):={
local j,v,L,P;
u:=evalf(u);
P:=point(u,couleur=point_width_4+noir),legende(u,0,quadrant4);
if (n<=0 ) {return P;}
v:=f(u);
L:=segment(u,u+i*v,couleur=rouge),P;
L:=L,segment(u+i*v,v+i*v,couleur=rouge);
u:=v;
v:=f(u);
P:=point(u,couleur=point_width_4+rouge),legende(u,1,quadrant4);
L:=L,segment(u,u+i*u,couleur=ligne_tiret+rouge),P;
for (j:=2;j<=n;j++) {
L:=L,segment(u+i*u,u+i*v,couleur=rouge);
L:=L,segment(u+i*v,v+i*v,couleur=rouge);
u:=v;
v:=f(u);
P:=point(u,couleur=point_width_4+rouge),legende(u,j,quadrant4);
L:=L,segment(u,u+i*u,couleur=ligne_tiret+rouge),P;}
return plotfunc(f(x),x,couleur=vert),plotfunc(x,x,couleur=bleu),L;
};

Par exemple pour voir les premiers termes avec leur indice de :
u₀=2, si n ≥ 1, u_n=cos(u_n−1), on tape : plottoilegende(cos,2,5)
On peut aussi faire une animation qui montrera la progression de la construction. Pour cela on modifie la fonction plottoile en toile pour avoir dans LT, la progression du tracé. On remarquera que l’on met entre crochet les objets graphiques qui seront affichés simultanément lors de l’animation.
On tape :

toile(f,u,n):={
local j,v,L,P,LT;
u:=evalf(u);
P:=point(u,couleur=point_width_4+noir);
v:=f(u);
LT:=P;
L:=segment(u,u+i*v,couleur=rouge);
L:=L,segment(u+i*v,v+i*v,couleur=rouge);
u:=v;
v:=f(u);
P:=point(u,couleur=point_width_4+rouge);
L:=L,segment(u,u+i*u,couleur=ligne_tiret+rouge),P;
LT:=LT,[LT,L];
for (j:=2;j<=n;j++) {
L:=L,segment(u+i*u,u+i*v,couleur=rouge);
L:=L,segment(u+i*v,v+i*v,couleur=rouge);
u:=v;
v:=f(u);
P:=point(u,couleur=point_width_4+rouge);
L:=L,segment(u,u+i*u,couleur=ligne_tiret+rouge),P;
LT:=LT,[LT,L];
}
return LT;
};

Puis on anime la liste LT renvoyée par toile.

animtoile(f,u,n):={
local LT;
LT:=toile(f,u,n);
return plotfunc(f(x),x,couleur=vert),
       plotfunc(x,x,couleur=bleu),
       animation(LT);
};

Par exemple pour voir en animation les premiers termes de :
u₀=2, si n ≥ 1, u_n=cos(u_n−1), on tape : animtoile(cos,2,5)
On peut régler la vitesse d’animation avec Menu ->Animation (situé dans le pavé de boutons à droite de la fenêtre graphique).
On peut arrêter l’animation avec le bouton ▸| (à gauche de Menu) : il suffit alors, de cliquer dans la fenêtre graphique, pour que l’animation se déroule au pas à pas.

13.1.6  Les suites récurrentes définies par une fonction de plusieurs variables

Un exemple : la suite de Fibonnacci

Commençons par un exemple : la suite de Fibonnacci définie par :
u₀=a
u₁=b
u_n=u_n−1+u_n−2 pour n ≥ 2
On écrit pour avoir u_n:

fibon(a,b,n):={
local k,uk;
for (k:=2;k<=n;k++) {
    uk:=a+b;
    a:=b;
    b:=uk;
}
return uk;
}

On écrit pour avoir u₀,u₁...u_n:

fibona(a,b,n):={
local k,uk,res;
res:=a,b;
for (k:=2;k<=n;k++) {
    uk:=a+b;
    a:=b;
    b:=uk
    res:=res,uk;
}
return res;
}

On écrit pour avoir u_c,u_c+1...u_n pour c ≥ 0 :

fibonac(a,b,c,n):={
local k,uk,res;
for (k:=2;k<c;k++) {
    uk:=a+b;
    a:=b;
    b:=uk
};
if c>1 res:=NULL else 
    if c==0 {res:=a,b;c:=2;} else 
        if c==1 {res:=b;c:=2};
for (k:=c;k<=n;k++) {
    uk:=a+b;
    a:=b;
    b:=uk
    res:=res,uk;
}
return res;
}

Remarque On peut bien sûr écrire un programme récursif qui donne la valeur de u_n. Mais cela n’est pas efficace car on calcule plusieurs fois le même terme. Car par exemple pour calculer u₅ on doit calculer u₃ et u₄ et pour calculer u₄, il faudra à nouveau calculer u₃ et u₂, donc u₃ sera calculé 2 fois et u₂ sera calculé 3 fois. Pour s’en rendre compte on peut imprimer les valeurs des variables pour chacun des appels récursifs.

On a :
u₀=u0
u₁=u1
u_n=u_n−1+u_n−2 pour n ≥ 2
Donc

fibonr(u0,u1,n):={
if (n==0) {print(u0,u1,n);return u0;}
if (n==1) {print(u0,u1,n)return u1;}
print(u0,u1,n)
return fibonr(u0,u1,n-2)+fibonr(u0,u1,n-1);
}:;

On peut aussi compter le nombre de fois que la fonction a èté appelée c’est à dire le nombre de print du programme fibonr précédent:

fibona(u0,u1,n):={
if (n==0) {return [u0,1];}
if (n==1) {return [u1,1];}
print(u0,u1,n)
return fibona(u0,u1,n-2)+fibona(u0,u1,n-1)+[0,1];
}

On tape : fibona(1,1,6)
On obtient : [13,25]
On a :
1 appel avec n=6,
1 appel avec n=5,
2 appels avec n=4,
3 appels avec n=3,
5 appels avec n=2,
8 appels avec n=1,
5 appels avec n=0,
On remarque que le nombre d’appels est une suite de Fibonacci et on voit que pour calculer u_n on doit calculer u₂ fibon(1,1,n-2 fois !

Suites récurrentes definies par une fonction de m variables

On suppose maintenant que la suite est définie par une relation de récurrence definie par une fonction f de m variables : pour définir la suite on se donne les m premiers termes :
u₀,u₁,..,u_m−1 et la relation :
u_n=f(u_n−m,u_n−m+1,..,u_n−1) pour n≥ m.
On veut calculer u_n, et on suppose que les valeurs de u₀,u₁,..,u_m−1 sont dans la liste l0.
On écrit :

urec(f,n,l0):={
local s,k,uk;
s:=size(l0);
l0:=op(l0);
for (k:=s;k<=n;k++) {
    uk:=f(l0);
    l0:=tail(l0),uk;
}
return uk;
}

On utilise op au début, pour transformer la liste l0 en une séquence et tail(l0) pour enlever le premier élément et ainsi l0:=tail(l0),uk est une séquence qui a toujours s éléments.
On peut aussi considérer que le paramètre l contient toutes les variables à savoir l=f,n,u₀,..,u_m−1 . On écrit mais c’est inutilement compliqué (!) :

urecs(l):={
local f,n,s,k,uk;
f:=l[0];
n:=l[1];
l:=tail(tail(l));
s:=size(l);
//f est une fonction de s variables
for (k:=s;k<=n;k++) {
    uk:=f(l);
    l:=tail(l),uk;
}
return uk;
}

Pour avoir tous les termes u_k de la suite pour k allant de 0 à n, On considère que le paramètre l contient toutes les variables à savoir l=f,n,u₀,..,u_m−1.
On écrit :

urec_termes(l):={
local f,n,s,k,uk,lres;
f:=l[0];
n:=l[2];
l:=tail(tail(tail(l)));
s:=size(l);
//f est une fonction de s variables
lres:=l;
for (k:=s;k<=n;k++) {
    uk:=f(l);
    lres:=lres,uk;
    l:=tail(l),uk;
}
return lres;
}

Par exemple on définit :
f(x,y):=normal(x+y)
On tape :
urec_termes(f,5,1,1)
On obtient la suite de Fibonacci :
1,1,2,3,5,8
On tape :
urec_termes(f,5,1,(sqrt(5)+1)/2)
On obtient :
1,(sqrt(5)+1)/2,(sqrt(5)+3)/2,sqrt(5)+2,
(3*sqrt(5)+7)/2,(5*sqrt(5)+11)/2
On tape, pour vérifier que l’on a obtenu la suite géométrique de raison (sqrt(5)+1)/2 :
seq(normal(((sqrt(5)+1)/2)^k),k=0..5)
On obtient :
1,(sqrt(5)+1)/2,(sqrt(5)+3)/2,sqrt(5)+2,
(3*sqrt(5)+7)/2,(5*sqrt(5)+11)/2

Pour avoir tous les termes u_k de la suite pour k allant de k0 à n, On considère que le paramètre l contient toutes les variables à savoir l=f,k0,n,u₀,..,u_m−1.
On écrit :

urec_termekn(l):={
local f,n,s,k,uk,k0,lres;
f:=l[0];
k0:=l[1];
n:=l[2];
l:=tail(tail(tail(l)));
s:=size(l);
//f est une fonction de s variables
for (k:=s;k<k0;k++) {
    uk:=f(l);
    l:=tail(l),uk;
};
if k0>1 res:=NULL else 
  if k0==0 {res:=a,b;k0:=2;} else 
    if k0==1 {res:=b;k0:=2};
for (k:=k0;k<=n;k++) {
    uk:=f(l);
    lres:=lres,uk;
    l:=tail(l),uk;
}
return lres;
}

Par exemple on définit :
f(x,y):=normal(x+y)
On tape :
urec_termekn(f,5,10,1,1)
On obtient la suite de Fibonacci :
8,13,21,34,55,89
On tape :
urec_termes(f,5,9,1,(sqrt(5)+1)/2)
On obtient :
5*sqrt(5)+11)/2,4*sqrt(5)+9,(13*sqrt(5)+29)/2,
(21*sqrt(5)+47)/2,17*sqrt(5)+38

13.2  Les séries

Soit u_n une suite de réels telle que la série ∑_k=0^∞u_k converge vers S. On veut ici, calculer une valeur appochée de cette somme. Si la série converge rapidement, il suffit de calculer ∑_k=0ⁿ u_k pour n assez grand, sinon il faut procéder à une accélération de convergence, en construisant une série de même somme et convergeant plus rapidement.

13.2.1  Les sommes partielles

On écrit :

sum_serie(f,n0,n):={
local s,k;
//un=f(n) ou f est une fonction de 1 variable
s:=0;
for (k:=n0;k<=n;k++) {
    s:=s+evalf(f(k));
}
return s;
}

Il est plus précis de faire le calcul de la somme en commençant par les plus petits termes, on écrit :

serie_sum(f,n0,n):={
local s,k;
//un=f(n) ou f est une fonction de 1 variable
s:=0;
for (k:=n;k>=n0;k--) {
    s:=s+evalf(f(k));
}
return s;
}

On peut avoir aussi besoin de la suite des sommes partielles : par exemple pour les séries alternées deux sommes partielles successives encadrent la somme de la série.
On écrit en utilisant un paramètre supplémentaire alt pour repérer les séries alternées de la forme u_n=altⁿ*f(n) :

sums_serie(f,n0,n,alt):={
local ls,s,k;
//un=(alt)^n*f(n) ou f est une fonction de 1 variable
s:=0;
ls:=[];
if (alt<0){
  if (irem(n0,2)==0) {alt:=-alt;}
  for (k:=n0;k<=n;k++) {
    s:=s+evalf(alt*f(k));
    alt:=-alt;
    ls:=concat(ls,s);
  }
}
else {
  for (k:=n0;k<=n;k++) {
    s:=s+evalf(alt*f(k));
    ls:=concat(ls,s);
  } 
}
return ls;
}

13.2.2  Un exemple simple : une approximation de e

On va calculer la somme :
∑_k=0ⁿ1/k! en commencant par les plus petits termes :

vale0(n):={
local S,k;
S:=0;
for (k:=n;k>=0;k--){
S:=S+1/k!;
}
return S;
}
:;

On tape :
SS:=vale0(22)
iquo(numer(SS)*10^22,denom(SS))
On obtient :
27182818284590452353602
On va calculer la somme :
S=∑_k=0ⁿ1/k! en calculant k! au fur et à mesure dans f

vale1(n):={
local S,f,k;
f:=1;
S:=1;
for (k:=1;k<=n;k++){
f:=f*k;
S:=S+1/f;
}
return S;
}
:;

On tape :
S:=vale1(22)
iquo(numer(S)*10^22,denom(S))
On obtient :
27182818284590452353602
On va calculer la somme :
S=∑_k=0ⁿ1/k! en calculant son numérateur p et son dénominateur f à chaque étape : on a p/(k−1)!+1/k!=(p*k+1)/k! et k!=f*k. Le résultat obtenu en cherchant le quotient de 10ⁿ*p par f donnera les n chiffres significatils de e i.e e=vale(n)*10⁻ⁿ lorsque n≥ 22 car 1./23!<4e−23

vale(n):={
local p,f,k;
f:=1;
p:=1;
for (k:=1;k<=n;k++){
f:=f*k;
p:=p*k+1;
}
return iquo(p*10^n,f);
}
:;

On tape :
vale(22)
On obtient :
27182818284590452353602

13.2.3  Exemple d’accélération de convergence des séries à termes positifs

On suppose que u_n=f(n) et que f(n) admet un développement limité à tous les ordres par rapport à 1/n.
On suppose que u_k∼ a/k^p et on pose :
v_k=u_k−a/(k+1)(k+2)...(k+p)
On a alors, v_k=O(1/k^p+1) et on connait :
∑_k=0^∞a/(k+1)(k+2)...(k+p)
En effet :
a/(k+1)(k+2)...(k+p)=a/p−1(1/(k+1)(k+2)...(k+p−1)−a/(k+2)(k+3)...(k+p)) donc
∑_k=0^∞a/(k+1)(k+2)...(k+p)=a/p−1(1/1· 2·...· (p−1))=a/(p−1)(p−1)! et,
∑_k=k0^∞a/(k+1)(k+2)...(k+p)=a/(p−1)(k0+1)(k0+2)..(k0+p−1)
On a :
∑_k=0^∞u_k=a/(p−1)(p−1)!+∑_k=0^∞v_k
On peut ensuite continuer à appliquer la même méthode à v_k.
Exercice
Utiliser cette méthode pour calculer numériquement : ∑_k=0^∞1/(k+1)².
On va faire "à la main" trois accélérations successives.
On pose :
u_k =1/(k+1)²

• 1-ière accélération :
  u_k =1/(k+1)²∼ 1/k² donc on pose
  v_k=u_k −1/(k+1)(k+2), et donc
  v_k=1/(k+1)²(k+2)
  puisque 1/(k+1)(k+2)=1/(k+1)−1/(k+2), on a :
  ∑_k=0^∞1/(k+1)(k+2)=1 donc
  ∑_k=0^∞u_k=1+∑_k=0^∞v_k
• 2-ième accélération :
  v_k=1/(k+1)²(k+2) ∼ 1/k³ donc on pose
  w_k=v_k −1/(k+1)(k+2)(k+3), et donc
  w_k=2/(k+1)²(k+2)(k+3)
  puisque 2/(k+1)(k+2)(k+3)=1/(k+1)(k+2)−1/(k+2)(k+3), on a :
  ∑_k=0^∞1/(k+1)(k+2)(k+3)=1/2*2!=1/4 donc
  ∑_k=0^∞u_k=1+1/2· 2!+∑_k=0^∞w_k
• 3-ième accélération :
  w_k=2/(k+1)²(k+2)(k+3)∼ 2/k⁴ donc on pose
  t_k=w_k −2/(k+1)(k+2)(k+3)(k+4), et donc
  t_k=6/(k+1)²(k+2)(k+3)(k+4)
  puisque 3/(k+1)(k+2)(k+3)(k+4)=1/(k+1)(k+2)(k+3)−1/(k+2)(k+3)(k+4), on a :
  ∑_k=0^∞2/(k+1)(k+2)(k+3)(k+4)=2/3*3!=1/18 donc
  ∑_k=0^∞u_k=1+1/2* 2!+2/3*3!+∑_k=0^∞t_k

On tape :
u(k):=1/(k+1)^2
On tape :
v(k):=1/((k+1)^2*(k+2))
On tape :
w(k):=2/((k+1)^2*(k+2)*(k+3))
On tape :
t(k):=6/((k+1)^2*(k+2)(k+3)(k+4))
On compare π²/6 et les valeurs obtenues pour n=200, car on sait que :
S=∑_k=0^∞1/(k+1)²=π²/6≃ 1.64493406685
On tape :
serie_sum(u,0,200)
ou
evalf(sum(1/(k+1)^2,k=0..200))
On obtient S à 5*10⁻³ prés (1 décimale exacte) :
1.63997129788
On tape :
1+serie_sum(v,0,200)
ou
evalf(1+sum(1/((k+1)^2*(k+2)),k=0..200))
On obtient S à 1.25*10⁻⁵ prés (4 décimales exactes) :
1.64492179293
On tape :
1+1/4+serie_sum(w,0,200)
ou
evalf(1+1/4+sum(2/((k+1)^2*(k+2)*(k+3)),k=0..200))
On obtient S à 8.3*10⁻⁸ prés (5 décimales exactes) :
1.64493398626
On tape :
1+1/4+1/9+serie_sum(t,0,200)
ou
evalf(1+1/4+1/9+sum(6/((k+1)^2*(k+2)*(k+3)*(k+4)),k=0..200))
On obtient S à 9.2*10⁻¹⁰ prés (8 décimales exactes) :
1.64493406596
Les erreurs
Si on compare la somme ∑_k=n+1^∞1/(k+1)² à une intégrale on a :

Ou encore, on peut aussi remarquer que :

puisque 1/k(k+1)=1/k−1/k+1.
Au bout de la p-ième accéleration on calcule la somme de :
u_k^(p)=p!/(k+1)²(k+2)...(k+p+1) et on a :

Et puisque :

On a :

Donc

On vérifie (π²/6≃ 1.64493406685) :
1.63997129788<π²/6<1.63997129788+1/201=1.64494642226
1.64492179293<π²/6<1.64492179293+1/(2*201²)=1.64493416886
1.64493398626<π²/6<1.64493398626+2/(3*201³)=1.64493406836
1.64493406596<π²/6<1.64493406596+6/(4*201⁴)=1.64493406688
Le programme
On peut écrire un programme qui va demander le nombre d’accélérations pour calculer ∑_k=0^∞1/(k+1)²

serie_sumacc(n,acc):={
local p,l,j,k,ls,sf,sg,gk,fact;
ls:=[];
//calcul sans acceleration
sf:=0.0;
for (k:=n;k>=0;k--) {
   sf:=sf+1/(k+1)^2;
}
ls:=[sf];
sf:=0.0;
fact:=1;
for (p:=1;p<=acc;p++){
  //calcul de 1+1/4+..+1/p^2, le terme a rajouter
  sf:=sf+evalf(1/p^2);
  //calcul de p!
  fact:=fact*(p);
  //calcul de sg, somme(de 0 a n) de la serie acceleree p fois
  sg:=0.0;
  for (k:=0;k<=n;k++) {
    gk:=1/(k+1)^2;
    //calcul du k-ieme terme gk de la serie acceleree p fois (sans p!)
    for (j:=1;j<=p;j++) {
      gk:=evalf(gk/(k+j+1));
    }
    sg:=sg+gk;
  }
  ls:=concat(ls,sf+fact*sg);
}
return(ls);
}

13.3  Méthodes d’accélération de convergence des séries alternées

13.3.1  Un exemple d’accélération de convergence des séries alternées

Un premier exemple

On suppose que u_k=(−1)^kf(k) avec f(k) tend vers zéro quand k tend vers +∞ et f décroissante de ℝ⁺ dans ℝ⁺.
On pose :
g(x)=1/2(f(x)−f(x+1)) donc
v_k=(−1)^kf(k)−f(k+1)/2=u_k+u_k+1/2=(−1)^k g(k)
On a :
∑_k=0ⁿ v_k=1/2(∑_k=0ⁿ u_k+∑_k=0ⁿ u_k+1)
donc,
∑_k=0ⁿ v_k=1/2(∑_k=0ⁿ u_k+∑_k=1ⁿ⁺¹ u_k)
donc,
∑_k=0ⁿ v_k=u₀/2+u_n+1/2+∑_k=0ⁿ u_k
Puisque f(k) tend vers zéro quand k tend vers +∞, g(k)=1/2(f(k)−f(k+1)) tend aussi vers zéro quand k tend vers +∞.
Si la fonction f est convexe (f″(x)>0), la série ∑_k=0^∞v_k vérifie aussi le théorème des séries alternées.
En effet, pour x>0 on a :
g(x)=1/2(f(x)−f(x+1)) ≥ 0 puisque f décroissante sur ℝ⁺
g′(x)=1/2(f′(x)−f′(x+1)) < 0 puisque f″(x)>0, f′ est négative et croissante sur ℝ⁺
donc g est décroissante de ℝ⁺ dans ℝ⁺ et g(k) tend vers zéro quand k tend vers +∞.
Conclusion : La série ∑_k=0^∞v_k est une série alternée de somme S+u₀/2.
Si de plus, f′(x)/f(x) tend vers zéro quand x tend vers l’infini, la série ∑_k=0^∞v_k converge plus rapidement que ∑_k=0^∞u_k, puisque il existe c, x<c<x+1 d’après le th des accroissements finis tel que:
0<g(x)=1/2(f(x)−f(x+1))=−1/2f′(c)
on a donc, puisque f′ est négative et croissante:
0<g(x)<−1/2f′(x)=o(f(x)).

Un exercice

Utiliser cette méthode pour calculer numériquement : ∑_k=0^∞(−1)^k/k+1.
Toutes les dérivées de f(x)=1/(x+1) ont un signe constant sur [0;+∞[ et tendent vers zéro à l’infini, ces dérivées sont donc monotones et on peut donc faire plusieurs accélérations successives.
On va faire "à la main " trois accélérations successives.
On pose :
u_k =(−1)^k/(k+1)

• 1-ière accélération :
  v_k=(−1)^k(1/2(k+1)−1/2(k+2)), et donc
  v_k=(−1)^k(1/2(k+1)(k+2))
  ∑_k=0^∞u_k=1/2+∑_k=0^∞v_k
• 2-ième accélération :
  w_k= (−1)^k(1/4(k+1)(k+2)−1/4(k+2)(k+3)), et donc
  w_k= (−1)^k(1/2(k+1)(k+2)(k+3))
  et comme v₀/2=1/8 on a :
  ∑_k=0^∞u_k=1/2+1/8+∑_k=0^∞w_k
• 3-ième accélération :
  t_k=(−1)^k(1/4(k+1)(k+2)(k+3)−1/4(k+2)(k+3)(k+4)), et donc
  t_k=(−1)^k(3/4(k+1)(k+2)(k+3)(k+4))
  et comme w₀/2=1/24 on a :
  ∑_k=0^∞u_k=1/2+1/8+1/24+∑_k=0^∞t_k

On tape :
u(k):=(-1)^k/(k+1)
On tape :
v(k):=(-1)^k/(2*(k+1)*(k+2))
On tape :
w(k):=(-1)^k/(2*(k+1)*(k+2)*(k+3))
On tape :
t(k):=(-1)^k*3/(4*(k+1)*(k+2)*(k+3)*(k+4))
On compare ln(2) et les valeurs obtenues pour n=200, car on sait que :
S=∑_k=0^∞(−1)^k1/(k+1)=ln(2)≃ 0.69314718056
On tape :
serie_sum(u,0,200)
On obtient S à 5*10⁻³ prés (2 décimales exactes) :
0.69562855486
On tape :
1/2+serie_sum(v,0,200)
On obtient S à 1.23*10⁻⁵ prés (4 décimales exactes) :
0.693153307335
On tape :
1/2+1/8+serie_sum(w,0,200)
On obtient S à 6.1*10⁻⁸ prés (8 décimales exactes) :
0.693147210666
On tape :
1/2+1/8+1/24+serie_sum(t,0,200)
On obtient S à 4.6*10⁻¹⁰ prés (10 décimales exactes) :
0.693147180781
Les erreurs
Le reste d’une série alternée est du signe de son premier terme et la valeur absolue du reste est inférieure à la valeur absolue de son premier terme :
|∑_k=n+1^∞(−1)^k1/(k+1)|<1/(n+2)
Au bout de la p-ième accéleration on calcule la somme de :
u_k^(p)=(−1)^kp!/2^p(k+1)(k+2)...(k+p+1) et on a :

On vérifie (ln(2)≃ 0.69314718055995) :
0.69562855486<ln(2)<0.69562855486+1/202=0.70057904991
0.693153307335<ln(2)<0.693153307335+1/(2*202²)=0.693165561036
0.693147210666<ln(2)<0.693147210666+2/(4*202³)=0.693147271328
0.693147180781<ln(2)<0.693147180781+6/(8*202⁴)=0.693147181231.

Le programme

On peut écrire un programme qui va demander le nombre p d’accélérations.
Si u_k^(p) désigne le k-ième terme de la série accélérée p fois, on a :
∑_k=0^∞(−1)^k/(k+1)=∑_k=0^pu₀^(k−1)/2+∑_k=0^∞u_k^(p)
avec
u_k^(p)=(−1)^k p!/2^p(k+1)...(k+p+1)
On choisit de multiplier seulement à la fin par p!/2^p et de ne calculer que la somme des n premiers termes :

On met cette somme dans la variable sg, pour cela on calcule (−1)^k/(k+1)...(k+p+1) que l’on met dans la variable gk :
au début sg=0 et gk=1/(p+1)! (c’est la valeur pour k=0)
puis, on ajoute gk à la somme sg, ensuite on calcule (−1)¹1!/(p+2)! que l’on met dans gk (c’est la valeur pour k=1) etc...
La variable sf sert au début à calculer ∑_k=0ⁿ (−1)^k/(k+1) puis,
sf sert à calculer la somme à rajouter ∑_k=0^pu₀^(k−1)/2 (qui vaut 1/2+1/8+1/24 pour p=3 accélerations).
Dans le programme, on utilise la variable fact pour calculer (p+1)! et la variable fact2 pour calculer p!/2^p.
On écrit :

seriealt_sumacc(n,acc):={
local l,j,k,ls,sf,sg,gk,fact,fact2,alt,t0,p;
//calcul sans acceleration
sf:=0.0;
alt:=1;
for (k:=n;k>=0;k--) {
sf:=sf+alt/(k+1);
alt:=-alt;
}
if (alt==1) {
ls:=[-sf];} 
else {
ls:=[sf];
}
t0:=0.5;
// sf maintenant est la somme a rajouter
sf:=0.0;
fact:=1;fact2:=1;
for (p:=1;p<=acc;p++){
  sf:=sf+fact2*t0;
  //calcul de p+1! et de p!/2^p
  fact:=fact*(p+1);
  fact2:=fact2*p/2;
//sg, somme(de k=0 a n) de la serie gk acceleree p fois
  sg:=0.0;
//terme d'indice 0 (ds gk) de la serie acceleree p fois 
//(sans p!/2^p=fact2)
  gk:=1/fact;
//on conserve gk/2 dans t0 car il faut rajouter t0 
//au prochain sf 
  t0:=gk/2;
  sg:=sg+gk;
  alt:=-1;
  for (k:=1;k<=n;k++) {
     gk:=1/(k+1);
//terme d'indice k (ds gk) de la serie acceleree p fois 
//(sans p!/2^p=fact2)
     for (j:=1;j<=p;j++) {
       gk:=evalf(gk/(k+j+1));
     }
     sg:=sg+alt*gk;
     alt:=-alt;
  } 
ls:=concat(ls,sf+fact2*sg);
}
return(ls);
}

On met ce programme dans un niveau éditeur de programmes (que l’on ouvre avec Alt+p), puis on le teste et le valide avec OK et on tape dans une ligne de commandes :
seriealt_sumacc(200,3)
On obtient :
[0.69562855486,0.693153307335,0.693147210666,0.693147180781] On tape :
seriealt_sumacc(100,4)
On obtient :
[0.698073169409,0.693171208625,0.693147412699,
0.693147183892,0.693147180623]

13.3.2  La transformation d’Euler pour les series alternées

La transformation d’Euler

On cherche une approximation de :
∑_n=0^∞(−1)ⁿ*u(n)=sum((-1)^n*u(n),n,0,infinity) lorsque u(n) tend vers 0 en décroissant.
On pose : Delta(u))(n)=u(n+1)-u(n) et
delta(u,p,n)=(Delta@@p(u))(n)
On a :
delta(u,2,n)=u(n+2)-2*u(n+1)+u(n)
delta(u,3,n)=u(n+3)-3*u(n+2)+3*u(n+1)-u(n)
delta(u,p,N)=u(n+p)-comb(p,1)*u(n+p-1)+comb(p,2)*u(n+p-2)+
         ....+(-1)^p*u(n)
c’est à dire :
delta(u,p,n)=sum((-1)^(p-j)*comb(p,j)*u(n+j),j,0,p)
La transformation d’Euler consiste à écrire :
sum((-1)^n*u(n),n,N,infinity)
sous la forme :
(-1)^N*sum((-1)^p*delta(u,p,N)/2^(p+1),p,0,infinity)
Pour prouver cette égalité il suffit de développer la dernière expression et de chercher le coefficient de u(N+k) dans la somme :
∑_p=0^∞,(−1)^p*delta(u,p,N)/2^p+1
Le coefficient de u(N+k) est :
s(k)=(-1)^k*sum(comb(k+p,p)/2^(k+p+1),p,0,infinity)
et cette somme vaut (-1)^k quelque soit k entier.
En effet par récurrence :
pour k=0, comb(k+p,p)=1 et
sum(1/2^(p+1),p,0,infinity)=1/2+1/4+...1/2^n+...=1
On a de plus :
- pour p=0, comb(k+p,p)=comb(k+1+p,p)=1
- pour p>0, comb(k+p,p)=comb(k+1+p,p)-comb(k+1+p-1,p-1)
donc
s(k)=(-1)^k*sum(comb(k+1+p,p)/2^(k+p+1),p,0,infinity)-
(-1)^k*sum(comb(k+1+p-1,p-1)/2^(k+1+p-1+1),p,1,infinity)=
-2*s(k+1)-
(-1)^k*sum(comb(k+1+p,p)/2^(k+1+p+1),p,0,infinity)=
-2*s(k+1)+s(k+1)=-s(k+1).
donc si s(k)=(-1)^k alors s(k+1)=(-1)^(k+1).
La transformation d’Euler permet une accélération de convergence car la série :
sum((-1)^p*delta(u,p,N)/2^(p+1),p,0,infinity)
converge plus rapidement.

Le programme

On définit, tout d’abord, la fonction delta :

delta(u,p,n):={
  local val,k,s;
  val:=0;
  s:=1;
  for (k:=p;k>=0;k--) {
  val:=val+comb(p,k)*u(n+k)*s;
  s:=s*-1;
  }
  return val;
};

On écrit la transforpmation d’Euler :
trans_euler(u,N,M) qui approche
sum((-1)^n*u(n),n,0,infinity) et vaut :
sum((-1)^n*u(n),n,0,N-1)+
(-1)^N*sum((-1/2)^p*delta(u,p,N)/2,p,0,M).

trans_euler(u,N,M):={
  local S,T,k,s;
  S:=0;
  s:=1;
  for (k:=0;k<N;k++) {
  S:=S+u(k)*s; 
  s:=s*-1;
  }
  T:=0;
  s:=s*1/2;
for (k:=0;k<=M;k++) {
  T:=T+delta(u,k,N)*s; 
  s:=s*-1/2;
};
  return evalf(normal(S+T));
};

Par exemple pour u(n)=1/(n+1) avec 20 digits, on tape :
u(n):=1/(n+1);
DIGITS:=20;
trans_euler(u,10,20);
On obtient :
0.693147180559945056511
trans_euler(u,9,21);
On obtient :
0.693147180559945594072
On remarque que l’on a 16 decimales exactes car on a :
evalf(ln(2))=0.693147180559945309415

13.3.3  Autre approximation d’une série alternée

La méthode présentée dans cette section est trés largement inspirée par le texte "Sommation de séries alternées" de l’épreuve de modélisation de l’agrégation de mathématiques (session 2006).

Le problème

On veut évaluer la somme S de la série alternée : S=∑_n=0^∞(−1)ⁿa_n
avec (a_n)_{n≥ 0} est une suite de nombres positifs qui tend vers 0 en décroissant.
On suppose que l’on a pour n≥ 0 :
a−n=∫₀¹xⁿdµ
où µ est une mesure positive sur [0,1].
C’est en particulier le cas si a_n=A(n) avec A fonction indefiniment dérivable pour laquelle les kième dérivées A^(k) sont telles que (−1)^k*A^(k)(x) soit positif pour x≥ 0 pour tout k≥ 0.

Le théorème

Théorème :
Soit P_n une suite de polynômes de degré n vérifiant P_n(−1)≠ 0.
À P_n, on associe les coefficients c_n,k pour 0 ≤ k<n définis par :
P_n(−1)−P_n(x)/1+x=∑_k=0ⁿ⁻¹c_n,kx^k
et le coefficient d_n défini par :
d_n=P_n(−1)
Soient S=∑_k=0^∞(−1)^ka_k et S_n=1/d_n∑_k=0ⁿ⁻¹c_n,ka^k
Alors :
|S−S_n|≥ sup_{x ∈ [0,1]} |P_n(x)|/|d_n|S
On a, en effet, avec l’hypothèse faite sur les a_k :
S=∫₀¹ 1/1+xdµ et
S−S_n=∫₀¹P_n(x)/d_n(1+x)dµ

Le choix des polynômes P_n

Pour calculer S il reste à choisir la suite des polynômes P_n.
On peut choisir :

• P_n(x)=(1−x)ⁿ
  on aura une convergence en 2⁻ⁿ car d_n=2ⁿ et sup_{x ∈ [0,1]} |P_n(x)|=1.
  On a :
  P₀(x)=1
  P₁(x)=1−x
  P_n+1(x)=P_n(x)*(1−x) si n≥ 0
  d_n=2ⁿ
  et la formule explicite de P_n :
  P_n(x)=∑_k=0ⁿ(−1)^k C_n^kx^k=∑_k=0ⁿ p_n,k x^k si n≥ 0
  donc les coefficients p_n,k vérifient :
  p_n,0=1
  p_n,k=p_n,k−1*(k−1−n)/k pour 1 ≤ k<n
  On a :
  d_n−P_n(x)=(1+x)∑_k=0ⁿ⁻¹c_n,kx^k=c_n,0∑_k=1ⁿ(c_n,k−1+c_n,k)x^k
  donc
  c_n,0=d_n−p_n,0=d_n−1
  c_n,k=−c_n,k−1−p_n,k pour 1 ≤ k<n
• P_n(x)=x^q(1−x)^2q et n=3*q
  on aura une convergence en 3⁻ⁿ car
  |d_n|=2^2q et
  sup_{x ∈ [0,1]} |P_n(x|)=P_n(1/3)=2^2q*3⁻ⁿ.
  On a :
  P₀(x)=1
  P₃(x)=x(1−x)²
  P_n+3(x)=P_n(x)*x*(1−x)² si n≥ 0
  d_n=(−1)^q*2^2q
  et la formule explicite de P_n :
  P_n(x)=∑_k=0^2q(−1)^k C_2q^kx^k+q=∑_k=qⁿ(−1)^k−q C_2q^k−qx^k=
  ∑_k=0ⁿ p_n,k x^k si n≥ 0
  donc les coefficients p_n,k vérifient :
  p_n,k=0 si k<q
  p_n,q=1
  et comme C_n^p=C_n^p−1*(n−p+1)/p
  p_n,k=(−1)^k−q C_2q^k−q si q ≤ k ≤ n
  p_n,k=−p_n,k−1*(k−1−n)/(k−q) si q < k ≤ n
  On a :
  d_n−P_n(x)=(1+x)∑_k=0ⁿ⁻¹c_n,kx^k=c_n,0+∑_k=1ⁿ(c_n,k−1+c_n,k)x^k
  donc
  c_n,0=d_n−p_n,0=d_n−1
  c_n,k=−c_n,k−1−p_n,k pour 1 ≤ k<n
• Si le ploynôme P_n défini P_n(sin(t)²)=cos(2nt)
  P_n est défini à partir du ploynôme T_n de Chebyshev (T_n(cos(t))=cos(nt)) :
  P_n(sin(t)²)=cos(2nt)T_n(cos(2t))
  on a donc puisque 1−2sin(t)²=cos(2t) :
  T_n(1−2sin(t)²)=cos(2nt) et
  P_n(x)=T_n(1−2x)
  on aura une convergence meilleure que dans les cas précédents car la convergence est en :
  2/((3+√8)ⁿ+(3−√8)ⁿ)≃ 2/((3+√8)ⁿ)≃ 2/(5.8)ⁿ
  car
  d_n=((3+√8)ⁿ+(3−√8)ⁿ)/2 et sup_{x ∈ [0,1]} |P_n(x)|=1.
  Donc :
  P₀(x)=1
  P₁(x)=1−2x
  P_n+2(x)=2(1−2x)P_n+1(x)−P_n(x) si n≥ 0
  d_n=((3+√8)ⁿ+(3−√8)ⁿ)/2
  et la formule explicite de P_n :
  P_n(x)=∑_k=0ⁿ (−1)^k n/n+kC_n+k^2k2^2kx^k=
  ∑_k=0ⁿ p_n,k x^k si n≥ 0
  donc les coefficients p_n,k vérifient :
  p_n,0=1 p_n,k= p_k−1,n(k−1+n)(k−1−n)/((k−1/2)(k)) pour 1 ≤ k<n On a :
  d_n−P_n(x)=(1+x)∑_k=0ⁿ⁻¹c_n,kx^k=c_n,0∑_k=1ⁿ(c_n,k−1+c_n,k)x^k donc
  c_n,0=d_n−p_0,n c_n,k=−c_n,k−1−p_n,k pour 1 ≤ k<n

Les formules de récurrences et le programme pour le polynôme Chebyshev

• Les formules de récurrences
  On va calculer les coefficients c_n,k de proche en proche pour n fixé.
  On pose :
  p:=1;
  d:=((3+sqrt(8))^n+(3-sqrt(8))^n)/2;
  c:=d-p;
  Le premier terme de S_n :
  S:=a(0)*c;
  puis, pour k:= 1 jusque k:= n-1 on calcule p_n,k et c_n,k :
  p:=p*(k+n-1)*(k-n-1)/(k-1/2)/k;
  c:=-p-c;
  On ajoute le kième terme de S_n :
  S:=S+a(k)*c;
  
• Le programme
  

  //n=nombres de termes et a fonction definissant a(n)
  //S_n(P_n) =seriealt(n,a) 
  //S_n(P_n) approche sum((-1)^k*a(k),k,0,+infinity)
  //avec P_n=poly de chebyshev
  seriealt1(n,a):={
  local k,d,c,p,S;
  d:=((3+sqrt(8))^n+(3-sqrt(8))^n)/2;
  p:=1;
  c:=d-p;
  S:=a(0)*c;
  for (k:=1;k<n;k++) {
  p:=p*(k+n-1)*(k-n-1)/(k-1/2)/k;
  c:=-p-c;
  S:=S+a(k)*c;
  }
  return evalf(S/d);
  };
  

Les formules et le programme pour le polynôme P_n(x)=(1−x)ⁿ

• Les formules de récurrences
  On va calculer les coefficients c_n,k de proche en proche pour n fixé.
  On pose :
  p:=1;
  d:=2^n;
  c:=d-p;
  Le premier terme de S_n :
  S:=a(0)*c;
  puis, pour k:= 1 jusque k:= n-1 on calcule p_n,k et c_n,k :
  p:=p*(k-n-1)/k;
  c:=-p-c;
  On ajoute le kième terme de S_n :
  S:=S+a(k)*c;
• Le programme

  //n=nombres de termes et a fonction definissant a(n)
  //S_n(P_n) =seriealt(n,a) 
  //S_n(P_n) approche sum((-1)^k*a(k),k,0,+infinity)
  //avec P_n(x)=poly (1-x)^n
  seriealt2(n,a):={
  local k,d,c,p,S;
  d:=2^n;
  p:=1;
  c:=d-p;
  S:=a(0)*c;
  for (k:=1;k<n;k++) {
  p:=p*(k-n-1)/k;
  c:=-p-c;
  S:=S+a(k)*c;
  }
  return evalf(S/d);
  };
  

Les formules et le programme pour le polynôme P_3q(x)=x^q(1−x)^2q

• Les formules de récurrences
  On va calculer les coefficients c_n,k de proche en proche pour n fixé.
  On pose :
  p:=0; si 0 ≤ k<q
  p:=1; si k=q
  d:=(-1)^q*2^2q;
  c:=d-p;
  Le premier terme de S_n :
  S:=a(0)*c;
  pour k:=1 jusque k:=q-1 on a p_n,k=0 et on calcule c_n,k (c:=-p-c;) et on ajoute le kième terme de S_n :
  S:=S+a(k)*c;
  puis, pour k:=q on a p_n,q=1 et on calcule c_n,q (c:=-p-c;) et on ajoute le qième terme de S_n :
  S:=S+a(q)*c;
  puis, pour k:=q+1 jusque k:= n-1 on calcule p_n,k et c_n,k :
  p:=p*(k-n-1)/k;
  c:=-p-c;
  On ajoute le kième terme de S_n :
  S:=S+a(k)*c;
  
• Le programme
  

  //n=nombres de termes et a fonction definissant a(n)
  //S_n(P_n) =seriealt(n,a) 
  //S_n(P_n) approche sum((-1)^k*a(k),k,0,+infinity)
  //avec P_n=poly de chebyshev
  seriealt3(n,a):={
  local k,d,c,p,q,S;
  q:=ceil(n/3);
  n:=3*q;
  d:=(-1)^q*2^(2*q);
  p:=0;
  c:=d-p;
  S:=a(0)*c;
  for (k:=1;k<q;k++) {
  c:=-p-c;
  S:=S+a(k)*c;
  }
  p:=1;
  c:=-c-p;
  S:=S+a(q)*c;
  for (k:=q+1;k<n;k++) {
  p:=p*(k-n-1)/(k-q);
  c:=-p-c;
  S:=S+a(k)*c;
  }
  return evalf(S/d);
  };
  

Les essais

On choisit n=20.
On tape :
evalf(2/(3+sqrt(8))^20,2^-20,3^-21)=
9.77243031253e-16,9.53674316406e-07,9.55990663597e-11
On a donc pour n=20 une approximation en 10⁻¹⁵ pour Chebyshev, en 10⁻⁶ pour (1−x)²⁰ et en 10⁻¹⁰ pour x⁷(1−x)¹4 :
On choisit dans la suite Digits:=20
Pour calculer une approximation de π/4.
On a :
sum((-1)^n/(2*n+1),n,0,+infinity)=pi/4
On tape :
b(n):=1/(2*n+1)
seriealt1(20,b);evalf(pi/4)
On obtient :
0.785398163397448309926, 0.785398163397448309615
On tape :
seriealt2(20,b);seriealt3(20,b);
On obtient :
0.785397981918786731599, 0.785398163413201025973
Pour calculer une approximation de ln(2).
On a :
sum((-1)^n/(n+1),n,0,+infinity)=ln(2)
On tape :
a(n):=1/(n+1)
seriealt1(30,a);evalf(ln(2))
On obtient :
0.693147180559945311245, 0.693147180559945309415
On tape :
seriealt2(20,a);seriealt3(20,a);
On obtient :
0.693147137051028936275, 0.693147180577738915258

13.3.4  Transformation d’une série en série alternée

On a l’identité formelle :
∑_{n ≥ 1}a_n=∑_{m ≥ 1}(−1)^m−1b_m avec
b_m=∑_{k ≥ 0} 2^k a_2^km.
En effet, si n0 est un entier il existe un entier p0 et un entier impair m0 uniques vérifiant n0=2^p0*m0 . Dans la somme ∑_{m ≥ 1}(−1)^m−1∑_{k ≥ 0} 2^k a_2^km on cherche le coefficient de a_n0, on a soit :
k=0 et m=n0=m0*2^p0, soit
k=1 et m=m0*2^p0−1, soit
................ soit
k=p0 et m=m0.
On remarquera que toutes les valeurs, sauf la dernière, de m sont paires, donc les différentes valeurs de (₁)^m−1 sont (-1) sauf la dernière qui vaut +1.
∑_{m ≥ 1}∑_{k ≥ 0} (−1)^m−12^k a_2^km=
∑_{n0 ≥ 1}a_n0*(∑_k=0^p0−1(−1)*2^k+2^p0)=
∑_{n ≥ 1}a_n puisque 2^p0−∑_k=0^p0−12^k=1

Application au calcul de ∑_n=0^∞1/n^s

Prenons comme exemple la série de terme général a_n=1/n^s avec s>1.
pour s=2
si a(n)=1/n^2
on a :
2^k*a(2^k*m)=1/(2^k*m^2)
b(m)=1/m^2*sum(1/2^k,k,0,+infinity)=2/m^2
pour s quelconque
si a(n)=1/n^s
On a 2^k*a(2^k*m)=1/(2^(k*(s-1))*m^s)
b(m)=1/m^s*sum((1/2^(s-1))^k,k,0,+infinity)
Donc :
b(m)=2^(s-1)/((2^(s-1)-1)*m^s)
pour s=2
b(m):=2/(m^2)
pour s=4
b(m):=8/(7*m^4)
On a :
sum((-1)^(m-1)*b(m),1,+infinity)=
sum((-1)^(m)*b(m+1),0,+infinity)
On choisit encore Digits:=20
pour s=2, ∑_n=1^∞1/n²=π²/6
On tape :
t2(m):=2/(m+1)^2
seriealt1(20,t2),evalf(pi^2/6)
On obtient :
1.64493406684822645248, 1.64493406684822643645
On tape :
seriealt2(20,t2);seriealt3(20,t2);
On obtient :
1.64493374613777534516, 1.64493406688805599300
pour s=4, ∑_n=1^∞1/n⁴=π⁴/90
On tape :
t4(m):=8/(7*(m+1)^4)
seriealt1(20,t4);evalf(pi^4/90)
On obtient :
1.08232323371113822384, 1.08232323371113819149
On tape :
seriealt2(20,t4);seriealt3(20,t4);
On obtient :
1.08232265198912440013, 1.08232323371697925335

Application au calcul de la constante d’Euler

Pour calculer une approximation de la constante d’Euler, (voir aussi 15.2.3 et 13.5.2).
γ=-psi(1).
On a :
-psi(1)=sum((-1)^n*ln(n)/n,n,1,+infinity)/ln(2)+ln(2)/2
et
sum((-1)^n*ln(n)/n,n,1,+infinity)=
-sum((-1)^n*ln(n+1)/(n+1),n,0,+infinity)
c(n):=log(n+1)/(n+1)
-seriealt1(20,c)/ln(2)+ln(2)/2;-evalf(psi(1),0)
On obtient :
0.577215664901532859864, 0.57721566490153
On tape :
-seriealt2(20,c)/ln(2)+ln(2)/2,-seriealt3(20,c)/ln(2)+ln(2)/2
On obtient :
0.577215550220266823551, 0.577215664918305723256 On tape : Digits:=24;
evalf(euler_gamma)
On obtient : 0.5772156649015328606065119

13.4  Polynômes de Bernstein

13.4.1  Définition et théorème

Définition
Le nième polynômes de Bernstein associé à f continuesur [0,1] est :

Théorème
Si f est continue la suite B_n(f) converge uniformément vers f dans I=[0,1].

13.4.2  Le programme

bernstein(f,n,t) approche uniformement f continue sur [0,1].
On tape :

bernstein(f,n,t):={
retourne sum(comb(n,p)*f(p/n)*(1-t)^(n-p)*t^p,p=0..n);
  }:;

bernab(f,n,t,a,b) approche uniformement f continue sur [0,1].
On tape :

bernab(f,n,t,a,b):={
 retourne  sum(comb(n,p)*f(a*(1-p/n)+b*p/n)*(b-t)^(n-p)*(t-a)^p/(b-a)^n,p=0..n)
   }:;

On tape :
plotfunc([bernstein(sin,12,x),sin(x)],x)
On obtient :
Un graphe proche de sin(x) sur 0,1
On tape :
plotfunc([bernstein(sin,12,x,-pi/2,pi/2),sin(x)],x)
On obtient :
Un graphe proche de sin(x) sur -pi/2,pi/2

13.5  Développements asymptotiques et séries divergentes

Un développement asymptotique est une généralisation d’un développement de Taylor, par exemple lorsque le point de développement est en l’infini. De nombreuses fonctions ayant une limite en l’infini admettent un développement asymptotique en l’infini, mais ces développements sont souvent des séries qui semblent commencer par converger mais sont divergentes. Ce type de développement s’avère néanmoins très utile lorsqu’on n’a pas besoin d’une trop grande précision sur la valeur de la fonction.

13.5.1  Un exemple:la fonction exponentielle intégrale

Nous allons illustrer ce type de développement sur un exemple, la fonction exponentielle intégrale, définie à une constante près par

On peut montrer que l’intégrale existe bien, car l’intégrand est positif et inférieur à e^−t (qui admet −e^−t comme primitive, cette primitive ayant une limite en +∞).
Pour trouver le développement asymptotique de f en +∞, on effectue des intégrations par parties répétées, en intégrant l’exponentielle et en dérivant la fraction rationnelle :

où

Le développement en séries est divergent puisque pour x>0 fixé et n tendant vers l’infini

mais si x est grand, au début la série semble converger, de manière très rapide :

On peut utiliser S(x) comme valeur approchée de f(x) pour x grand si on sait majorer R(x) par un nombre suffisamment petit. On a

On retrouve une majoration du type de celle des séries alternées, l’erreur est inférieure à la valeur absolue du dernier terme sommé. Pour x fixé assez grand, il faut donc de trouver un rang n, s’il en existe un, tel que n!/xⁿ<є où є est la précision relative que l’on s’est fixée. Par exemple, si x≥ 100, n=12 convient pour є=12!/100¹²=5e−16 (à peu près la précision relative d’un “double”).

13.5.2  Le calcul approché de la constante d’Euler γ

Pour d’autres méthodes concernant le calcul approché de la constante d’Euler voir aussi 13.3.4 et 15.2.3.
On peut montrer que

existe (par exemple en cherchant un équivalent de u_n+1−u_n qui vaut −1/2n²) et on définit γ comme sa limite. Malheureusement, la convergence est très lente et cette définition n’est pas applicable pour obtenir la valeur de γ avec une très grande précision. Il y a un lien entre γ et la fonction exponentielle intégrale (définie par f(x)=∫_x^+∞ e^−t/t dt) , plus précisément lorsque x→ 0, f(x) admet −ln(x) comme singularité, plus précisément f(x)+ln(x) admet un développement en séries (de rayon de convergence +∞), car :

Que vaut la constante du membre de droite :

Il se trouve que C=−γ (voir plus bas une démonstration condensée) et donc :

Pour obtenir une valeur approchée de γ, il suffit donc de prendre un x assez grand pour pouvoir calculer f(x) par son développement asymptotique à la précision requise (f(x)S(x)+R(x) avec S(x)=e^−x(1/x−1/x² + 2/x³+...+(−1)ⁿ n!/xⁿ⁺¹) et R(x)=− ∫_x^+∞ (−1)ⁿ n!e^−t/tⁿ⁺¹ dt et |R(x)|≤ n!e^−x/xⁿ⁺¹), puis de calculer l’intégrale du membre de droite par le développement en séries en x=0 (en utilisant une précision intermédiaire plus grande puisque ce développement en séries va sembler diverger au début avant de converger pour n suffisamment grand).
Exemple1 : on pose x=13.
On calcule f(13) par (1) avec n=13 et une erreur absolue inférieure à e⁻¹³ 13!/13¹⁴≤ 3.6e−12.
On a en effet pour x=x₀ si v_n=e^−x₀n!/x₀ⁿ⁺¹ :
v_n/v_n−1≤ 1 équivalent à n/x₀≤ 1équivalent à n≤ x₀.
Donc si x=x₀=13 on calcule f(13) avec :
f(13)≃ ∑_n=0¹³ e⁻¹³(−1)ⁿ n!/xⁿ⁺¹ Ou bien, on tape :
Digits:=2; puis
exp(-13)*n!/13.^(n+1))$(n=0..20) renvoie :
1.7e-07,1.3e-08,2.1e-09,4.7e-10,1.5e-10,5.6e-11,2.6e-11,
1.4e-11,8.6e-12,5.9e-12,4.6e-12,3.9e-12,3.6e-12, 3.6e-12,
3.8e-12,4.4e-12,5.5e-12,7.1e-12,9.9e-12,1.4e-11,2.2e-11
donc

puis on remplace dans (3), avec

dont on obtient une valeur approchée, en faisant la somme jusqu’au rang 49, le reste de cette somme R₅₀ est positif et est inférieur à 13.^51/51/51!) qui est de l’ordre de 8.2e-12. On a en effet si v_n=13ⁿ⁺¹/(n+1)(n+1)! :
|R_N|=∑_n=N+1^∞v_n<v_n+1=13^N+2/(N+2)(N+2)! et |R₄₉|<8.2e−12.

La somme argument de evalf étant exacte, il n’y a pas de problèmes de perte de précision.
On obtient finalement comme valeur approchée de γ

On choisit alors 12 chiffres significatif et on tape :
Digits:=12; f13:=exp(-13.)*evalf(sum((-1)^n*n!/13^(n+1),n=0..13)) I13:=evalf(sum((-1)^*13^(n+1)/(n+1)/(n+1)!,n=0..49))
La constante d’Euler vaut donc à 1.2e-11 près :
-f13-ln(13.)+I13
On obtient :
0.577215664897
On tape :
evalf(euler_gamma)
On obtient :
0.5772156649018
soit 0.57721566489 avec une erreur inférieure à 1.2e-11.
Exemple2 : on pose x=40.
On tape :
r40:=(exp(-40.)*40!/40.^41) on obtient r40 inférieur à 7.2e-36
On choisit alors 36 chiffres significatif et on tape :
Digits:=36;
f40:=exp(-40.)*evalf(sum((-1)^n*n!/40^(n+1),n=0..40))
puisque :
40.^168/168./168! est inférieur à 3.3e-36, on tape :
I40:=evalf(sum((-1)^n*40^(n+1)/(n+1)/(n+1)!,n=0..166))
La constante d’Euler vaut donc à (7.2+3.3)e-36 près:
-f40-ln(40.)+I40
On obtient avec une erreur inférieure à 1.1e-35:
0.5772156649015328606065120900824024285
On tape :
evalf(euler_gamma)
On obtient :
0.5772156649015328606065120900824024308
Remarques La somme argument de evalf étant exacte, il n’y a pas de problèmes de perte de précision, on peut aussi faire les calculs intermédiaires en arithmétique approchée, lorsque x=13 on doit alors prendre 4 chiffres significatifs de plus (pour tenir compte de la valeur du plus grand terme sommé dans la série v_n=13ⁿ⁺¹/(n+1)(n+1)! qui est v₁₀=13^11/11/11!\simeq 4.08e+03).
On a en effet :
v_n/v_n−1=13*n/(n+1)²>1 si n²+11n+1>0 i.e. n≤ 10 et
v_n/v_n−1=13*n/(n+1)²<1 si n≥ 11
On tape avec des calculs intermédiaires en arithmétique approchée :

On obtient dans ce cas comme valeur approchée de γ : 0.57721566489675213
Bien entendu, cette méthode avec des calculs intermédiaires en arithmétique approchée est surtout intéressante si on veut calculer un grand nombre de décimales de la constante d’Euler c’est à dire quand on prend x=x0 tres grand, sinon on peut par exemple appliquer la méthode d’accélération de Richardson (cf 15.2.3) à la suite convergente (2) qui définit γ.
On peut calculer π de la même manière avec le développement en séries et asymptotique de la fonction sinus intégral (on remplace exponentielle par sinus dans la définition de f, voir plus bas une démonstration condensée) et l’égalité

Calcul de C (et preuve de (4)) :
Pour cela on effectue une intégration par parties, cette fois en intégrant 1/t et en dérivant l’exponentielle (moins 1 dans la première intégrale).

Pour calculer cette intégrale, on utilise l’égalité (qui se démontre par récurrence en faisant une intégration par parties) :

On va à nouveau intégrer par parties, on intègre un facteur 1 et on dérive l’intégrand, on simplifie, puis on intègre t et on dérive l’autre terme, puis t²/2, etc.

où

Pour déterminer I_n on fait le changement de variables t=nu

Or en faisant le même changement de variables t=nu :

Donc

Lorsque n tend vers l’infini, on peut montrer que I_n → 0, en effet les intégrales sont équivalentes à leur valeur sur un petit intervalle autour de u=1, point où l’argument de l’exponentielle est maximal, et comme l’intégrand du numérateur a une amplitude ln(u) qui s’annule en u=1, il devient négligeable devant le dénominateur. Finalement on a bien C=−γ.
On peut remarquer qu’en faisant le même calcul que C mais en remplacant e^−t par e^{−α t} pour ℜ(α)>0, donne limI_n=−ln(α) (car le point critique où la dérivée de la phase s’annule est alors 1/α). Ceci peut aussi se vérifier pour α réel en faisant le changement de variables α t=u

En faisant tendre α vers −i, −ln(α) tend vers ln(i)=iπ/2 et on obtient

dont la partie imaginaire nous donne (4), et la partie réelle une autre identité sur γ faisant intervenir la fonction cosinus intégral.

13.6  Solution de f(x)=0 par la méthode de Newton

Dans Xcas, il existe déjà une fonction qui calcule la valeur approchée r d’une solution de f(x)=0 par la méthode de Newton, qui est : newton.

13.6.1  La méthode de Newton

Soit f deux fois dérivable ayant un zéro et un seul r dans l’intervalle [a ; b]. Supposons de plus que f′ et f″ ont un signe constant sur [a ; b]. La méthode de Newton consiste à approcher r par l’abscisse x₁ du point commun à Ox et à la tangente en un point M₀ du graphe de f. Si M₀ a pour coordonées (x₀,f(x₀)) (x₀ ∈ [a ; b]), la tangente en M₀ a pour équation :
y=f(x₀)+f′(x₀)*(x−x₀) et donc on a :

On peut alors réitérer le processus, et on obtient une suite x_n qui converge vers r soit par valeurs supérieures, si f′*f″>0 sur [a ; b] (i.e. si f′(r)>0 et si f est convexe (f″>0 sur [a ; b]) ou si f′(r)<0 et si f est concave (f″<0 sur [a ; b])) soit par valeurs inférieures, si f′*f″<0 sur [a ; b] (i.e. si f′(r)<0 et si f est convexe (f″>0 sur [a ; b]) ou si f′(r)>0 et si f est concave (f″<0 sur [a ; b])).
On fait le dessin en tapant :

f(x):=x*x-2;
x0:=5/2;
G:=plotfunc(f(x));
T0:=tangent(G,x0);
Ox:=droite(0,1);
M1:=inter(T0,Ox)[0];
x1:=affixe(M1)
segment(point(x1,0),point(x1,f(x1)));
T1:=tangent(G,x1);
M2:=inter(T1,droite(0,1))[0]
x2:=affixe(M2):
segment(point(x2,0),point(x2,f(x2)));

ou encore pour faire le dessin de la méthode de Newton pour la fonction f en partant du point de coordonnèes (a,f(a)) et obtenir p nouveaux points.

plotnewton(f,a,p):={
local L,P,m,j,b;
L:=plotfunc(f(x),x,affichage=vert);
L:=L,point(a,couleur=point_width_4+rouge);
for (j:=1;j<=p;j++) {
b:=f(a);
L:=L,segment(a,a+i*b,couleur=ligne_tiret+rouge);
m:=function_diff(f)(a);
L:=L,plotfunc(b+(x-a)*m,x);
if (m==0){return "pente nulle"}
a:=a-f(a)/m;
P:=point(a,couleur=point_width_4+rouge);
L:=L,P;
}
return affixe(P),L;
};

On tape : plotnewton(sq-2,4,2)
pour obtenir les termes x₀,x₁,x₂ de la suite de Newton qui converge vers √2 et où x₀=4 :

On remarquera que : plotnewton(sq-2,4,2)[0] renvoie :
113/72 ≃ 1.56944444444
On tape : plotnewton(ln-1,5,2)
pour obtenir les termes x₀,x₁,x₂ de la suite de Newton qui converge vers e et où x₀=5 :

On peut aussi faire une animation, pour cela, on tape :

newtonsuite(f,a,p):={
local L,P,m,j,b,LT;
P:=point(a,couleur=point_width_4+rouge);
LT:=P;
f1:=function_diff(f);
for (j:=1;j<=p;j++) {
b:=f(a);
L:=L,segment(a,a+i*b,couleur=ligne_tiret+rouge);
m:=f1(a);
L:=L,plotfunc(b+(x-a)*m,x);
if (m==0){return "pente nulle"}
a:=a-f(a)/m;
P:=point(a,couleur=point_width_4+rouge);
LT:=LT,[LT,L,P];
}
print(affixe(P));
return LT;
};
animnewton(f,a,p):={
local LT;
LT:=newtonsuite(f,a,p);
return plotfunc(f(x),x,affichage=vert),animation(LT);
};

On tape : animnewton(sq-2,4,3)
Puis, on écrit la fonction newton_rac qui renvoie la valeur approchée à eps près de la racine de f(x)=0 on commençant l’itération avec x0.
On remarquera que le paramètre f est une fonction et donc, que sa dérivée est la fonction g:=function_diff(f).
On cherche une valeur approchée donc il faut écrire :
x0:=evalf(x0-f(x0)/g(x0))
car si on ne met pas evalf, les calculs de l’itération se feront excactement et seront vite compliqués.

newton_rac(f,x0,eps):={
local x1,h,g;
g:=function\_diff(f)
x0:=evalf(x0-f(x0)/g(x0));
x1:=x0-f(x0)/g(x0);
if (x1>x0) {h:=eps;} else {h:=-eps;}
while (f(x1)*f(x1+h)>0){
   x1:=x1-f(x1)/g(x1);
}
return x1;
}

13.6.2  Exercices sur la méthode de Newton

1. L’énoncé
   1. Étudier rapidement les variations de f(x)=xexp(x)+0.2 pour montrer que l’équation f(x)=0 a deux solutions a et b qui vérifient a<−1<b<0
   2. Calculer à l’aide d’un programme par la méthode de Newton, les valeurs approchées de a et de b obtenues après 5 itérations.
   3. Modifier votre programme pour avoir des valeurs de a et b avec une précision de eps (par exemple de 1e−6).
   4. Écrire un programme qui dessine le graphe de la fonction implicite xexp(x)+yexp(y)=0 lorsque x≥−5.

   La solution avec Xcas

   1. Pour avoir les variations de f(x)=xexp(x)+0.2 on peut calculer la dérivée et faire le graphe de f.
      On tape : f(x):=x*exp(x)+0.2
      factor(diff(f(x)))
      On obtient : (x+1)*exp(x)
      On tape : G:=plotfunc(f(x),x=-3..1)
      On obtient :
      

      Pour montrer que l’équation f(x)=0 a deux solutions a et b qui vérifient a<−1<b<0, on calcule, d’après le graphe f(−3),f(−2),f(−1),f(0).
      On tape : f(x)$(x=-3..0)
      On obtient :
      0.0506387948964,-0.0706705664732,-0.167879441171,0.2
      donc puisque f est continue, d’après le théorème des valeurs intermédiaires on a : −3<a<−2 et −1<b<0
      

   2. La méthode de Newton consiste à itérer la fonction h définie par h(x)=x−f(x)/g(x) où g est la derivée de f. Pour trouver a, on va commencer en x₀=−2 (car la fonction est concave et décroissante sur [−3;−2] et les x_n seront des valeurs approchées de a par excés) et pour trouver b, on va commencer en x₀=0 (car la fonction est convexe et croissante sur [−1;0] et les x_n seront des valeurs approchées de b par excés). On tape :

      Newtonvaleur(x0):={
        local j,f,g,h;
        f(x):=x*exp(x)+0.2;
        g(x):=(x+1)*exp(x);
        h(x):=x-f(x)/g(x);
        pour j de 1 jusque 5 faire
          x0:=h(x0);
        fpour;
        retourne x0;
      }:;
      

      Remarque
      On peut remplacer g(x):=(x+1)*exp(x); par g:=function_diff(f);
      On tape pour avoir la valeur de x₅ qui approche a :
      Newtonvaleur(-2.)
      On obtient :
      -2.54264135777
      On tape pour avoir la valeur de x₅ qui approche b :
      Newtonvaleur(0)
      On obtient :
      -0.259171101819
   3. Si on veut avoir une valeur approchée de a (resp b) à eps prés, il faut avoir un x_j qui vérifie x_j−eps≤ a<x_j (resp x_j−eps ≤ b<x_j) c’est à dire f(x_j−eps)>0 (resp f(x_j−eps)<0) i.e. on doit avoir dans les 2 cas, f(x_j−eps)*f(x0)≤ 0.
      On tape :

      Newtonvalpres(x0,eps):={
        local j,g,h,t,s;
        f(x):=x*exp(x)+0.2;
        g(x):=(x+1)*exp(x);
        h(x):=x-f(x)/g(x);
        j:=0;
        t:=x0-eps;
        //s:=ifte(f(x0)>0,1,-1);
        s:=sign(f(x0));
        tantque s*f(t)>0 faire
          x0:=h(x0);
          t:=x0-eps;
          j:=j+1;
        ftantque;
        print(j);
        retourne t,x0;
      }:;
      

      On tape pour avoir la valeur de x_j qui donne un encadrement de a à 1e−6 prés :
      Newtonvalpres(-2.,1e-6)
      On obtient pour j=3 :
      -2.54264235686,-2.54264135686
      On tape pour avoir la valeur de x_j qui donne un encadrement de a à 1e−10 prés :
      Newtonvalpres(-2.,1e-10)
      On obtient pour j=4 :
      -2.54264135787,-2.54264135777
      On tape pour avoir la valeur de x_j qui donne un encadrement de b à 1e−6 prés :
      Newtonvalpres(0,1e-6)
      On obtient pour j=4 :
      -0.259172101477,-0.259171101477 On tape pour avoir la valeur de x_j qui donne un encadrement de b à 1e−10 prés :
      Newtonvalpres(0,1e-10)
      On obtient pour j=5 :
      -0.259171101919,-0.259171101819
   4. yexp(y)+a=0 a une solution si −exp(−1)+a ≤ 0 i.e si a<1/e Cette fonction est définie pour des x tels que xexp(x)=a<1/e. On résout donc l’équation xexp(x)−1/e=0:
      On modifie le programme en :
      

      Newtonvaleura(x0,a):={
        local j,f,g,h;
        f(x,a):=x*exp(x)+a;
        g(x):=(x+1)*exp(x);
        h(x):=x-f(x,a)/g(x);
        pour j de 1 jusque 5 faire
          x0:=h(x0);
        fpour;
        retourne x0;
      }:;
      

      lorsque a=−1/e, on a f(0,a)=a<0 et f(1,a)=e−1/e>0. On tape :
      Newtonvaleura(1,-1./e)
      On obtient :
      0.278464542823
      Pour x≤ 0, a=x*exp(x)≤ 0 donc l’équation en y, y*exp(y)+x*exp(x)=0 a une seule solution positive (y*exp(y)+a vaut a≤ 0 pour y=0 et vaut (x(exp(x)²−1)/(exp(x))>0 pour y=−x). On peut donc l’obtenir par la méthode de Newton : on démarre avec y₀=−a car pour x>0 la courbe de f(x,a)=x*exp(x)+a se trouve au dessus de sa tangente en x=0 (puisque f"(x)>0 pour x>0) et que cette tangente d’équation y=x+a coupe l’axe des x en x=−a.
      Pour 0<x<0.278464542823 l’équation en y, y*exp(y)+x*exp(x)=0 a deux solutions que l’on peut obtenir par la méthode de Newton : on démarre soit par x₀=0 soit par x₀=−2.
      On tape (on choisit de prendre x₅ comme valeur approchée) :

      Newtonimplicit():={
        local j,f,g,h,a,xj,y0,y,L;
        g(y):=(y+1)*exp(y);
        f(y,a):=y*exp(y)+a;
        pour xj de -4 jusque 0 pas 0.1 faire
        a:=xj*exp(xj);
        h(y):=y-f(y,a)/g(y);
        y0:=-a;
        pour j de 1 jusque 5 faire
          y0:=h(y0);
        fpour;
        L:=L,point(xj+i*y0);
      fpour;
      L:=L,point(0.28-i);
      pour xj de 0.01 jusque 0.28 pas 0.02 faire
        a:=xj*exp(xj);
        h(y):=y-f(y,a)/g(y);
        y0:=0;
        pour j de 1 jusque 5 faire
          y0:=h(y0);
        fpour;
        L:=L,point(xj+i*y0);
        y0:=-2;
        pour j de 1 jusque 5 faire
          y0:=h(y0);
        fpour;
        L:=L,point(xj+i*y0)
      fpour;
      retourne L;
      }:;
      

      On tape Newtonimplicit()
      On obtient :
      

      On peut aussi vouloir calculer y à eps -près. Mais attention lorsqu’on part de y0=−2 on obtient une valeur soit par défaut, soit par excés selon le signe de f(−2,a) (si f(−2,a)>0 ce sera par excès car pour x=−2 on a un point d’inflexion).
      On tape alors :

      Newtonimpl(eps):={
        local j,f,g,h,a,xj,y0,y,t,s,L;
        g(y):=(y+1)*exp(y);
        f(y,a):=y*exp(y)+a;
        L:=NULL;
        pour xj de -5 jusque 0 pas 0.05 faire
        a:=evalf(xj*exp(xj));
        h(y):=y-f(y,a)/g(y);
        y0:=-a;
        t:=y0-eps;
        s:=sign(f(y0,a));
        tantque s*f(t,a)>0 faire
          y0:=h(y0);
          t:=y0-eps;
        ftantque;
        L:=L,point(xj+i*y0);
      fpour;
      L:=L,point(0.28-i);
      pour xj de 0.01 jusque 0.28 pas 0.02 faire
        a:=evalf(xj*exp(xj));
        h(y):=y-f(y,a)/g(y);
        y0:=0.;
        t:=y0-eps;
        s:=sign(f(y0,a));
        tantque s*f(t,a)>0 faire
          y0:=h(y0);
          t:=y0-eps;
        ftantque;
        L:=L,point(xj+i*y0);
      fpour;
        pour xj de 0.01 jusque 0.28 pas 0.02 faire
          a:=evalf(xj*exp(xj));
        h(y):=y-f(y,a)/g(y);
        y0:=-2.;
        s:=sign(f(y0,a));
        si s>0  alors eps:=-abs(eps); fsi; 
        t:=y0-eps;
        tantque s*f(t,a)>0 faire
          y0:=h(y0);
          t:=y0-eps;
        ftantque;
        L:=L,point(xj+i*y0)
      fpour;
      retourne L;
      }:;
      

      On tape Newtonimpl(0.01)
      On peut vérifier en tapant :
      plotimplicit(x*exp(x)+y*exp(y),[x,y])

2. L’énoncé
   Donner la valeur approchée de cos(x)=x pour x∈ [0;1] obtenue en partant de x₀=0 après 4, puis après 10 itérations lorsque :
   1. On applique la méthode du point fixe à f(x)=cos(x).
   2. On applique la méthode de Newton.
   3. On applique la méthode du Δ² d’Aitken.
   4. On applique la méthode de Steffensen.
   Quelle méthode vous semble la meilleure ? Expliquez pourquoi.

   La solution avec Xcas
   On configure Xcas avec 20 digits.

   1. fa fonction f(x)=cos(x) est sin(1)-contractante sur [0;1], car d’après le théorème des accroissements finis :
      il existe c∈ [0;1] tel que pour tout x₁∈ [0;1] et tout x₂∈ [0;1] on a cos(x₁)−cos(x₂)=(x₁−x₂)sin(c) donc
      |cos(x₁)−cos(x₂)|=|x₁−x₂||sin(c)|≤ |x₁−x₂|sin(1).
      On tape :

      ptfixecos(x0,n):={
      local j,f;
      f(x):=cos(x);
      pour j de 1 jusque n faire
        x0:=f(evalf(x0));
      fpour;
      retourne x0;
      }:;
      

      On tape : ptfixecos(0,4)
      On obtient : 0.654289790497779149974
      On tape : ptfixecos(0,10)
      On obtient : 0.731404042422509858293
      
   2. On pose F(x)=cos(x)−x et on a F′(x)=−sin(x)−1, on va donc itérer la fonction g(x)=x−F(x)/F′(x)=(x*sin(x)+cos(x))/(sin(x)+1).
      On tape :

      Newtoncos(x0,n):={
      local j,g,F,dF;
      F(x):=cos(x)-x;
      dF:=function_diff(F);
      //g(x):=(x*sin(x)+cos(x))/(sin(x)+1);
      g(x):=normal(x-F(x)/dF(x));
      pour j de 1 jusque n faire
        x0:=g(evalf(x0));
      fpour;
      retourne x0;
      }:;
      

      On tape : Newtoncos(0,4)
      On obtient : 0.739085133385283969762
      On tape : Newtoncos(0,10)
      On obtient : 0.739085133215160641651
      
   3. La méthode du Δ² d’Aitken consiste à transformer la suite des itérées du point fixe par la fonction :
      gs(x)=x−(f(x)−x)*(f(x)−x)/(f(f(x))−2f(x)+x) avec f(x)=cos(x).
      On tape :

      Aitkencos(x0,n):={
      local j,gs,f,y0;
      f(x):=cos(x);
      gs(x):=x-(f(x)-x)*(f(x)-x)/(f(f(x))-2f(x)+x);
      pour j de 1 jusque n faire
        x0:=f(evalf(x0));
        y0:=gs(x0);
      fpour;
      print(x0);
      retourne y0;
      }:;
      

      On tape : Aitkencos(0,4)
      On obtient : 0.738050421371663847259
      On tape : Aitkencos(0,10)
      On obtient : 0.739076383318955862683
      
   4. La méthode de Steffenson consiste à itérer la fonction :
      gs(x)=x−(f(x)−x)*(f(x)−x)/(f(f(x))−2f(x)+x) avec f(x)=cos(x).
      On tape :

      Steffensencos(x0,n):={
      local j,gs,f;
      f(x):=cos(x);
      gs(x):=x-(f(x)-x)*(f(x)-x)/(f(f(x))-2f(x)+x);
      pour j de 1 jusque n faire
        x0:=gs(evalf(x0));
      fpour;
      retourne x0;
      }:;
      

      On tape : Steffensencos(0,4)
      On obtient : 0.739085133215160534355
      On tape : Steffensencos(0,10)
      On obtient : 0.739085133215160641651
      

   Les méthodes de Newton et de Steffensen sont plus performantes car ce sont des méthodes d’ordre 2 (la fonction que l’on iére a une dérivée nulle au point solution de f(x)=cos(x)=x).
   Même avec Digits:=30 on a :
   Steffensencos(0,6)=Steffensencos(0,10)=
   Newtoncos(0,6)=Newtoncos(0,10)=
   0.7390851332151606416553120876735

3. L’énoncé
   Dans un problème de baccalauréat, on se propose de calculer des valeurs approchées des solutions de exp(−x)cos(x)=x sur [−π/2;π/2].
   1. Déterminer le nombre de solutions.
   2. Donner un algorithme de calcul et écrire le programme correspondant.
   3. Donner un encadrement de la solution ?
   4. Dessiner les points de coordonnèes t,x qui vérifient :
      exp(−x)cos(x)−x+t=0 pour t ∈ [−π/2;π/2] et x ∈ [−π/2;π/2]
   La solution avec Xcas
   
   1. On tape : f(x:)=exp(-x)*cos(x)-x
      f1:=function_diff(f);f2:=normal(diff(f1(x))) On trouve :
      f1(x)=(-cos(x)-sin(x))*exp(-x)-1 et
      f2=2*exp(-x)*sin(x)
      On a :
      f2=0 en x=0 le point(0,1) est donc un point d’inflexion et
      sur [−π/2;0], on a f2<0 et sur [0;π/2], on a f2>0
      On a :
      f1(x)=exp(π/2)−1>0 pour x=−π/2
      f1(0)=−2
      f1(x)=−exp(−π/2)−1 <0 pour x=π/2
      f1 s’annule donc pour x=a<0 et donc f est croissante puis décroissante et comme on a f(−pi/2)=pi/2, f(0)=1 et f(pi/2)=−pi/2 et f continue, f s’annule en un seul point x=b>0 Pour vérifier, on tape :
      G:=plotfunc(f(x),x=-pi/2..pi/2);tangente(G,0)
      
   2. On peut appliquer la méthode de Newton en partant de x₀=0.0, la suite x_n=x_n−1−f(x_n−1)/f1(x_n−1) va donner une valeur approchée par défaut de b car sur [0;b] f est positive décroissante et convexe.

      On tape :
      Newton0(4)
      On obtient avec 22 digits:
      0.51775736368245829829471
      

   3. On tape :

      Newton0(n):={
      local j,f,f1,g,x0;
      f(x):=exp(-x)*cos(x)-x;
      f1:=function_diff(f);
      g(x):=normal(x-f(x)/f1(x));
      x0:=0.0;
      pour j de 1 jusque n faire
      x0:=g(x0)
      fpour;
      retourne x0;
      }:;
      

   4. Pour avoir un encadrement de la solution à eps près, on continue l’itération tant que la valeur de f(x_n+eps) reste strictement positive.
      On tape

      Newtoneps(n,eps):={
        local j,f,f1,g,x0;
        f(x):=exp(-x)*cos(x)-x;
        f1:=function_diff(f);
        g(x):=normal(x-f(x)/f1(x));
        x0:=0.0;
        j:=0;
        tantque f(x0+eps)>0 faire
          x0:=g(x0);
          j:=j+1;
        ftantque;
        print(j);
        retourne x0,x0+eps;
      }
      :;
      

      On tape :
      Newtoneps(0,1e-20)
      On obtient :
      0.51775736368245829832277,0.51775736368245829833277
      
   5. Il faut voir que si on remplace f(x) par f(x)+t la suite définie par :x₀=0 et x_n+=x_n−(f(x_n)+t)/f1(x_n) approche la solution en x de f(x)+t=0 par excès si t<−1 et par défaut si t>−1 car pour x=0 la fonction f(x)+t a un point d’inflexion qui est le point (0;1+t). De plus f(x)+t>1+t pour x<0 et f(x)+t<1+t pour x>0. Donc si 1+t>0 la solution sera positive et si 1+t<0 la solution sera negative
      On tape :

      Newtonimpl():={
        local j,f,f1,g,x0,t,a,L;
        a:=evalf(pi/2);
        f(x):=exp(-x)*cos(x)-x;
        f1:=function_diff(f);
        L:=NULL;
          pour t de -a jusque -1 pas 0.1 faire
            g(x):=normal(x-(f(x)+t)/f1(x));
            x0:=0.0;
            tantque f(x0-0.01)+t<0 faire
              x0:=g(x0);
            ftantque;
            L:=L,point(t,x0);
          fpour;
          pour t de -1 jusque a pas 0.1 faire
            g(x):=normal(x-(f(x)+t)/f1(x));
            x0:=0.0;
            tantque f(x0+0.01)+t>0 faire
              x0:=g(x0);
            ftantque;
            L:=L,point(t,x0);
          fpour;
          return L;
        }
      :;
      

      On tape :
      Newtonimpl()
      On obtient :
      

13.6.3  La méthode de Newton avec préfacteur

Lorsqu’on part d’une valeur x₀ trop éloignée de la racine de f(x) (si par exemple |f(x₀)| est grand), on a intérêt à utiliser un préfacteur pour se rapprocher plus vite de la solution de f(x)=0.
Posons n(x)=−f(x)/f′(x), on a alors :

donc

ce qui veut dire que :
f(x₀+h*n(x₀))=f(x₀)(1−h)+h· є(h) avec є(h) tend vers 0 quand h tend vers 0.
Donc, il existe h₀ vérifiant :
|f(x₀+h₀*n(x₀))|<|f(x₀)|
Remarque : Il faut minimiser |f(x₀+h₀*n(x₀))|. or plus h₀ est proche de 1 et plus |f(x₀)*(1−h₀)| sera petit. Par exemple, on prendra le plus grand h₀, dans la liste [1,3/4,(3/4)²,...] qui vérifie |f(x₀+h₀*n(x₀))|<|f(x₀)|
Pour cette valeur de h₀, x₀+h₀*n(x₀) est probablement plus proche de la racine que x₀ : on dit que h₀ est le préfacteur de la méthode de Newton.
On va choisir par exemple au début h₀=1, et on regarde si |f(x₀+n(x₀))|<|f(x₀)|, si ce n’est pas le cas on prend h₀=(3/4) et on regarde si |f(x₀+3/4*n(x₀))|<|f(x₀)|, si ce n’est pas le cas on prend h₀=(3/4)² etc...
On change de préfacteurs à chaque étape jusqu’à ce que : abs(f(x1))−abs(f(x0))<0 sans préfacteur, on continue alors l’itération sans préfacteur, c’est à dire avec la méthode de Newton normale. On écrit donc :

newton_prefacts(f,x0,eps):={
local x1,h,h0,prefact,niter;
//prefact est egal par ex a 3/4
h0:=1.0;
niter:=0;
prefact:=0.75;
x1:=x0-h0*f(x0)/function_diff(f)(x0);
while (abs(f(x1))-abs(f(x0))>0) {
   h0:=h0*prefact;
   x1:=x0-h0*f(x0)/function_diff(f)(x0);
} 
h:=eps;
while (h0!=1 and niter<100){
   x0:=x1;
   x1:=x1-h0*f(x1)/function_diff(f)(x1);
   while (abs(f(x1))-abs(f(x0))>0) {
      h0:=h0*prefact;
      x1:=x0-h0*f(x0)/function_diff(f)(x0);
   } 
   while (abs(f(x1))-abs(f(x0)<0 and h0!=1)) {
      h0:=h0/prefact;
      x1:=x0-h0*f(x0)/function_diff(f)(x0);
   } 
niter:=niter+1;   
}
while (f(x1-h)*f(x1+h)>0  and niter<200){
   x0:=x1;
   x1:=x1-f(x1)/function_diff(f)(x1);
   niter:=niter+1;
}

if  (niter<200) {return x1;} else {return "pas trouve";} 
}

On définit la fonction f par f(x):=x²-2 et on met ce programme dans un niveau éditeur de programmes (que l’on ouvre avec Alt+p), puis on le teste et on le valide avec OK.
On tape :
newton_prefacts(f,100,1e-10)
On obtient :
1.41421356237 On tape :
newton_prefacts(f,3,1e-5)
On obtient :
1.41421378005

13.7  Trouver un encadrement de x0 lorsque f(x0) est minimum

Soit f une fonction définie sur [a;b]. On suppose que f est unimodale sur [a;b], c’est à dite que f a un seul extremum sur [a;b]. On suppose de plus que cet extremum est un minimum (sinon on remplacera f par −f.) On se propose de trouver un encadrement à eps près de la valeur pour laquelle f est minimum.

13.7.1  Déscription du principe de la méthode

On partage [a;b] en trois morceaux en considérant c et d vérifiant : a<c<d<b.
On calcule f(c) et f(d) et on les compare.
Puisque f a un seul minimum sur [a;b] elle décroit, passe par son minimum, puis f croit. Selon les trois cas possibles on a :

• f(c)<f(d)
  dans ce cas, la valeur pour laquelle f est minimum n’appartient pas à [d;b]
• f(c)>f(d)
  dans ce cas, la valeur pour laquelle f est minimum n’appartient pas à [a;c]
• f(c)=f(d)
  dans ce cas, la valeur pour laquelle f est minimum n’appartient pas à [d;b] ni à [a;c]

Ainsi, l’intervalle de recherche a diminué et on peut recommencer le processus. Pour que l’algorithme soit performant, on veut que l’intervalle de recherche diminue rapidement et que le nombre de valeurs de f à calculer soit le plus petit possible. Pour cela comment doit-on choisir c et d ?

13.7.2  Déscription de 2 méthodes

On fait presque une dichotomie

On choisit c et d proche de a+b/2 par exemple :
c=a+b−eps/2 et d=a+b+eps/2 pour eps donné. Dans ce cas, à chaque étape l’intervalle diminue presque de moitié mais on doit calculer, à chaque étape, deux valeurs de f.

On utilise la suite de Fibonacci

Comment faire pour que l’une des valeurs de f déjà calculée serve à l’étape suivante ?
La solution se trouve dans la suite de Fibonacci, suite définie par :
u₀=1, u₁=2, u_n=u_n−2+u_n−1 dont les premiers termes sont : 1,2,3,5,8,13,21,34,55,89... Par exemple si on partage [a;b] en 89 parties égales si l=(b−a)/89, on choisit c=a+34*l et d=a+55*l et ainsi on a :
c−a=34*l, d−c=21*l, b−d=34*l (car 89=55+34 et 34+21=55 puisque 21,34,55,89 sont des termes consécutifs de la suite de Fibonacci).
On calcule f(c) et f(d) puis on réduit l’intervalle en un intervalle de longueur (b−a)*55/89, par exemple si l’intervalle suivant est [a;d] et, si on recommence le processus, le point c est le futur point d.
Donc à chaque étape il suffit de calculer une seule valeur de f pour passer de l’intervalle [a;b] (proportionnel à u_n) à l’intervalle [a;d] ou [c;b] (proportionnel à u_n−1). Il y a bien sûr le cas f(c)=f(d) où il faut à l’étape suivante calculer deux valeurs de f, mais dans ce cas on gagne 3 étapes car on passe de l’intervalle [a;b] (proportionnel à u_n) à l’intervalle [c;d] (proportionnel à u_n−3).
Selon la valeur eps de la longueur de l’encadrement, on calcule k:=ceil(2*(b−a)/eps); et la première valeur t=u_n de la suite de Fibonacci supérieure à k. il faut alors diviser l’intervalle [a;b] en t parties égales. On applique alors plusieurs fois le processus et on s’arrête quand n=1, c’est à dire quand l’intervalle a été réduit à un intervalle de longueur 2*(b−a)/t qui est, grace au choix de t (t>k>2*(b−a)/eps) inférieur à eps.

13.7.3  Traduction Xcas de l’algorithme avec Fibonacci

//f(x):=2*x^4-10*x^3-4*x^2+100
//fibomin(f,1,5,0.000001)
//g(x):=2*x^4-10*x^3+4*x^2+100
//fibomin(g,1,5,1e-20)
//calcul la valeur du min d'une fonction ayant 
//un seul extrema sur [a,b]
fibomin(f,a,b,eps):={
  local c,d,F,k,n,t,g,h,l,fc,fd;
  if (a>b) {c:=a;a:=b;b:=c;}
  k:=ceil(2*(b-a)/eps);
  F:=1,2;
  n:=1;
  g:=1;
  t:=2;
  //construction de F=la suite de Fibonacci
  //h,g,t sont 3 termes consecutifs de F
  while (t<k) {
    n:=n+1;
    h:=g;
    g:=t;
    t:=h+g;
    F:=F,t;
  }
  l:=(b-a)/t;
  c:=a+h*l;
  d:=a+g*l;
  fc:=f(c);
  fd:=f(d);
  //on itere le processus et on s'arrete qd n=1
  while (n>1) {
    if (fc>fd) {
      a:=a+h*l;
      fc:=fd;
      t:=h;
      h:=g-h;
      g:=t;
      fd:=f(a+g*l);
      n:=n-1;
    }else{
      if (fc<fd) {
 b:=a+g*l;
 t:=h;
 h:=g-h;
 g:=t;
 fd:=fc;
 fc:=f(a+h*l);
 n:=n-1;
      }else{
 a:=a+h*l;
 b:=b-h*l;
 t:=g-h;
 g:=h-t;
 h:=t-g;
 fc:=f(a+h*l);
 fd:=f(a+g*l);
 n:=n-3;
      }
    }
  }
  return [a,b];
}

On tape :
f(x):=x^4-10
fibomin(f,-1,1,1e-10)
On obtient :
[(-1)/53316291173,1/53316291173]
On tape :
g(x):=2*x^4-10*x^3-4*x^2+100
fibomin(g,1,5,1e-10)
On obtient :
[86267571271/21566892818,86267571273/21566892818] On tape :
h(x):=2*x^4-10*x^3+4*x^2+100
fibomin(h,1,5,1e-10)
On obtient :
[74644573011/21566892818,74644573013/21566892818]