Why does using double instead of float gives a wrong result in this double integration code?-CodePudding

I found the following code in this page to compute a double integral. Whenever I run it with all variables being declared as float, it gives the right result for the example integral, which is 3.91905. However, if I just change all float variables to double, the program gives a completely wrong result (2.461486) for this integral.

Could you help me undestanding why this happens? I expected to have a better result using double precision, but that's evidently not the case here.

Below is the code pasted from the aforementioned website.

// 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;
}

Thanks in advance for your help!

CodePudding user response：

Due to numeric imprecisions, this line:

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

Evaluates to 5 when the types of uy, ly and k are float. When the type is double, it yields 4.

You may use std::round((uy - ly) / k) or a different formula (I haven't checked the mathematical correctness of the whole program).