Fehler in Berechnungen

[Christian (Gast)]

Hallo,

ich habe eine Funktion zum Umwandeln von GPS-Koordinaten ins UTM-Format 
geschrieben. Als ich sie zuerst auf dem PC in VB ausprobiert hatte, 
kamen exakte Ergebnisse raus.
Jetzt hab ich sie aber in C für einen ATmega32 geschrieben und der 
Rechtswert lag um ca. 1,5m daneben, was nicht so schlimm wäre. Beim 
Hochwert stimmten aber nur die ersten drei Stellen, was einem Fehler von 
über 1000m entspricht. Ich habe sowohl in VB als auch in C 
Gleitkommazahlen einfacher Genauigkeit verwendet.
Beim Hochwert ist schon die Meridianbogenlänge so ungenau.
Woran kann dieser große Fehler liegen? sind die float-Variablen zu 
ungenau?

Hier der Code:

  const float A  = 6378137.0;    // grosse WGS84-Halbachse
  const float B  = 6356752.3;    // kleine WGS84-Halbachse
  const float K0 = 0.9996;    // UTM Streckfaktor
  const float L0 = 0.1570796;    // UTM-Zentralmeridian fuer Zone 32 (Bogenmass)
  const float e2 = 0.0067395;    // (pow(A, 2) - pow(B, 2)) / pow(A, 2); 
  const float n  = 0.0016792;    // (A - B) / (A + B);

  float Lon = 7 * (M_PI / 180);   // Geographische Laenge als Bogenmass
  float Lat = 49 * (M_PI / 180);  // Geographische Breite als Bogenmass
  float G;            // Meridianbogenlaenge
  float Hw;            // Hochwert
  float Rw;            // Rechtswert

  // Parameter zur Berechnung der UTM-Koordinaten
  float nu2 = (e2 / (1 - e2)) * pow(cos(Lat), 2);
  float t   = tan(Lat);
  float N   = A / (sqrt(1 - e2 * pow(sin(Lat), 2)));
  float l   = Lon - L0;
  
  // Meridianbogenlaenge berechnen
  G = (A / (1 + n)) * ((1 + pow(n, 2) / 4) * Lat - (3 / 2) * (n - pow(n, 3) / 8) * sin(2 * Lat) + (15 / 16) * pow(n, 2) * sin(4 * Lat) - (35 / 48) * pow(n, 3) * sin(6 * Lat));
  
  // Hochwert berechnen
  Hw = K0 * (G + ((N * t) / 2) * pow(l, 2) * pow(cos(Lat), 2) + ((N * t) / 24) * (5 - pow(t, 2) + 9 * nu2) * pow(l, 4) * pow(cos(Lat), 4) + ((N * t) / 720) * (61 - 58 * pow(t, 2) + pow(t, 4) + 270 * nu2 - 330 * pow(t, 2) * nu2) * pow(l, 6) * pow(cos(Lat), 6));
  
  // Rechtswert berechnen
  Rw = 500000 + K0 * (N * l * cos(Lat) + (N / 6) * (1 - pow(t, 2) + nu2) * pow(l, 3) * pow(cos(Lat), 3) + (N / 120) * (5 - 18 * pow(t, 2) + pow(t, 4) + 14 * nu2 - 58 * pow(t, 2) * nu2) * pow(l, 5) * pow(cos(Lat), 5) + (N / 5040) * (61 - 479 * pow(t, 2) + 179 * pow(t, 4)) * pow(l, 7) * pow(cos(Lat), 7));

Vielen Dank für eure Hilfe!

Grüße, CHristian

[Rolf Magnus (Gast)]

> Woran kann dieser große Fehler liegen? sind die float-Variablen zu
> ungenau?

Höchstwahrscheinlich. Kann es sein, daß in deinem VB-Code die 
Berechnungen mit doppelter Genauigkeit durchgeführt werden und damit 
auch die Zwischenwerte genauer sind? Dann wird erst beim Schreiben in 
die Variable die Genauigkeit reduziert. Beim AVR dagegen sind die 
Berechnungen immer mit einfacher Genauigkeit. Sind die Ergebnisse des 
C-Code genauer, wenn du ihn auf dem PC ausführst?

[Karl H. (kbuchegg)]

Christian schrieb:

> Woran kann dieser große Fehler liegen? sind die float-Variablen zu
> ungenau?

Als erstes ersetzt du mal alle pow Aufrufe
Das ist doch Unsinn pow(n,2) über Logarithmen rechnen zu lassen, wenn es 
ein simples n*n auch tut.

>   G = (A / (1 + n)) * ((1 + pow(n, 2) / 4) * Lat - (3 / 2)

Aha.
Dir ist bewusst, dass 3 / 2 ein Ergebnis von 1 liefert und keineswegs 
1.5?

http://www.mikrocontroller.net/articles/FAQ#Datentypen_in_Operationen

Ditto für noch viele andere Divisionen mit konstanten Zahlen, die in 
deinem Code vorkommen.

Bei deinen Berechnungen wirst du allerdings nie ein einigermassen 
genaues Ergebnis erhalten. 6 signifikante Stellen sind nun mal nicht 
viel. Aber ersetz erst mal die offensichtlichen Fehler, dann sollte es 
schon genauer werden.

Und wenn du dann auch noch gleiche Berechnungen vorziehst und in 
Hilfsvariablen ablegst, dann kannst du sogar die Laufzeit noch ein wenig 
drücken. sin(Lat) und cos(Lat) sind zb solche Kandidaten.

[Christian (Gast)]

Hallo,

jetzt hab ich mein genaues Ergebnis... ich hatte wirklich übersehen, 
dass die ganzzahligen Divisionen ja gar kein float ergeben.

Ab welcher Potenz lohnt es sich eigentlich, pow anstatt von 
Multiplizieren zu verwenden?

Vielen Dank und Grüße,
Christian

[Karl H. (kbuchegg)]

Christian schrieb:

> Ab welcher Potenz lohnt es sich eigentlich, pow anstatt von
> Multiplizieren zu verwenden?

Musst du ausprobieren.
Du kannst allerdings davon ausgehen, dass pow diese Berechnung hier 
durchführt

      exp( exponent * ln( mantisse ) )

eine andere Chance hat pow gar nicht um nicht-ganzzahlige Exponenten zu 
berechnen.

weder ln noch exp sind aber 'billige' Funktionen. In ihnen stecken 
selbst wieder eine ganze Menge Multiplikationen.

Benutze pow wenn du musst!
Praktisch alle Potenzen die in technischen Formeln vorkommen, können 
einfacher berechnet werden, indem man sich selbst eine Potenz dafür 
macht. In technischen Formenl kommen häufig Quadrate, 3te Potenzen und 
ganz selten 4te Potenzen vor. Bei all diesen bist du mit Ausschreiben 
(oder Makro oder Inline-Funktion) besser bedient (sowohl was Genauigkeit 
als auch Geschwindigkeit angeht) als mit pow.

[Andreas F. (aferber)]

Karl heinz Buchegger schrieb:
> Bei all diesen bist du mit Ausschreiben
> (oder Makro oder Inline-Funktion) besser bedient (sowohl was Genauigkeit
> als auch Geschwindigkeit angeht) als mit pow.

Und wie so oft kommt es auf den Compiler an ;-)

Gerade ausprobiert, gcc 4.3.2 für x86 mit -O2 macht tatsächlich aus 
"pow(a, 2.0)" eine einfache Multiplikation:

   3:   dd 45 08                fldl   0x8(%ebp)
    return pow(a, 2.0);
   6:   d8 c8                   fmul   %st(0),%st

Bei 3.0 ruft er aber die Library-Funktion auf.

Andreas

[Karl H. (kbuchegg)]

Andreas Ferber schrieb:

> Gerade ausprobiert, gcc 4.3.2 für x86 mit -O2 macht tatsächlich aus
> "pow(a, 2.0)" eine einfache Multiplikation:

Jetzt schämm ich mich.
Dabei bin ich es doch, der dauernd predigt: Lass den Optimizer seinen 
Job machen. :-)