CODEWITHSUNDEEP
c++
Twigg
Twigg

Reputation: 13

C++ Trapezoidal Integration Function Returning Negative Numbers when it shouldn't

I am using the following function written in C++, whose purpose is to take the integral of one array of data (y) with respect to another (x)

// Define function to perform numerical integration by the trapezoidal rule
double trapz (double xptr[], double yptr[], int Npoints)
{
    // The trapzDiagFile object and associated output file are how I monitor what data the for loop actually sees.
    std::ofstream trapzDiagFile;
    trapzDiagFile.open("trapzDiagFile.txt",std::ofstream::out | std::ofstream::trunc);

    double buffer = 0.0;
    for (int n = 0; n < (Npoints - 1); n++)
    {
        buffer += 0.5 * (yptr[n+1] + yptr[n]) * (xptr[n+1] - xptr[n]);
        trapzDiagFile << xptr[n] << "," << yptr[n] << std::endl;
    }

    trapzDiagFile.close();
    return buffer;
}

I validated this function for the simple case where x contains 100 uniformly spaced points from 0 to 1, and y = x^2, and it returned 0.33334, as it should.
But when I use it for a different data set, it returns -3.431, which makes absolutely no sense. If you look in the attached image file, the integral I am referring to is the area under the curve between the dashed vertical lines.
It's definitely a positive number.
Moreover, I used the native trapz command in MATLAB on the same set of numbers and that returned 1.4376.
In addition, I translated the above C++ trapz function into MATLAB, line for line as closely as possible, and again got 1.4376.
I feel like there's something C++ related I'm not seeing here. If it is relevant, I am using minGW-w64.

Apologies for the vagueness of this post. If I knew more about what kind of issue I am seeing, it would be easier to be concise about it.

Plot of the dataset for which the trapz function (my homemade C++ version) returns -3.431:

enter image description here

Upvotes: 1

Views: 760

Answers (1)

Owen Delahoy
Owen Delahoy

Reputation: 805

Please check the value of xptr[Npoints - 1]. It may be less than xptr[Npoints - 2], and was not included in the values that you output.

Upvotes: 1

Related Questions