解决浮点问题C

概述我开发了一个科学应用(模拟染色体在细胞核中移动).染色体被划分成使用4×4旋转矩阵围绕随机轴旋转的小碎片. 问题是仿真执行了数百亿次旋转,因此浮点舍入误差叠加并呈指数增长,因此随着时间的流逝,碎片往往会“浮起”并脱离染色体的其余部分. 我用C双精度.软件在CPU上运行,但将被移植到CUDA,模拟可以持续1个月最多. 我不知道我可以以某种方式重新归一化染色体,因为所有的片段都链接在一起(你可以看到它 我开发了一个科学应用(模拟染色体在细胞核中移动).染色体被划分成使用4×4旋转矩阵围绕随机轴旋转的小碎片.

问题是仿真执行了数百亿次旋转,因此浮点舍入误差叠加并呈指数增长,因此随着时间的流逝,碎片往往会“浮起”并脱离染色体的其余部分.

我用C双精度.软件在cpu上运行,但将被移植到CUDA,模拟可以持续1个月最多.

我不知道我可以以某种方式重新归一化染色体,因为所有的片段都链接在一起(你可以看到它是一个双重链接的列表),但我认为这将是最好的想法,如果可能的话.

你有什么建议吗 ？我觉得有点迷失了.

非常感谢你,

H.

编辑：
添加了一个简化的示例代码.
你可以假定所有的矩阵数学都是古典的实现.

// Rotate 1000000 timesfor (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;    }}

解决方法 您需要为您的系统找到一些约束,并努力将其保持在一定的合理范围内.我做了一系列的分子碰撞模拟,在这些系统中,总能量是保守的,所以每一步我都仔细检查系统的总能量,如果它变化一些阈值,那么我知道我的时间步长很差(太大或太小),我选择一个新的时间步骤并重新运行.这样我就可以实时跟踪系统发生了什么.

对于这种模拟,我不知道你有多少保存的数量,但是如果你有一个,你可以尝试保持这个常数.记住,让你的时间步长更小,并不总是提高精度,你需要用你所拥有的精度来优化步长.我已经进行了数周模拟运行了几个星期的cpu时间,保存的数量总是在10分之1的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,0) + (dx,dy) -> (x0,y0) (dx1,dy1)t2 (x0,dy1) + (dx,y0) (dx2,dy2)t3 (x0,dy2) + (dx,y0) (dx3,dy3)

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

现在简单的翻译像我的例子,你可能不会赢得很多.但是为了旋转,这可以是一个救命.只需保留与每个染色体相关联的欧拉角度的矩阵,并更新而不是染色体的实际位置.至少这样他们不会漂浮.

总结

以上是内存溢出为你收集整理的解决浮点问题C全部内容，希望文章能够帮你解决解决浮点问题C所遇到的程序开发问题。

如果觉得内存溢出网站内容还不错，欢迎将内存溢出网站推荐给程序员好友。