c - 双浮点精度

Question

我正在用高斯-塞德尔方法求解拉普拉斯方程，但在某些区域，它呈现出类似平台的现象。形式上，即通过数值分析，这样的区域不应该存在，即使梯度几乎为零。

我被迫相信 double 不足以执行算术运算，需要使用大量数字库(降低性能，因为现在它将由软件完成)。或者，我应该以不同的顺序执行操作，旨在保留小数的一些意义。

例子

Cell (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.9999999999999966。

代码

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

Answer 1

因为你正在解决:

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

相反，您可以解决等效问题:

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

其中，phi' = phi - 1。

请记住根据 phi' 应用边界条件。

最后在解收敛后，可以得到解为phi = 1 + phi'。

我在这里假设边界值接近 1。

我没试过这个，但我认为数字将以浮点符号的有效数字表示，因此截断错误会减少。