c - 双浮点精度

我用高斯-赛德尔方法解拉普拉斯方程，但在某些地区，它显示出一个类似高原的方面。形式上，即通过数值分析，即使梯度几乎为零，也不应该存在这些区域。
我不得不相信，双精度不足以执行算法，需要使用大量的库（扼杀性能，因为现在它将由软件完成）。或者，我应该按照不同的顺序进行操作，目的是保留小数点的一些意义。
例子
单元（13、14、0）正在由7点网格（在3D中）更新，其邻域为：

(12,14,0)=  0.9999999999999936; // (x-)
(14,14,0)=  0.9999999999999969; // (x+)
(13,13,0)=  0.9999999999999938; // (y-)
(13,15,0)=  1.0000000000000000; // (y+)
(13,14,-1)= 1.0000000000000000; // (z-)
(13,14,1)=  0.9999999999999959; // (z+)

因此，单元（13,14,0）的新值将被评估为：

p_new = (0.9999999999999936 + 0.9999999999999969 + 0.9999999999999938 + 1.0000000000000000 + 1.0000000000000000 + 0.9999999999999959) / 6.0 ;

这导致p_new是1.0000000000000000，而它应该是0.99999999999666。
代码

#include <stdio.h>

int main()
{
    double ad_neighboor[6] = {0.9999999999999936, 0.9999999999999969,
                              0.9999999999999938, 1.0000000000000000,
                              1.0000000000000000, 0.9999999999999959};

    double d_denom = 6.0;

    unsigned int i_xBackward=0;
    unsigned int i_xForward=1;

    unsigned int i_yBackward=2;
    unsigned int i_yForward=3;

    unsigned int i_zBackward=4;
    unsigned int i_zForward=5;

    double d_newPotential = (ad_neighboor[i_xForward] + ad_neighboor[i_xBackward] +
                             ad_neighboor[i_yForward] + ad_neighboor[i_yBackward] +
                             ad_neighboor[i_zForward] + ad_neighboor[i_zBackward] ) / d_denom;

    printf("%.16f\n", d_newPotential);
}

最佳答案

既然你在解决：

d²(phi)/dx² + d²(phi)/dy² = 0

相反，你可以解决同样的问题：

d²(phi')/dx² + d²(phi')/dy² = 0

其中，phi' = phi - 1。
记住应用phi'中的边界条件。
最后，在解收敛后，可以得到phi = 1 + phi'的解。
我假设这里的边界值接近1。
我还没有试过，但我认为这些数字将用它们的有效数字用浮点数表示，从而减少截断误差。

关于c - 双浮点精度，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/48464876/