la fonction evalf semble donner des résultats pardoxaux !

[majak]

je voudrais faire évaluer les termes de la suite récurrente :
u(n+1)=u(n)*(n+1) -1 en partant de u(0)=e-2

Cette suite décroit vers 0
Ceci est vrai pour une valeur de u(0) rigoureusement égale à e=exp(1)
et devient totalement faux pour toute autre valeur de u(0),
aussi proche soit-elle de e. XCAS confirme plus ou moins bien
ces prévisions ce qui m'a amené à poster ce mail sur le forum de XCAS :

f:=exp(1)-2;
pour j de 2 jusque 30
faire f:=developper(f*j-1);
fpour;
evalf(f,20);

Les termes de la suite définie dans la boucle pour .......fpour
ci-dessus décroissent vers 0.
Leur évaluation en changeant la valeur finale de l'indice de boucle est
très surprenante.
Tout semble confirmer que la limite est 0, alors que de manière
aléatoire est dépendante
de la précision demandée dans l'évaluation décimale certains termes
deviennent très grands
en valeur absolue ( positivement ou négativement ), voir les exemples
ci-dessous :

1) la valeur donnée par l'exemple ci-dessus est de l'ordre de 10^18;
alors que les trois termes suivants ou précédents sont évalués égaux
à 0 à 10^-15 près.

2)Par contre le résultat est évaluer nul comme prévu à 10^-20 près

f:=exp(1)-2;
pour j de 2 jusque 30
faire f:=developper(f*j-1);
fpour;
evalf(f,20);

3) mais maintenant, c'est pour l'indice final 34 que rien ne va plus,
il est évaluer nul à 10^-15 près , mais devient -10^19 à 10^-20 près.

Est ce que quelqu'un peut m'expliquer ce qui me semble être paradoxal.
Ceci peut venir de ma mauvaise interprétation de la fonction evalf(),
et comment programmer quelque chose, rétablissant des résultats
plus cohérents.

Jacques MAROT

[parisse]

c'est normal, par exemple pour j=34, on obtient sans evalf
295232799039604140847618609643520000000*exp(1)-802525952794456998100845719437321898885
et on fait alors la différence de 2 réels de l'ordre de 0.8e39, comme ils sont très proches, on perd toute précision, donc 20 chiffres significatifs et on se retrouve avec une différence sans chiffre significatif correct dans la mantisse et avec un exposant de 1e19 (39-20). Bien sur il pourrait se produire par hasard qu'il y ait quelques chiffres non significatifs de la mantisse nuls.

[majak]

Merci de la réponse que je soupçonnais, car c'est comme cela que fonctionne toute évaluation avec une calculette,
mais comment expliquer que pour j=35, 36, 37 on retrouve bien le 0 attendu !

Jacques MAROT

[parisse]

je pense que la précision est un peu meilleure que prévue "par hasard" (il faudrait regarder plus précisément la représentation en base 2 pour comprendre, comme on peut le faire pour comprendre pourquoi 0.3-3*0.1 ne renvoie pas 0).

[majak]

Pour mieux comprendre le problème sur lequel je viens de tomber,
j'ai été amené à programmer la fonction suivante et et les résultats
aberrants que j'avais rencontrés viennent de ce qui est pour moi
une bizarrerie. Je ne suis qu'un utilisateur en phase de découverte
de XCAS , est-ce que cei viendrait de la représentation en système
binaire des nombres entiers lors de leur évaluation:


fb(n):={
local u,j;
u:=1;
pour j de 1 jusque n faire u:=j*u+1; fpour;
}

Le calcul de fb(23) par XCAS donne l'entier :70273067330330098091156 (résultat exact !)
Mais son évaluation avec une mantisse à 15 chiffres donne : 0.7027306733033007e23
dont le dernier chiffre me semble aberrant.

Merci d'éclairer mes tourments !

Jacques MAROT

[alb]

Bonjour,

je voudrais faire évaluer les termes de la suite récurrente :
u(n+1)=u(n)*(n+1) -1 en partant de u(0)=e-2

je n'ai pas bien compris si cette suite est convergente vers 0 pour u(0)=e-2 ou pour u(0)=e

f:=exp(1)-2;
pour j de 2 jusque 30
faire f:=developper(f*j-1);
fpour;
evalf(f,20);

Je ne sais pas si la fonction ci-dessous peut faire avancer mais il semblerait qu'il faut éviter d'évaluer la forme simplifiée de u(n)

Code : Tout sélectionner

suite_expo1(n):={
  local T,L,j,d;
  L:=NULL; 
  T:=e-2;
  pour j de 2 jusque n faire
    T:=(j)*T-1;
  fpour;
  pour d de 10 jusque 40 faire
    L:=L,[evalf(T,d),evalf(simplify(T),d)];
  fpour; 
  return[L];
}
:;

j'ai essayé suite_expo1(100) puis suite_expo1(200)
suite_expo1(1000) donne un crash

[alb]

Je suis allé un peu trop vite,les approximations sont incohérentes dans les deux colonnes
Quelle est l'idée de la méthode prouvant la convergence de u(n) ?

[majak]

Il s'agit bien de la suite u définie par récurrence : u(n+1)=(n+1)u(n) -1 .
Cette suite diverge pour toute valeur de u(1) autre que e-1,
auquel cas elle converge vers 0. Ce qu'aucune calculatrice
ou tableur pourrait laisser croire. Il faut jouer fin avec xcas
pour le vérifier.

2 méthodes pour prouver sa convergence vers 0

1)u(n) est l'intégrale de 0 à 1 de exp(x).(1-x)^n
la formule de récurrence vient d'une intégration par partie.
Sa convergence vers 0 est une application classique de
la convergence dominée ou monotone.

2)On peut prouver par récurrence que u(n) est égal au produit de n!
par la somme pour k variant de (n+1) à l' infini des ( 1/k! )
Ceci prouve que u(n) =1/(n+1) + v(n) où v(n) décroit vers 0.
mais xcas évalue u(n) égal à 0 pour des valeurs de n
beaucoup trop petites. Sauf de manière "alétoire" pour
certaines valeurs de n où le résultat est très éloigné de 0.
Il faut demander une évaluation de u(n) avec une mantisse
dont le nombre de chiffres est au moins égal à celui de n!
pour avoir un résultat cohérent, effectivement supérieur
est proche de 1/(n+1) .

En fait u(n) = a(n).e -b(n) ; où la fraction b(n)/a(n) est le rationnel
égal à la somme partielle d'indice n du développement classique
de e=exp(1), c'est à dire la somme des ( 1/ k!) pour k variant de 0 à n.

[majak]

ERREUR !!!
Cette suite diverge pour toute valeur de u(1) autre que :
***************
* u(1)=e-2 *
***************

[alb]

Avez-vous essayé de calculer u(n) = a(n).e-b(n) puis evalf(1+ln(a(n))-evalf(b(n)) ?
Pour écrire le programme il faudrait pouvoir renvoyer les arguments de A-B et je ne retrouve pas la fonction correspondante.

[majak]

je ne comprends pas trop la question, je pense qu'il s'agit plutôt d'évaluer :
1+ ln[a(n)]- ln[b(n)].

on a a(n)= n!
et b(n) est donné par la boucle définissant
la fonction suivante :

fb(n):=
{local u,j;
u:=1;
pour j de 1 jusque n
faire u:=j*u + 1;
fpour;
}

Pour moi mon problème est résolu ,
si on demande

p=ln(u(n))/ln(10) +1 :
evalf( u(n) , p) ;

donne un résultat cohérent !

[parisse]

bonjour,

je ne peux pas regarder plus en détails la question avant mercredi, mais en général il faut savoir en ce qui concerne le calcul numérique que les problèmes de précision se produisent si on essaie de soustraire deux nombres très proches (ou additionner deux nombres presque opposés bien sur). Le fait de travailler avec xcas ou avec des nombres flottants "habituels" (double du langage C par exemple) ne change rien (sauf pour la précision qu'on peut demander plus grande dans xcas ou tout logiciel utilisant des flottants multi-précision), et c'est le même problème sur une calculatrice (il peut y avoir des différences suite à la base utilisée, par exemple base 10 sur certaines calculatrices ou sur des logiciels de calcul formel comme maple au lieu de la base 2). Donc si vous arrivez à reformuler le problème de convergence en terme d'un calcul ne faisant pas apparaitre de différence de nombres proches, vous n'aurez plus de "mauvaises surprises numériques".

[alb]

il s'agissait bien sûr de 1+evalf(ln(a(n)))-evalf(b(n))
Je propose sur ce principe la fonction suivante, mais c'est sans garantie:

Code : Tout sélectionner

Suite(u1,n,d):={
  //par ex u1=e-2
  local T,L,j,N,D,F;
  T:=a+u1-e;N:=ln(2);D:=1;F:=0;
  L:=[1,simplify(T),D-N];
  pour j de 1 jusque n faire
    T:=(j+1)*T-1;
    F:=solve(T,a)[0];
    N:=evalf(ln(numer(F)),d);
    D:=1+evalf(ln(denom(F)),d);
    L:=L,[j+1,simplify(T),D-N]; 
  fpour;
  return [L];
}
:;

Il faudrait voir si Suite(e-2,50,50) ou Suite(e,50,50) est cohérent avec les résultats mathématiques.
J'ai également deux questions plus générales:
1)faut-il déclarer a en local ?
2)peut-on forcer solve(6a-16,a)[0] à renvoyer 16/6, cad à ne pas simplifier ?

[parisse]

Code : Tout sélectionner

Suite(u1,n,d):={
  //par ex u1=e-2
  local T,L,j,N,D,F;
  T:=a+u1-e;N:=ln(2);D:=1;F:=0;
  L:=[1,simplify(T),D-N];
  pour j de 1 jusque n faire
    T:=(j+1)*T-1;
    F:=solve(T,a)[0];
    N:=evalf(ln(numer(F)),d);
    D:=1+evalf(ln(denom(F)),d);
    L:=L,[j+1,simplify(T),D-N]; 
  fpour;
  return [L];
}
:;

Il faudrait voir si Suite(e-2,50,50) ou Suite(e,50,50) est cohérent avec les résultats mathématiques.
J'ai également deux questions plus générales:
1)faut-il déclarer a en local ?
2)peut-on forcer solve(6a-16,a)[0] à renvoyer 16/6, cad à ne pas simplifier ?

Sur les questions générales:
On peut déclarer a en local, mais il faut alors purger a (ou utiliser l'instruction assume(a,symbolic), sinon certaines instructions de calcul formel ne vont pas fonctionner car quand on declare une variable locale elle est initialisée à 0. Si a n'est pas déclarée en local, un warning apparait, il faut alors prendre garde que a n'ait pas une valeur globale.
On ne peut pas empêcher la simplification de fractions d'entiers (ou d'entiers de Gauss).

[alb]

oubliez mes remarques précédentes, le comportement de lna-lnb ne donne aucun renseignement sur a-b