将float改为double后二重积分C++代码结果错误的原因排查

我找到一段计算二重积分的C++代码，当所有变量声明为float时，示例积分能得到正确结果3.91905；但仅将所有float改为double后，结果变为错误的2.461486。原本预期double精度会带来更优结果，求解问题原因。代码如下：

// C++ program to calculate
// double integral value

#include <bits/stdc++.h>
using namespace std;

// Change the function according to your need
float givenFunction(float x, float y)
{
    return pow(pow(x, 4) + pow(y, 5), 0.5);
}

// Function to find the double integral value
float doubleIntegral(float h, float k,
                    float lx, float ux,
                    float ly, float uy)
{
    int nx, ny;

    // z stores the table
    // ax[] stores the integral wrt y
    // for all x points considered
    float z[50][50], ax[50], answer;

    // Calculating the number of points
    // in x and y integral
    nx = (ux - lx) / h + 1;
    ny = (uy - ly) / k + 1;

    // Calculating the values of the table
    for (int i = 0; i < nx; ++i) {
        for (int j = 0; j < ny; ++j) {
            z[i][j] = givenFunction(lx + i * h,
                                    ly + j * k);
        }
    }

    // Calculating the integral value
    // wrt y at each point for x
    for (int i = 0; i < nx; ++i) {
        ax[i] = 0;
        for (int j = 0; j < ny; ++j) {
            if (j == 0 || j == ny - 1)
                ax[i] += z[i][j];
            else if (j % 2 == 0)
                ax[i] += 2 * z[i][j];
            else
                ax[i] += 4 * z[i][j];
        }
        ax[i] *= (k / 3);
    }

    answer = 0;

    // Calculating the final integral value
    // using the integral obtained in the above step
    for (int i = 0; i < nx; ++i) {
        if (i == 0 || i == nx - 1)
            answer += ax[i];
        else if (i % 2 == 0)
            answer += 2 * ax[i];
        else
            answer += 4 * ax[i];
    }
    answer *= (h / 3);

    return answer;
}

// Driver Code
int main()
{
    // lx and ux are upper and lower limit of x integral
    // ly and uy are upper and lower limit of y integral
    // h is the step size for integration wrt x
    // k is the step size for integration wrt y
    float h, k, lx, ux, ly, uy;

    lx = 2.3, ux = 2.5, ly = 3.7,
    uy = 4.3, h = 0.1, k = 0.15;

    printf("%f", doubleIntegral(h, k, lx, ux, ly, uy));
    return 0;
}

原因解析 #

• 整数转换的截断误差：代码中nx = (ux - lx)/h +1和ny = (uy - ly)/k +1是核心问题。以原参数为例，(2.5-2.3)/0.1理论值为2，加1后nx应为3。但0.1是二进制无法精确表示的十进制小数：

  • float精度较低，计算后结果刚好近似为2.0，转int后得到正确的3；
  • double精度更高，计算结果会是接近2但略小于2的数值（比如1.9999999999999996），直接转int会被截断为1，加1后nx=2，采样点数量少了一个。
    采样点数量错误会完全打乱辛普森积分的权重分配（首尾点、偶次点、奇次点的系数规则），最终导致结果严重偏离正确值。

• 精度“反直觉”的本质：不是double精度差，而是它更精准地反映了二进制浮点数无法精确表示部分十进制小数的特性；float的低精度反而“碰巧”让计算结果凑成了预期的整数，掩盖了代码的逻辑漏洞。

修复方案 #

• 对浮点数运算结果做四舍五入后再转整数：

  nx = round((ux - lx)/h) + 1;
  ny = round((uy - ly)/k) + 1;
  

• 或者根据参数手动计算采样点数量（原示例中nx=3，ny=5），避免浮点数运算误差。

内容的提问来源于stack exchange，提问作者Gabriel Dante