为什么在此双积分代码中使用 double 而不是 float 会给出错误的结果？

Question

我在此页面中找到了以下代码来计算二重积分。每当我在所有变量都声明为 的情况下运行它时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;
}

在此先感谢您的帮助！

Answer 1

由于数字不精确，这一行：

ny = (uy - ly) / k + 1;  // 'ny' is an int.

uy当、ly和的类型k为 时，计算结果为 5 float。当类型为 时double，它产生 4。

您可以使用std::round((uy - ly) / k)或不同的公式（我没有检查整个程序的数学正确性）。