comment faire un développement limité qui converge

**Rniamo** · 16/10/2006, 22h37

Bonjour à tous,
en cours j'ai abordé le C et en TD on a du faire la fonction cosinus via un dl (développemnt limité).
J'ai voulu tester ma fonction sur un PC (c'est plus parlant quand même) et quelle ne fut pas ma surprise quand pour cos(33) j'ai obtenu plus de 3000!!!...il me semblait que |cos(x)|<=1!

Donc voici ma question : comment corriger ceci (une bonne idée serait la bienvenue).

Je joins mon code pour la fonction cos (basée sur le fait que cos(x)=somme ((-1^n)*(x^2n)/factoriel(2n)) :

double cosinus(double x)
{
double resultat=0;double signe=1;double terme=1;long n=0;
while (resultat+terme!=resultat)
{
resultat=resultat+signe*terme;
signe=-signe;n=n+1;
terme=terme*((x*x)/(2*n*(2*n-1)));
}
return resultat;
}

ma fonction est bonne certes mais [-24;24] est le plus grand intervalle où je peux avoir 6 décimales éxactes...un peu juste pour un supercalculateur...

**Médinoc** · 16/10/2006, 22h39

Tu fais une comparaison d'égalité sur des nombres à virgule flottante, c'est très mauvais.
Tu dois tester si leur différence est inférieure à un certain seuil de tolérance...

**Emmanuel Delahaye** · 16/10/2006, 22h43

Envoyé par Médinoc

Tu fais une comparaison d'égalité sur des nombres à virgule flottante, c'est très mauvais.
Tu dois tester si leur différence est inférieure à un certain seuil de tolérance...

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
#include <float.h>
 
   if (fabs(a - b) < DBL_EPSILON)
   {
      ...

mat.M · 16/10/2006, 22h49

Aaarghh attention en C/C++ les angles sont en radians
Il faut les convertir donc en divisant par 180 et multiplier par PI ( ou vice versa )

Ce qui ne semble pas le cas dans le code donné

Tiens petite fonction de rotation ( Merci Mr Petzold )

#define TWOPI ( 2* 3.14159)

void RotatePoints(POINT pt[], int iNum, int iAngle)
{
int i ;
POINT ptTemp ;

for (i = 0 ; i < iNum ; i++)
{
ptTemp.x = (int) (pt[i].x * cos (TWOPI * iAngle / 360) +
pt[i].y * sin (TWOPI * iAngle / 360)) ;

ptTemp.y = (int) (pt[i].y * cos (TWOPI * iAngle / 360) -
pt[i].x * sin (TWOPI * iAngle / 360)) ;

pt[i] = ptTemp ;
}
}

**Rniamo** · 17/10/2006, 00h06

excusez moi mais

-> vous vous trompez, je n'ai pas besoin de multiplier par 180 ; la fonction marche mais pour des petites valeurs (donc le 180 ne sert à rien...et surtout un dl est défini en radian...question d'homogéinité).

-> la question n'était pas sur le test (ok il est nul mais je vais le corriger).

Je repose donc ma question : pourquoi ma fonction dépasse le 1 majorant de la fonction
x -> cos(x)
?

**Jean-Marc.Bourguet** · 17/10/2006, 09h04

Envoyé par Médinoc

Tu fais une comparaison d'égalité sur des nombres à virgule flottante, c'est très mauvais.
Tu dois tester si leur différence est inférieure à un certain seuil de tolérance...

Les séries convergentes sont un des cas où un test d'égalité a du sens.

**Jean-Marc.Bourguet** · 17/10/2006, 09h07

Envoyé par Emmanuel Delahaye

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
#include <float.h>
 
   if (fabs(a - b) < DBL_EPSILON)
   {
      ...

Même si on n'était pas dans un cas où l'égalité à du sens, ce test est beaucoup plus mauvais qu'un test d'égalité. Il revient à un test d'égalité pour des valeurs absolues plus grandes que 2, il est complètement absurde pour des valeurs inférieures à DBL_EPSILON.

**Emmanuel Delahaye** · 17/10/2006, 09h14

Envoyé par Jean-Marc.Bourguet

Même si on n'était pas dans un cas où l'égalité à du sens, ce test est beaucoup plus mauvais qu'un test d'égalité. Il revient à un test d'égalité pour des valeurs absolues plus grandes que 2, il est complètement absurde pour des valeurs inférieures à DBL_EPSILON.

C'est ma fête en ce moment...

J'ai pas compris pourquoi. Les paramètres et retour de fabs() ne sont pas de type double ? (j'aurais peut être dû vérifier, je me suis fié au 'f' comme 'double' !)

Ca me parait correct :

http://man.developpez.com/man3/fabs.3.php
EDIT :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 
#include <stdio.h>
#include <math.h>
#include <float.h>
 
int is_equalf_ed (double a, double b)
{
   return fabs (a - b) < DBL_EPSILON;
}
 
int main (void)
{
   double n;
   int i = 0;
 
#define MAX 0.001
#define MIN -MAX
 
   for (n = MIN; n < MAX; n += 0.00001)
   {
      if (n == 0)
      {
         printf ("%f (%e) (n==0)\n", n, n);
      }
 
      if (is_equalf_ed (n, 0))
      {
         printf ("%f (%e) (is_nulf_ed)\n", n, n);
      }
 
      i++;
   }
   printf ("%d values tested from %e to %e\n", i, MIN, MAX);
   return 0;
}

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
0.000000 (1.927847e-018) (is_equalf_ed)
200 values tested from -1.000000e-003 to 1.000000e-003
 
Press ENTER to continue.

**Jean-Marc.Bourguet** · 17/10/2006, 09h14

Le problème est que les doubles ont une précision relative fixe. Pour que la série que tu emploies converge avec 6 chiffres significatifs, il faut que tous les calculs soient fait avec une précision absolue supérieure à ces 6 chiffres. Or, plus ton argument est grand, plus les résultats intermédiaires sont élevés. En calculant cos(24) tu as des résultats intermédiaires de l'ordre de 10^9, ce qui te permet d'avoir encore 6 chiffres significatifs sur le résultat; avec des nombres plus grands, tes résultats intermédiaires sont plus élevés également et tu commences à perdre en précision.

Pour profiter au mieux de la précision des doubles, il faut d'abord réduire ton argument pour faire un calcul d'un sinus entre -PI/2 et PI/2.

**Jean-Marc.Bourguet** · 17/10/2006, 09h48

Envoyé par Emmanuel Delahaye

J'ai pas compris pourquoi

Le problème sur les flottants, c'est que je suis à la limite de mes compétances. Il faudrait un Gaby ou un Laurent Deniau pour vérifier que je ne dise pas de grosses bétises.

DBL_EPSILON, c'est le plus petit nombre tel que 1+DBL_EPSILON != 1, ce qui veut dire -- en supposant des flottants binaires -- que 2+DBL_EPSILON == 2.
Donc pour a plus grand que deux fabs(a-b) < DBL_EPSILON, c'est la même chose que fabs(a-b) == 0.

Si on part de l'autre coté, dans les nombres beaucoup plus petit que 1, faire fabs(a-b) < DBL_EPSILON, c'est prendre une marge relativement très grande. Donc il faut faire quelque chose comme:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

fabs(a-b) <= f * (fabs(a) + fabs(b))

pour avoir une borne qui soit relative à la taille des nombres comparés.

Quelle valeur prendre pour f? La ce n'est pas possible de se déterminer sans connaître le problème exact. Dans le cas présent d'une série convergente, f == 0 me ne semble pas être mauvais, s'il doit être non nul, le _EPSILON pour le type est un strict minimum et il me semble douteux que ce soit la bonne valeur: elle ne permet une différence que sur le bit le plus faible.

**Rniamo** · 17/10/2006, 16h54

Merci Jean-Marc.Bourguet mais j'ai l'impression que vous vous focalisez sur la condition d'arrêt du terme (cela marche sans problème...surtout pour les petites valeurs)...maintenant que j'y pense j'ai peut-être mal exprimé mon problème.

Mon problème est donc de faire un dl qui soit pas trop mauvais.
Par exemple la fonction cosinus du C (de math.h) va facilement jusqu'à 1 million (j'ai pas testé au-delà)...la mienne s'arrête à 24...

En ce qui concerne le "décalage de x" (c'esy à dire que si on a à calculer cos(2*pi+1), on calcule cos(1), ceci nécessite une valeur éxacte de pi (difficile en informatique) ; j'ai testé avec une valeur à 20 décimales...résultat peu concluant.

Voice, j'espère avoir était plu clair...

**millie** · 17/10/2006, 17h01

J'espère que ce que je vais dire n'a pas été dit.

Ca dépasse 1, moi je dis, pourquoi ça ne dépasserait pas, tu approches cos par un polynome qui diverge. Si tu utilises un polynome de n eme degré, alors la méthode que tu utilises approcheras assez bien sur un intervalle, mais comme ça diverge, ça n'ira jamais.

Cela dit, pourquoi ne pas ramener toutes tes valeurs entre - pi et pi (vu que cos(2pi +x ) = cos(x)). Ainsi, tu as besoin d'un développement à un ordre beaucoup moins important pour avoir une très bonne approximation.

ceci nécessite une valeur éxacte de pi (difficile en informatique)

Ca sera beaucoup plus concluant que de ne pas se ramener dans l'intervalle.

edit Ca aurait été un post à faire plutôt dans le forum Algorithmique.

**Rniamo** · 17/10/2006, 17h07

euh...millie, pour rapporter x entre -pi et pi c'est raté parce qu'il faut pi...

et puis en math on prouve facilement que la somme (infinie) des (-1^n)*(x^2n)/(2n)! converge vers cos...donc vu que mon "infini" est suffisemment grand (j'ai testé les puissances de 10 entre 100 et 1000000) les terme suivant (en 1/1000...) sont négligeables...

donc ceci ne m'aide pas vraiment...mais merci quand même

**Jean-Marc.Bourguet** · 17/10/2006, 17h08

Envoyé par Rniamo

Merci Jean-Marc.Bourguet mais j'ai l'impression que vous vous focalisez sur la condition d'arrêt du terme (cela marche sans problème...surtout pour les petites valeurs)...maintenant que j'y pense j'ai peut-être mal exprimé mon problème.

Je ne me focalise pas dessus, j'ai écrit que le test d'équalité n'était pas mauvais dans le cas présent et j'ai essayé d'expliquer quelques erreurs dans les alternatives proposées.

Pour me répéter: pour que ton calcul soit précis, il faut travailler avec une précision absolue constante, c'est à dire toujours le même nombre de chiffres après la virgule. Or les flottants ont une précision relative constante: ils ont toujours le même nombre de chiffres significatifs. Plus tes arguments sont élevés, plus ils ont besoin de chiffres avant la virgule, perdant ceux après.

En ce qui concerne le "décalage de x" (c'esy à dire que si on a à calculer cos(2*pi+1), on calcule cos(1), ceci nécessite une valeur éxacte de pi (difficile en informatique) ;

C'est pourtant comme ça qu'il faut faire. Te rends-tu compte à quel point ta série converge lentement pour des arguments plus grand que 1?

j'ai testé avec une valeur à 20 décimales...résultat peu concluant.

J'en doute. Pour arriver à tester avec une valeur qui a plus de 15 décimales, tu aurais dû employer des techniques que visiblement tu ne connais pas. Si tu les comprenais, tu ne poserais pas les questions que tu poses.

**millie** · 17/10/2006, 17h08

Tu peux donner une approximation de pi à 10^-20 près (ou plus). Je ne vois pas le problème...

Ta série converge, mais il faudrait au moins que tu saches elle converge comment. Indication : ça converge vraiment lentement pour les grandes valeurs de x.

donc vu que mon "infini" est suffisemment grand

Ton infini ne sera justement jamais assez grand car : suffisament grand différent de infini.

Pour tout n, Somme((-1^k)*(x^2k)/(2k)!, k = 0 à n) est divergente...

**Rniamo** · 17/10/2006, 17h09

...j'ai oublié de préciser que en me ramenant dans [-pi;pi] l'erreur restée monumentale (sur le 3eme chiffre à peu prés pour des nombre de l'ordre de 400).

**Elijha** · 17/10/2006, 17h09

terme=terme*((x*x)/(2*n*(2*n-1)));

Attention à l'écriture, terme_1 est différent de terme_2

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
terme_1 = terme_1*((x*x)/(2*n*(2*(n-1)))) ;

terme_2 = terme_2*((x*x)/(2*n*((2*n)-1))) ;

Faire très attention avec les paranthèses !

**Jean-Marc.Bourguet** · 17/10/2006, 17h13

Envoyé par millie

Tu peux donner une approximation de pi à 10^-20 près (ou plus). Je ne vois pas le problème...

La mettre dans un double qui a 15 chiffres significatifs.

**Rniamo** · 17/10/2006, 17h14

je sais Elijha, j'y fait attention mais en C la multiplication est prioritaire sur la'addition donc c'est bon ... je veux bien faire (2n-2)!x^2/(2n)!=x*x/((2*n-1)*2*n)