Minimising cumulative floating point arithmetic error

Question

I have a 2D convex polygon in 3D space and a function to measure the area of the polygon. ``` public double area() { if (vertices.size() >= 3) { double area = 0; Vector3 origin = vertices.get(0); Vector3 prev = vertices.get(1).clone(); prev.sub(origin); for (int i = 2; i < vertices.size(); i++) { Vector3 current = vertices.get(i).clone(); current.sub(origin); Vector3 cross = prev.cross(current); area += cross.magnitude(); prev = current; } area /= 2; return area; } else { return 0; } } ``` To test that this method works at all orientations of the polygon I had my program rotate it a little bit each iteration and calculate the area. Like so... ``` 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 + ":	" + f.area()); ``` My method seems correct when testing with a 20x20 square. However the rotate method (a method in the Vector3 class) seems to introduce some error…

Accepted Answer

You'll always have some cumulative error with repeated floating point trig operations — that's just how they work. To deal with it, you basically have two options:

Just ignore it. Note that, in your example, after 20,000 iterations(!) the area is still accurate down to 13 decimal places. That's not bad, considering that doubles can only store about 16 decimal places to begin with.

Indeed, plotting your graph, the area of your square seems to be going down more or less linearly: This makes sense, assuming that the effective determinant of your approximate rotation matrix is about 1 − 3.417825 × 10-18, which is well within normal double precision floating point error range of one. If that's the case, the area of your square would continue a very slow exponential decay towards zero, such that you'd need about two billion (2 × 109) 7.3 × 1014 iterations to get the area down to 399. Assuming 100 iterations per second, that's about seven and a half months 230 thousand years.

Edit: When I first calculated how long it would take for the area to reach 399, it seems I made a mistake and somehow managed to overestimate the decay rate by a factor of about 400,000(!). I've corrected the mistake above.

If you still feel you don't want any cumulative error, the answer is simple: don't iterate floating point rotations. Instead, have your object store its current orientation in a member variable, and use that information to always rotate the object from its original orientation to its current one.

This is simple in 2D, since you just have to store an angle. In 3D, I'd suggest storing either a quaternion or a matrix, and occasionally rescaling it so that its norm / determinant stays approximately one (and, if you're using a matrix to represent the orientation of a rigid body, that it remains approximately orthogonal).

Of course, this approach won't eliminate cumulative error in the orientation of the object, but the rescaling does ensure that the volume, area and/or shape of the object won't be affected.

Minimising cumulative floating point arithmetic error

Solution

Problem