解决浮点舍入问题C++

Question

我开发了一个科学应用程序（模拟在细胞核中移动的染色体）。染色体被分成小片段，使用4x4旋转矩阵围绕随机轴旋转。

问题是模拟执行了数千亿次的旋转，因此浮点舍入误差会堆积并呈指数级增长，所以碎片会随着时间流逝而“飘离”并与染色体的其余部分分离。

我使用C++的双精度。目前在CPU上软运行，但将被移植到CUDA，并且仿真最多可持续1个月。因为所有的片段都被链接在一起（你可以看到它是一个双向链表），但我认为如果可能的话，这将是最好的想法。

你有什么建议吗？我感觉有点失落。

非常感谢你，

H.

编辑： 增加了一个简单的示例代码。 你可以假设所有的矩阵数学都是经典的实现。

// Rotate 1000000 times 
for (int i = 0; i < 1000000; ++i) 
{ 
    // Pick a random section start 
    int istart = rand() % chromosome->length; 

    // Pick the end 20 segments further (cyclic) 
    int iend = (istart + 20) % chromosome->length; 

    // Build rotation axis 
    Vector4 axis = chromosome->segments[istart].position - chromosome->segments[iend].position; 
    axis.normalize(); 

    // Build rotation matrix and translation vector 
    Matrix4 rotm(axis, rand()/float(RAND_MAX)); 
    Vector4 oldpos = chromosome->segments[istart].position; 

    // Rotate each segment between istart and iend using rotm 
    for (int j = (istart + 1) % chromosome->length; j != iend; ++j, j %= chromosome->length) 
    { 
     chromosome->segments[j].position -= oldpos; 
     chromosome->segments[j].position.transform(rotm); 
     chromosome->segments[j].position += oldpos; 
    } 
} 

Answer 1

您需要为您的系统找到一些约束，并努力将其保持在合理范围内。我做了大量的分子碰撞模拟，并且在这些系统中总能量是守恒的，所以每一步我都要仔细检查系统的总能量，如果它的变化达到某个阈值，那么我知道我的时间步长选择的很差（太大或太小），我选择一个新的时间步并重新运行。这样我可以实时跟踪系统发生的情况。

对于这个模拟，我不知道你有多少守恒量，但是如果你有一个守恒量，你可以试着保持不变。请记住，缩短时间步长并不总是会提高准确度，您需要根据精确度来优化步长。我已经进行了数周的CPU时间数值模拟，并且守恒数量总是在10^8的1个部分内，所以有可能，你只需要玩一些。另外，正如Tomalak所说，也许试着总是引用你的系统到开始时间，而不是上一步。因此，不是总是移动你的染色体，而是将染色体保存在其起始位置，并与它们一起存储一个转换矩阵，从而使你到达当前位置。当你计算你的新旋转时，只需修改变换矩阵。它可能看起来很愚蠢，但有时这很好，因为错误的平均值为0.

例如，假设我有一个位于（x，y）处的粒子和我计算的每个步骤（dx，dy）和移动粒子。分步的方式将做到这一点

t0 (x0,y0) 
t1 (x0,y0) + (dx,dy) -> (x1, y1) 
t2 (x1,y1) + (dx,dy) -> (x2, y2) 
t3 (x2,y2) + (dx,dy) -> (x3, y3) 
t4 (x3,30) + (dx,dy) -> (x4, y4) 
... 

如果你总是参考T0，你可以在任何时候，TN做到这一点

t0 (x0, y0) (0, 0) 
t1 (x0, y0) (0, 0) + (dx, dy) -> (x0, y0) (dx1, dy1) 
t2 (x0, y0) (dx1, dy1) + (dx, dy) -> (x0, y0) (dx2, dy2) 
t3 (x0, y0) (dx2, dy2) + (dx, dy) -> (x0, y0) (dx3, dy3) 

所以，要想让你的，你需要做的实际位置（ x0，y0）+（dxn，dyn）

现在为我的例子简单的翻译，你可能不会赢得很多。但是对于轮换而言，这可以成为一种生活救星。只需保留与每个染色体相关联的欧拉角的矩阵并更新该矩阵，而不是染色体的实际位置。至少这样他们不会飘走。

Answer 2

写您的公式，这样的时间步T数据不仅仅从浮点数据的时间步长T-1派生。尽量确保浮点错误的产生仅限于一个时间步。

在没有更具体的问题解决的情况下，很难说出更具体的内容。

Answer 3

我认为这取决于你使用的编译器。

的Visual Studio编译器支持/ fp的开关，它告诉浮点运算

可以read more about it的行为。基本上，/ fp：严格是最苛刻的模式

Answer 4

我想这取决于所需的精度，但您可以使用基于'整数'的浮点数。通过这种方法，您可以使用整数并为小数位数提供您自己的偏移量。

例如，与4个小数点精确，你将有

浮点值 - > int值 1.0000 - > 10000 1.0001 - > 10001 0.9999 - > 09999

你必须当你进行乘法和除法时要小心，当你应用精度偏移时要小心。其他方面，你可以很快得到溢出错误。

1.0001 * 1.0001变为10001 *一万分之一万〇一

Answer 5

问题描述相当含糊，所以这里有一些比较模糊的建议。

选项1：

寻找一些约束集，从而（1）应始终持有，（2），如果他们失败了，但只是刚刚，可以很容易地调整系统，使他们这样做，（3 ）如果他们全都成立，那么你的模拟不会太疯狂，以及（4）当系统开始疯狂时，约束开始失败，但只是轻微。例如，或许染色体相邻位之间的距离应该至多为d，对于某些d，并且如果少数距离略大于d，则可以（例如）从一端沿着染色体行走，修复通过将下一个片段及其所有后继片段移动到其前一个片段的距离太大。或者其他的东西。

然后经常检查约束条件以确保任何违规行为在被捕获时仍然很小;当你发现违规行为时，请纠正。 （你应该安排，当你修理东西时，你“不仅仅是满足”约束条件）。

如果检查限制的时间很长，那么当然你可以做到这一点。 （这样做可以使你做的修正更便宜，例如，如果这意味着任何违反总是很小的。）

选项2：

查找描述，使系统的状态的一个新途径这个问题不可能出现。例如，也许（我怀疑这一点），你可以为每个相邻的碎片对存储一个旋转矩阵，并强制它总是一个正交矩阵，然后让这些碎片的位置由这些旋转矩阵隐式确定。

方案3：

而不是你的约束思维的限制，提供一些小的“复原力”，这样，当某样东西脱节也容易被拉向后朝事情应该是这样。请注意，没有任何问题时，恢复力为零或至少可以忽略不计，以免它们比原始数字错误更严重地干扰您的结果。

Answer 6

如果我正确读取此代码，任何时候任何两个相邻染色体片段之间的距离都应该改变。在这种情况下，在主循环计算每对相邻点之间和主循环之后的距离之前，必要时移动每个点，以便与前一点具有适当的距离。

根据具体情况，您可能需要在主循环中多次强制执行此约束。

Answer 7

基本上，您需要避免来自这些（不精确的）矩阵运算符的错误累积，并且在大多数应用程序中有两种主要方法。

1. 而不是写的位置上多次操作的一些初始位置，你可以写出来，究竟该运营商将个运算后明确。例如，假设你有一个位置x，并且你正在增加一个值e（你不能完全表示）。比计算x + = e好得多;大量的时间将是计算x + EN;其中EN是一些更准确的方式来表示N次后的操作。你应该想一下，你是否有更精确地表示许多旋转动作的方法。
2. 稍微有些人为的是把你的新发现的点和项目从你的旋转中心预期的半径差异。这将保证它不会漂移（但不一定能保证旋转角度准确。）