[2.0] Erreur de calcul numérique de .net o_O

Version imprimable

Bonjour les gens :)
Je suis en prépa et dans le cadre de mon TIPE, je bosse sur la simulation numérique des équations de propagation d'ondes.
J'ai donc fait mon appli, et ca marche bien, sauf que ca a un comportement bizarre. Je suis donc remonté à la racine du problème et ai créé une appli console simulant le problème en une dimension.
L'équation différentielle discrétisée qui régit l'évolution du problème est :
t3[i] = 2 * t2[i] - t1[i] + pasTemps² * v² / pas²(t2[i + 1] + t2[i - 1] - 2 * t2[i])
Où t3 représente le tableau à l'instant t+dt, t2 celui à l'instant t et t1 celui à l'instant t-dt.
Le point intéressant de cette formule est qu'elle est symétrique pour i : il n'y a aucune raison pour que l'onde parte plutôt à droite ou à gauche. C'est pourtant ce qui se produit.

Voici mon code :
Code:

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 //les constantes, simplifiées au max pour pouvoir faire les calculs à la main const int nbre = 5; //sachant que deux n'interviendront pas dans les calculs mais seront toujours nuls const double pasTemps = 0.1; const double v = 1; const double pas = 1; double[] t1 = new double[nbre]; double[] t2 = new double[nbre]; double[] t3 = new double[nbre]; t1[nbre/2] = 1; t2[nbre/2] = 1;//initialisation for (int j = 0; j < 5; j++) { Console.WriteLine("Etape " +j); for (int i = 1; i <nbre-1; i++) //on laisse tomber les extrémités, qui vaudront toujours 0 { t3[i] = 2 * t2[i] - t1[i] + pasTemps * pasTemps * v * v / pas / pas * (t2[i + 1] + t2[i - 1] - 2 * t2[i]); Console.WriteLine(t3[i]); } t1 = t2; t2 = t3;//on échange et on recommence }
Voici les résultats :
Etape 0

0
0
0,01
0,98
0,01
0
0

Jusque là c'est normal
Etape 1

0
0,0001
0,029601
0,94079601
0,0292079601
0,000292079601
2,92079601E-06
Là par contre c'est pas normal du tout. D'une part c'est dissymétrique, d'autre part le calcul à la main ne donne pas les mêmes résultats (pour la partie en bas). Et en plus, la dernière valeur n'a aucune raison de ne pas être nulle, vu que ses voisines aux deux rangs précédents sont nulles
Le résultat attendu serait :

0
0,0001
0,029601
0,94079601
0,029601
0,0001
0

Je m'interroge donc, est ce que les problèmes viennent du calcul numérique fait par le runtime ? Pourquoi favoriser ainsi une direction ?

Si quelqu'un avait une idée du pourquoi du comment du problème, ca m'arrangerait bien, pasque ca fait deux soirs que je m'arrache les cheveux dessus :)
Smeuuh

11/11/2006, 00h36
Thomas Lebrun

Hum, curieux car en testant ton programme, j'obtiens les bonnes valeurs (si l'on ne prend que 2 chiffres après la virgule):

Citation:

Etape 0
0,01
0,98
0,01

Etape 1
0,0296
0,940796
0,02920796

Etape 2
0,03841596
0,9226563192
0,037850363992

Etape 3
0,046874203992
0,90505043849584
0,0461438610971184

Etape 4
0,0549872242971184
0,887960740579865
0,0541005912809747

Bonjour,

En fait, le vrai problème vient de l'affectation de tes tableaux.

Code:

t1 = t2; t2 = t3;//on échange et on recommence

Tu copies seulement les références
exemple :
t1 = t2 veut dire que t1 et t2 pointe sur la même référence d'objet
Les modifications que tu apporteras au tableau t2 seront répercuté au tableau t1 et réciproquement
Si tu veux copier le tableau il faut passer par Array.copy(...)

Autre chose
Les calculs par défaut de nombres décimaux ne sont pas précis.
Fait un test
0.1*0.1 ne donnera pas 0.01 mais 0.01000000000002

Si tu veux vraiment de la précision, il faut utiliser "decimal"
Voila ce qu'il faut faire :
J'ai décomposé le calcul pour une meilleur lisibilité

Code:

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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
 
//les constantes, simplifiées au max pour pouvoir faire les calculs à la main
const int nbre = 5; //sachant que deux n'interviendront pas dans les calculs mais seront toujours nuls
const decimal pasTemps = 0.1M;
const decimal v = 1;
const decimal pas = 1;
decimal[] t1 = new decimal[nbre];
decimal[] t2 = new decimal[nbre];
decimal[] t3 = new decimal[nbre];
// t3[i] = 2 * t2[i] - t1[i] + pasTemps² * v² / pas² * (t2[i + 1] + t2[i - 1] - 2 * t2[i])
decimal part1; // 2 * t2[i]
decimal part2; // pasTemps²
decimal part3; // v²
decimal part4; // pas²
decimal part5; // part2 * part3 (pasTemps² * v²)
decimal part6; // part5 / part4 (pasTemps² * v² / pas²)
decimal part7; // t2[i + 1] + t2[i - 1]
decimal part8; // part7 - part1 (t2[i + 1] + t2[i - 1] - 2 * t2[i])
decimal part9; // part6 * part8 ' pasTemps² * v² / pas² * (t2[i + 1] + t2[i - 1] - 2 * t2[i])
decimal part10; // part1 - t2[i] (2 * t2[i] - t1[i])
// t3[i] = part10 + part9
t1[nbre/2] = 1;
t2[nbre/2] = 1;//initialisation
for (int j = 0; j < 5; j++)
{
    Console.WriteLine("Etape " +j);
    for (int i = 1; i <nbre-1; i++) //on laisse tomber les extrémités, qui vaudront toujours 0
    {
        part1 = decimal.Multiply(2, t2[i]);
        part2 = decimal.Multiply(pasTemps, pasTemps);
        part3 = decimal.Multiply(v, v);
        part4 = decimal.Multiply(pas, pas);
        part5 = decimal.Multiply(part2, part3);
        part6 = decimal.Divide(part5, part4);
        part7 = decimal.Add(t2[i + 1], t2[i - 1]);
        part8 = decimal.Subtract(part7, part1);
        part9 = decimal.Multiply(part6, part8);
        part10 = decimal.Subtract(part1, t1[i]);
        t3[i] = decimal.Add(part10, part9);
        Console.WriteLine(t3[i]);
    }
    // Ici, il y a une erreur
    // tu passes les références
    // t1 = t2 => veut dire que t1 et t2 ont le même pointeur
    // donc les modifications de valeurs sur t2 seront affecté à t1 et reciproquement...
    //t1 = t2; t2 = t3;//on échange et on recommence
    // Tu peux faire comme ceci
    Array.Copy(t2,0,t1,0,nbre);
    Array.Copy(t3,0,t2,0,nbre);
    // ou comme ceci
    // for (int i = 0; i <nbre; i++) 
    // {
    //     t1(i) = t2(i);
    //     t2(i) = t3(i);
    // }
}
Console.Read();

P.S. : ton erreur au moment du transfert des tableaux est expliqué dans le code

11/11/2006, 12h27
Thomas Lebrun

Je suis tout de même curieux de voir que chez moi, cela fonctionne.... :?
11/11/2006, 13h51
Smeuuh

1 pièce(s) jointe(s)

Merci pour vos réponses :)

Merci Kelpan !
C'était en effet le problème. J'ai remplacé par une copie, tout marche impec :)
Les résultats avec l'erreur étaient étranges (propagation d'un seul côté, en hyperbole, l'image est attachée au post si vous voulez voir).
Concernant l'utilisation des décimaux, je pense que la précision actuelle me suffit. Mon application en 2D avec visualisation rend maintenant de très bons résultats :)
Par contre, y a-til une raison particulière pour utiliser decimal.Multiply(a,b) plutôt que a*b(voire (decimal)a*b si a est un double), et pour décomposer le calcul, mis à part la lisibilité ?

Thomas, si cela marche chez toi c'est que je me suis trompé, le code proposé n'est pas celui qui donne ces résultats. Il faut remplacer nbre = 5 par nbre = 7, et tu trouveras les mêmes résultats.

Meric encore :)
11/11/2006, 15h26
Kelpan
Citation:

Envoyé par Smeuuh

Par contre, y a-til une raison particulière pour utiliser decimal.Multiply(a,b) plutôt que a*b

Je l'ai dis dans mon post précédente, voici la différence :
Code:

1 2 3 0.1 * 0.1 = 0.01000000000000002 decimal.Multiply(0.1,0.1) = 0.01
La décomposition du calcul n'est pas une obligation.
Mais si tu ne le fais pas tu risque de te retrouver avec des imbrications qui risque d'être difficile à lire :

Code:

T3[i] = decimal.Multiply(decimal.Add(x,y),decimal.Substrac(z,decimal.Multiply(b,r)),decimal.Add(....)))
11/11/2006, 15h41
Smeuuh

Oui, si tu parles de multiplication de doubles.
Mais a*b (où a et b sont des décimaux) devrait n'être qu'un raccourci pour decimal.Multiply(a,b) non ?
De plus, le code :
double a = 0.1; double b = 0.1;
Console.WriteLine(a*b);
rend chez moi un beau 0.01 tout propre

J'ia mal compris quelque chose ?
11/11/2006, 16h24
Kelpan
Oui, désolé mon exemple n'est pas si visuel que ça.

Si tu te mets en mode déboggage, et que tu met le calcul en espion tu réaliseras le problème.

Sinon j'ai un exemple plus visuel :
Code:

1 2 3 Console.WriteLine(0.1 * 0.1 - 0.01); Console.WriteLine(decimal.Subtract(decimal.Multiply(0.1M,0.1M),0.01M));
11/11/2006, 17h04
Smeuuh

Effectivement :)

Cependant Console.WriteLine(0.1M * 0.1M - 0.01M); produit 0, donc c'est pas la peine de s'emmerder avec les decimal.Multiply du moment qu'on fait les casts qui vont bien, non ?

En tout cas merci, j'ignorais totalement tout ca ^^