最小化累积浮点算术错误

Question

我在3D空间中有一个2D凸多边形和一个测量多边形面积的函数.

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;
    }
}

为了测试这个方法适用于多边形的所有方向,我的程序每次迭代都会旋转一点并计算面积.像这样......

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());

使用20x20平方测试时,我的方法似乎是正确的.然而,旋转方法(Vector3类中的方法)似乎在多边形中的每个顶点的位置中引入了一些误差,这会影响面积计算.这是Vector3.rotate()方法

public void rotate(double xAngle, double yAngle, double zAngle) {
    double oldY = y;
    double oldZ = z;
    y = oldY * Math.cos(xAngle) - oldZ * Math.sin(xAngle);
    z = oldY * Math.sin(xAngle) + oldZ * Math.cos(xAngle);

    oldZ = z;
    double oldX = x;
    z = oldZ * Math.cos(yAngle) - oldX * Math.sin(yAngle);
    x = oldZ * Math.sin(yAngle) + oldX * Math.cos(yAngle);

    oldX = x;
    oldY = y;
    x = oldX * Math.cos(zAngle) - oldY * Math.sin(zAngle);
    y = oldX * Math.sin(zAngle) + oldY * Math.cos(zAngle);
}

这是我的程序的输出格式为"iteration:area":

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

由于这最初是针对物理引擎的,我想知道如何最小化累积误差,因为Vector3.rotate()方法将定期使用.

谢谢!

几个奇怪的笔记:

• 误差与旋转量成比例.即.每次迭代更大的旋转 - >每次迭代更大的错误.

• 将双精度数传递给旋转函数时比传递浮点数时出现的错误更多.

Answer 1

对于重复的浮点触发操作,您总会遇到一些累积错误 - 这就是它们的工作原理.要处理它,你基本上有两个选择:

1. 只是忽略它.请注意,在您的示例中,在20,000次迭代(!)之后,该区域仍然精确到13位小数.考虑到双打只能存储大约16位小数,这还不错.

   实际上,绘制图表时,方块的面积似乎或多或少线性下降:
   
   这是有道理的,假设近似旋转矩阵的有效行列式约为1 - 3.417825×10 ^-18,这完全在1的正常双精度浮点误差范围内.如果是这种情况,你的方块的面积将继续一个非常缓慢的指数衰减为零,这样你就需要二十亿(2×10 ⁹)7.3×10 ¹⁴次迭代使区域降至399.假设每秒100次迭代,那就是七个半月 23万年.

   编辑:当我第一次计算出该区域达到399所需的时间时,似乎我犯了一个错误,并且设法将衰减​​率高估了大约400,000(!).我纠正了上面的错误.

2. 如果您仍然觉得不需要任何累积错误,答案很简单:不要迭代浮点旋转.相反,让对象将其当前方向存储在成员变量中,并使用该信息始终将对象从其原始方向旋转到当前方向.

   这在2D中很简单,因为您只需要存储一个角度.在3D中,我建议存储四元数或矩阵,偶尔重新缩放它以使其范数/行列式保持大约一(并且,如果你使用矩阵来表示刚体的方向,它仍然存在近似正交).

   当然,这种方法不会消除对象方向的累积误差,但重新缩放确保不会影响对象的体积,面积和/或形状.