解决浮点舍入问题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,你可以这样做

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)

所以在任何时候,tn,要获得你的真实位置,你必须做(x0,y0)+(dxn,dyn)

现在像我的例子一样简单的翻译,你可能不会赢得很多.但对于轮换,这可以节省生命.只需保留一个与每个染色体相关的欧拉角矩阵,然后更新它而不是染色体的实际位置.至少这种方式他们不会飘走.