Calcul approchés ? evalf se trompe ?

[yvand]

Bonjour,

Pour un partiel j'ai eu recours à Xcas mais je ne comprend pourquoi le résultat qu'il donne est aberrant alors qu'il est censé travailler en exact.

L'objectif était de conjecturer la limite d'une suite d'intégrales :

Code : Tout sélectionner

I(n) := int(x^n*e^(x^2),x,0,1)

On peut facilement montrer que :

• I(n) vérifie la relation de récurrence I(n+2) = 1/2*exp(1) - (n+1)/2*I(n)
• I(n)est décroissante et minorée (car positive) donc convergente et que sa limite est nulle (à partir de la relation de récurrence par exemple).

Voici l'algorithme utilisé :

Code : Tout sélectionner

n := 1;
u := 1/2*exp(1)-1/2;
tantque n < 45 faire
  u := 1/2*exp(1)-(n+1)/2*u;
  n := n+2;
  afficher(n+": "+evalf(u));
ftantque;

J'affiche les valeurs de I(n) pour n=1,3,5,...,45.
Xcas est censé travailler en exact car je n'ai pas mis de virgules (1/2 et non 1./2 ou 0.5).
C'est lors de l'evalf que ça fourre. Le problème est que au-delà de 35, les résultats sont aberrants : la suite n'est pas décroissante et certains termes sont négatifs (I(39) par exemple). Or ceci est impossible en théorie !

Comment expliquer cela ?

[parisse]

C'est juste que evalf sans 2eme argument evalue l'expression avec des flottants double precision, ce qui ne suffit pas pour compenser les pertes de precision. En mettant par exemple evalf(.,50) au lieu de evalf on a bien une suite de valeurs numeriques decroissantes.
En effet l'erreur absolue sur I(n) est multipliee par (n+1)/2 losrqu'on applique la recurrence, et comme la valeur de I(n) est decroissante, l'erreur relative augmente donc par un facteur au moins (n+1)/2, a un moment on depasse la precision (10^(-15) environ pour des doubles), donc l'evaluation de l'expression (qui obeit a la meme regle que la recurrence comme on le voit en tapant u sur un niveau tout seul) donne n'importe quoi.
Donc ici il vaut mieux faire une evaluation numerique de l'integrale (en mettant une borne 1.0 par exemple)

[yvand]

Merci pour votre réponse rapide et les explications.

Quand on fait (calcul de I(43) avec borne sup de l'intégrale = 1) :

Code : Tout sélectionner

int(x^43*e^(x^2),x,0,1)

On obtient bien la valeur exacte de I(43).
Pour obtenir une valeur approchée en faisant evalf (sans 2e argument), Xcas se trompe car les deux nombres (affichés pour la valeur exacte : -9397653627525472270*exp(1)+25545471085854720000) sont très proches et que exp(1) est remplacé par une valeur approchée. C'est cela ?

[alb]

mais alors le warning n'est pas tres clair:
int(x^91*e^(x^2),x,0,1)
Erreur, la valeur numerique ne correspond pas a la valeur exacte de l'integrale, on renvoie les deux!
[-22003277881839850715808338506873781370572398889102040802*exp(1)+59811110432740097280981580747828857532191866880000000000,0.0289307266231]
Comment doit-on comprendre "Erreur" ? La valeur exacte est-elle fausse ? Que doit faire l'utilisateur ?
[Edit] je n'avais pas vu le post precedent

[parisse]

Oui, en fait ça devrait être qualifié de warning plutot que d'erreur: l'intégration formelle avec bornes effectue une intégration numérique de vérification et compare avec evalf sur la forme exacte, ici c'est un problème d'évaluation numérique de la formule exacte et pas une erreur de la valeur exacte. Je pourrais peut-etre réessayer de faire un evalf avec plus de décimales si les 2 valeurs ne correspondent pas pour éviter ça (enfin dans une certaine mesure).

[alb]

Quel test declenche-t-il le message d'erreur ?

[parisse]

C'est le test d'egalite numerique entre evalf(la valeur exacte) et le calcul de l'integrale par quadratures. J'affine le test en evaluant numeriquement avec 256 bits de mantisse, ca repoussera le warning.