Minimierung des kumulativen Fehlers bei der Gleitkomma-Arithmetik

Question

Ich habe ein 2D konvexes Polygon im 3D-Raum und eine Funktion zur Messung der Fläche des Polygons.

public double Fläche() {
    if (vertices.size() >= 3) {
        double Fläche = 0;
        Vector3 Ursprung = vertices.get(0);
        Vector3 vorher = vertices.get(1).clone();
        vorher.sub(Ursprung);
        for (int i = 2; i < vertices.size(); i++) {
            Vector3 aktuell = vertices.get(i).clone();
            aktuell.sub(Ursprung);
            Vector3 kreuz = vorher.cross(aktuell);
            Fläche += kreuz.magnitude();
            vorher = aktuell;
        }
        Fläche /= 2;
        return Fläche;
    } else {
        return 0;
    }
}

Um zu testen, ob diese Methode bei allen Ausrichtungen des Polygons funktioniert, habe ich mein Programm so gemacht, dass es bei jeder Iteration das Polygon ein wenig dreht und die Fläche berechnet. Wie folgt...

Face f = poly.getFaces().get(0);
        for (int i = 0; i < f.size(); i++) {
            Vector3 v = f.getVertex(i);
            v.rotate(0.1f, 0.2f, 0.3f);
        }
        if (blah % 1000 == 0)
            System.out.println(blah + ":\t" + f.area());

Meine Methode scheint korrekt zu sein, wenn ich mit einem 20x20 Quadrat teste. Allerdings scheint die rotate-Methode (eine Methode in der Vector3-Klasse) etwas Fehler in die Position jedes Vertex im Polygon einzuführen, was die Flächenberechnung beeinflusst. Hier ist die Vector3.rotate() Methode

public void rotate(double xWinkel, double yWinkel, double zWinkel) {
    double altY = y;
    double altZ = z;
    y = altY * Math.cos(xWinkel) - altZ * Math.sin(xWinkel);
    z = altY * Math.sin(xWinkel) + altZ * Math.cos(xWinkel);

    altZ = z;
    double altX = x;
    z = altZ * Math.cos(yWinkel) - altX * Math.sin(yWinkel);
    x = altZ * Math.sin(yWinkel) + altX * Math.cos(yWinkel);

    altX = x;
    altY = y;
    x = altX * Math.cos(zWinkel) - altY * Math.sin(zWinkel);
    y = altX * Math.sin(zWinkel) + altY * Math.cos(zWinkel);
}

Hier ist die Ausgabe für mein Programm im Format "Iteration: Fläche":

0:  400.0
1000:   399.9999999999981
2000:   399.99999999999744
3000:   399.9999999999959
4000:   399.9999999999924
5000:   399.9999999999912
6000:   399.99999999999187
7000:   399.9999999999892
8000:   399.9999999999868
9000:   399.99999999998664
10000:  399.99999999998386
11000:  399.99999999998283
12000:  399.99999999998215
13000:  399.9999999999805
14000:  399.99999999998016
15000:  399.99999999997897
16000:  399.9999999999782
17000:  399.99999999997715
18000:  399.99999999997726
19000:  399.9999999999769
20000:  399.99999999997584

Weil dies letztendlich für eine Physik-Engine gedacht ist, möchte ich wissen, wie ich den kumulativen Fehler minimieren kann, da die Vector3.rotate() Methode sehr regelmäßig verwendet wird.

Danke!

Ein paar seltsame Anmerkungen:

• Der Fehler ist proportional zur gedrehten Menge. d.h. größere Rotation pro Iteration -> größerer Fehler pro Iteration.

• Es gibt mehr Fehler, wenn Double-Werte an die rotate-Funktion übergeben werden als wenn Float-Werte übergeben werden.

Answer 1

Sie werden immer einige Kumulativfehler mit wiederholten Gleitkomma-Trig-Operationen haben - so funktionieren sie einfach. Um damit umzugehen, haben Sie grundsätzlich zwei Möglichkeiten:

1. Einfach ignorieren. Beachten Sie, dass im Beispiel nach 20.000 Iterationen(!) der Bereich noch genau bis auf 13 Dezimalstellen genau ist. Das ist nicht schlecht, wenn man bedenkt, dass Double nur etwa 16 Dezimalstellen speichern können.

   Tatsächlich scheint die Fläche Ihres Quadrats beim Plotten Ihres Graphen mehr oder weniger linear abzunehmen:
   
   Das ergibt Sinn, unter der Annahme, dass die effektive Determinante Ihrer ungefähren Rotationsmatrix etwa 1 − 3,417825 × 10-18 beträgt, was gut innerhalb des normalen Bereichs des Gleitkommafehlers von eins liegt. Wenn das der Fall ist, würde die Fläche Ihres Quadrates weiterhin exponentiell langsam gegen Null abnehmen, so dass Sie ungefähr zwei Milliarden (2 × 109) 7,3 × 1014 Iterationen benötigen, um die Fläche auf 399 zu reduzieren. Unter der Annahme von 100 Iterationen pro Sekunde sind das ungefähr sieben und einhalb Monate 230 Tausend Jahre.

   Bearbeiten: Als ich zuerst berechnet habe, wie lange es dauern würde, bis die Fläche 399 erreicht, scheint es, dass ich einen Fehler gemacht und es irgendwie geschafft habe, die Zersetzungsgeschwindigkeit um einen Faktor von etwa 400.000(!) zu überschätzen. Ich habe den oben genannten Fehler korrigiert.

2. Wenn Sie dennoch das Gefühl haben, dass Sie keinen Kumulativfehler möchten, ist die Antwort einfach: Iterieren Sie keine Gleitkomma-Rotationen. Speichern Sie stattdessen die aktuelle Ausrichtung Ihres Objekts in einer Member-Variablen und verwenden Sie diese Informationen, um das Objekt immer von seiner ursprünglichen Ausrichtung auf seine aktuelle zu drehen.

   Dies ist einfach in 2D, da Sie nur einen Winkel speichern müssen. In 3D würde ich vorschlagen, entweder einen Quaternion oder eine Matrix zu speichern und gelegentlich neu zu skalieren, damit seine Norm / Determinante ungefähr eins bleibt (und, wenn Sie eine Matrix verwenden, um die Ausrichtung eines starren Körpers darzustellen, dass sie ungefähr orthogonal bleibt).

   Natürlich beseitigt dieser Ansatz den Kumulativfehler in der Ausrichtung des Objekts nicht, aber das Neuskalieren gewährleistet, dass das Volumen, die Fläche und/oder Form des Objekts nicht beeinträchtigt werden.

Answer 2

Du sagst, es gibt kumulative Fehler, aber ich glaube nicht, dass es welche gibt (beachte, wie deine Ausgabe nicht immer nach unten geht) und der Rest des Fehlers entsteht nur durch Rundung und Verlust an Genauigkeit in einem Float.

Ich habe an einem 2D-Physik-Engine an der Universität gearbeitet (in Java) und fand double genauer zu sein (natürlich ist es so, siehe Oracles Datentypgrößen

Kurz gesagt, man wird dieses Verhalten nie loswerden, man muss nur die Genauigkeitsgrenzen akzeptieren

EDIT:

Wenn ich mir nun deine .area Funktion anschaue, gibt es möglicherweise kumulative Fehler aufgrund von

+= cross.magnitude

aber ich muss sagen, dass die ganze Funktion etwas seltsam aussieht. Warum muss sie die vorherigen Vertices kennen, um die aktuelle Fläche zu berechnen?