0

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!

Gabriel Dante
  • 392
  • 1
  • 12
  • 3
    Yeah don;t use that website to learn anything honestly. For starters **never** use ``#include `` – user438383 Oct 05 '22 at 09:26
  • Why isn't this include recommended, and neither is this website? – Gabriel Dante Oct 05 '22 at 09:28
  • The printf would need to use `%lf` too. Could it be that? – Joop Eggen Oct 05 '22 at 09:30
  • 2
    https://stackoverflow.com/questions/31816095/why-should-i-not-include-bits-stdc-h – VLL Oct 05 '22 at 09:30
  • @GabrielDante Why do you use `printf` in C++? – Daniel Langr Oct 05 '22 at 09:30
  • 1
    @JoopEggen It's the same thing as `%f`. – HolyBlackCat Oct 05 '22 at 09:32
  • @DanielLangr Perhaps you should ask the person who wrote the tutorial linked in the post. – VLL Oct 05 '22 at 09:33
  • @DanielLangr is there any reason I shouldn't? I use it simply for convenience (I'm more used to the syntax). I also tried printing the result using `std::cout`, but the error still persists. – Gabriel Dante Oct 05 '22 at 09:37
  • @HolyBlackCat thanks I am more at home with java. Is it not that a float arg takes 4 bytes on the stack, or are they pushed as double? How does printf know its vararg is float resp. double? – Joop Eggen Oct 05 '22 at 09:39
  • 1
    @JoopEggen Passing a `float` to a C-style variadic function converts it to a `double`. Similarly, passing integers smaller than an `int` converts them to `int`. Both `%f` and` %lf` thus expect a `double`. – HolyBlackCat Oct 05 '22 at 09:40
  • @GabrielDante Well, `printf` is a C library function. In C++, one should prefer what C++ itself provides. It's just a matter of good practice. Same as not using `#include ` and `using namespace std;`. – Daniel Langr Oct 05 '22 at 09:48

1 Answers1

7

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).

Bob__
  • 12,361
  • 3
  • 28
  • 42