Ist wirklich der round-off Error so gross????

Davood

Hallo C++ Experten!

in meinem Programm berechne ich eine Summierung und diese Summirung muss null werden. Aber der Code in C++ gibt einen sehr groben Wert (grob für mich) und zwar 6.50521E-19. Ich habe den Code mit Visaul C++ kompiliert. Das kann aber nicht sein so ein round-off error. Das probelm hatte ich nie in Fortran. Hier unten ist der sehr einfache Code, der die selbe Berechnungen macht: (wenn man diese Berechungen z.B. in Excell ausführt, kommt genau 0.0 als Endwert raus!!!!) Ist es ein algemeines Problem oder muss ich in Compiler eine Option einschalten??? Ich hoffe, dass jemand mir helfen kann!

#include <stdio.h>
#include <iostream.h>
#include <string.h>
#include <fstream.h>

double gravty = 9.80933;
double x1,x2,x3;
double y1,y2,y3;
double z1,z2,z3;

double h1,h2,h3;

double wse, dt, A, slopX;
double bedSlope1, bedSlope2, bedSlope3;
double pressurX1, pressurX2, pressurX3;
double ResVx = 0.0;

void main()
{
   x1 = 6.25e+1, x2 = 6.25e+1, x3 = 6.49999976e+1;
   y1 = 1.40000001e+1, y2 = 1.20000005e+1, y3 = 1.20000005e+1;
   z1 = 0.0, z2 = 0.0, z3 = 5.0e-1;

   wse = 0.8;
   dt = 0.001;

   h1 = wse - z1;
   h2 = wse - z2;
   h3 = wse - z3;

   A = 0.5 * ((x2*y3-x3*y2)+(x3*y1-x1*y3)+(x1*y2-x2*y1));
   slopX = 0.5 / A *(z1*(y2-y3)+z2*(y3-y1)+z3*(y1-y2));

   bedSlope1 = gravty * 0.5 * (h1 + h2) / 3.0;
   bedSlope2 = gravty * 0.5 * (h2 + h3) / 3.0;
   bedSlope3 = gravty * 0.5 * (h3 + h1) / 3.0;

   pressurX1 = 0.5 * gravty * 1./3. * (h1 * h1 + h2 * h2 + h1 *h2) *(y1 - y2);
   pressurX2 = 0.5 * gravty * 1./3. * (h2 * h2 + h3 * h3 + h2 *h3) *(y3 - y2);
   pressurX3 = 0.5 * gravty * 1./3. * (h3 * h3 + h1 * h1 + h3 *h1) *(y1 - y3);

   printf(" h1 = %15.12f , h2 = %15.12f , h3 = %15.12f \n", h1,h2,h3);
   printf(" A = %15.12f , slopX = %15.12f \n", A,slopX);
   printf(" bedSlope1 = %15.12f ,  bedSlope2 = %15.12f ,  bedSlope3 = %15.12f \n",  bedSlope1, bedSlope2, bedSlope3);
   printf(" pressurX1 = %15.12f ,  pressurX2 = %15.12f ,  pressurX3 = %15.12f \n",  pressurX1, pressurX2, pressurX3);

   ResVx -= dt * slopX * bedSlope1;
   ResVx += dt * pressurX1 / A;
   printf(" ResVx = %15.12e  \n", ResVx);   

   ResVx -= dt * slopX * bedSlope2;
   ResVx -= dt * pressurX2 / A;
   printf(" ResVx = %15.12e  \n", ResVx);   

   ResVx -= dt * slopX * bedSlope3;
   ResVx -= dt * pressurX3 / A;
   printf(" ResVx = %15.12e \n", ResVx);   
}

Bitte Code-Tags benutzen!

[ Dieser Beitrag wurde am 18.06.2003 um 13:47 Uhr von HumeSikkins editiert. ]

PeterTheMaster

krass. 1:1 fortran.
es gibt keine allgemeine moeglichkeit, den rundungsfehler loszuwerden. hoechstens eine bibliothek die symbolisch rechnet.
du koenntest im einzelfall versuchen, die reihenfolge der operationen vorteilhaft zu veraendern. als /3. am schluss, weil das mit sicherheit ungenau ist oder statt (gross+klein)-gross lieber (gross-gross)+klein oder solche sachen halt.
wieviel bit hatte denn dein fortran double?

Davood

Hallo Peter

danke für den Tipp. Ich werde ihn probieren. Aber was mich wundert ist die Grosse des Werts e-019. Wenn er z.B. e-100 waere konnte ich gut damit leben.
Laut Compaq Visaul Fortran Double nimmt 8 bytes von memory.

danke noch mal

PeterTheMaster

dann wundert es mich eigentlich, dass beim fortran double was anderes als hier rauskommt.
warum solle der exponent so hoch sein? wenn die groesste belegte stelle im einer oder zehntel bereich ist, hat er doch gar keine chance, so kleine fehler zu machen.
edit: vielleicht kann ja der fortran compiler die ausdruecke entsprechend optimieren, vergleich doch mal den asm output.

[ Dieser Beitrag wurde am 18.06.2003 um 20:13 Uhr von PeterTheMaster editiert. ]

Bashar

wenn in Fortran 0.0 rauskommt, kann das nicht vielleicht unter Umständen daran liegen, dass die Ausgabefunktion bereits rundet? Bei Excel kann ich mir das auf jeden Fall vorstellen.

Helium

Also am wichtigesten ist so wenig rechnen wie möglich. also aus

0.5 * gravty * 1./3. * ...

ein 1./6. * gravty * ... machen, etc.

Dann gibts noch verschiedene Arten, wie der Prozessor rundet (zumindest bei Intel-kompatiblen). Das kan man IMHO nur per ASM umstellen. Vielleicht ist das, was der MSVC verwendet, schlecht geeignet für deine rechnungen.

was ist an blah*10^-19 so gross?

ich weiss nicht ob das was hilft:
wenn ich die Variablen als long double deklariere dann erhalte ich haargenau 0.0000(..

float vergrössert den fehler auf 1.7e-10

OOOOOPS....

ich bin vorher einem Fehler aufgesessen. (Die printf funktion verlangt %lf für long double)

Der Fehler bei long double beträgt noch 2.1684e-19

Davood

Liebe Freunde,

danke für die Tipps. Was mich noch stört ist das so.
Ich reduziere die oprationen in folgendweise (die Konstante werde nich mehr in den Oprationen eingeschrieben!!!) Damit erhoe ich die Grosse der Werte, trotzdem ist der Fehler noch da. Der Wert von sum sollte eigentlich gleich wie Endergebniss von ResVx sein!!!

double gravty = 1.0;
double x1,x2,x3;
double y1,y2,y3;
double z1,z2,z3;

double h1,h2,h3;

double wse, dt, A, slopX;
double bedSlope1, bedSlope2, bedSlope3;
double pressurX1, pressurX2, pressurX3;
double ResVx = 0.0;

void main()
{
x1 = 6.25e+1, x2 = 6.25e+1, x3 = 6.49999976e+1;
y1 = 1.40000001e+1, y2 = 1.20000005e+1, y3 = 1.20000005e+1;
z1 = 0.0, z2 = 0.0, z3 = 0.5; // z3 = 5.0e-1

wse = 0.55;
cout << " Enter the Water Surface Elevation? wse > 0.55";
cin >> wse;
dt = 1.0;

h1 = wse - z1;
h2 = wse - z2;
h3 = wse - z3;

A = 0.5 * ((x2*y3-x3*y2)+(x3*y1-x1*y3)+(x1*y2-x2*y1));
slopX = (z1*(y2-y3)+z2*(y3-y1)+z3*(y1-y2));

double tmp = 1./3.;
bedSlope1 = h1;
bedSlope2 = h2;
bedSlope3 = h3;

pressurX1 = (h1 * h1 + h2 * h2 + h1 *h2) *(y1 - y2);
pressurX2 = (h2 * h2 + h3 * h3 + h2 *h3) *(y3 - y2);
pressurX3 = (h3 * h3 + h1 * h1 + h3 *h1) *(y1 - y3);

printf(" x1 = %19.15f , x2 = %19.15f , x3 = %19.15f \n", x1,x2,x3);
printf(" y1 = %19.15f , y2 = %19.15f , y3 = %19.15f \n", y1,y2,y3);
printf(" z1 = %19.15f , z2 = %19.15f , z3 = %19.15f \n", z1,z2,z3);
printf(" h1 = %19.15f , h2 = %19.15f , h3 = %19.15f \n", h1,h2,h3);
printf(" A = %19.15f , slopX = %19.15f \n", A,slopX);
printf(" bedSlope1 = %19.15f , bedSlope2 = %19.15f , bedSlope3 = %19.15f \n", bedSlope1, bedSlope2, bedSlope3);
printf(" pressurX1 = %19.15f , pressurX2 = %19.15f , pressurX3 = %19.15f \n", pressurX1, pressurX2, pressurX3);

double in1 = (-slopX * (h1+h2+h3));

double in2 = (+pressurX1 - pressurX2 - pressurX3);
// in2 = -(y3*(z2-z1)+y2*(z1-z3)+y1*(z3-z2))(z1+z2+z3-3wse);
double sum = (in1 + in2);

printf(" sum = %19.15e \n", sum);

ResVx -= dt * slopX * bedSlope1;
ResVx += dt * pressurX1;
printf(" ResVx = %19.15e \n", ResVx);

ResVx -= dt * slopX * bedSlope2;
ResVx -= dt * pressurX2;
printf(" ResVx = %19.15e \n", ResVx);

ResVx -= dt * slopX * bedSlope3;
ResVx -= dt * pressurX3;
printf(" ResVx = %19.15e \n", ResVx);
}

Walli

1.) Benutze doch bitte CodeTags
2.) Du solltest dir angewöhnen globale Variablen zu vermeiden indem du lokale benutzt und diese bei Erstellung auch initialisierst. Du programmierst C++ und nicht irgendein uralt-C oder FORTRAN. siehe hier

[ Dieser Beitrag wurde am 19.06.2003 um 18:10 Uhr von MaSTaH editiert. ]

Du sagtest, du hättest VC++. Ich weiß jetzt natürlich nicht welche Version (6 oder 7 (7 == .NET)).

Also bei Version 7 (.NET) klickst du mit Rechts auf dein Projekt und sagst "Eigenschaften". Da dann unter "C/C++"->Optimierung->Gleitkommakonsistenz kann man auswählen "Konsistenz verbessern". Ich habe keine Ahnung was das genau bewirkt, weil ich es bisher nicht benutzt habe, aber vielleicht hilft dir das ja.

Ansonsten sollte man ein paar Regeln bei diesen Berechnungen beachten:
-niemals mit zu kleinen Zahlen rechnen, weil die Genauigkeit beschränkt ist
-niemals mit zu großen Werten rechnen, weil dann ebenfalls die Genauigkeit flöten geht
-niemals ganz oft Werte addieren, wenn es eine Multiplikation auch täte (weil Zahlen nur ungenau repräsentiert werden können und sich der Fehler bei der Addition dann kumuliert, er bei EINER Multiplikation aber nur gering ausfällt)

Daher sollte man immer erst multiplizieren bevor man dividiert. Aber mit ewig großen Zahlen sollte man auch nicht rechnen.
Aber wahrscheinlich weißt du das ja schon alles.

Ach ja, bei Visual C++ 6.0 müsste es auch eine Einstellung geben, die die Genauigkeit erhöht (bzw. die sich so anhört als täte sie das

Jan.