Find the area between two functions

Question

I solve the following problem.

There are two functions f(x) and g(x). Each of them is divided into parts. The f(x) function consists of 'n' parts, and g(x) consists of 'm' parts. Each part is represented as a Second-Degree Polynomial Function (ax^2 + bx + c). We need to find an area between f(x) and g(x) with precision 10^(-6).

Input data:

• n, m (1 <= n, m <= 10^5);

• boundary points of f(x); //n+1 points

• polynomials of f(x); //n strings (a, b, c for ax^2 + bx + c)

• boundary points of g(x); //m+1 points

• polynomials of g(x); //m strings (a, b, c for ax^2 + bx + c)

The first and last points of the functions are equal.

Sample input:

1 1
0 1
1 -2 1
0 1
-1 2 1

Output:

1.3333333333

Here is my code:

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


long double integrate(int f[3], int g[3], long double left, long double right) {
    long double a = f[0] - g[0], b = f[1] - g[1], c = f[2] - g[2];
    long double up = (a * right * right * right) / 3 + (b * right * right) / 2 + c * right;
    long double down = (a * left * left * left) / 3 + (b * left * left) / 2 + c * left;
    return (up > down) ? up - down : down - up;
}

void solve(int f[3], int g[3], int left, int right, long double &sum) {

    long double a = f[0] - g[0], b = f[1] - g[1], c = f[2] - g[2];
    if (a == 0) {
        if (b != 0) {
            long double x = -c/b;
            if (x > left && x < right) {
                sum += integrate(f, g, left, x);
                sum += integrate(f, g, x, right);
            }
        } else {
            sum += integrate(f, g, left, right);
        }
        return;
    }

    long double discriminant = b * b - 4 * a * c;
    if (discriminant < 0) { sum += integrate(f, g, left, right); }
    else {
        long double q = b >= 0 ? (-b - sqrt(discriminant))/2 : (-b + sqrt(discriminant))/2;
        long double x1 = q / a, x2 = c / q;
        if (discriminant == 0.0) {
            if (x1 > left && x1 < right) {
                sum += integrate(f, g, left, x1);
                sum += integrate(f, g, x1, right);
            } else {
                sum += integrate(f, g, left, right);
            }
        } else {
            long double first = min(x1, x2), second = max(x1, x2);
            if (first > left && first < right) {
                sum += integrate(f, g, left, first);
                if (second > left && second < right) {
                    sum += integrate(f, g, first, second);
                    sum += integrate(f, g, second, right);
                } else {
                    sum += integrate(f, g, first, right);
                }
            } else if (second > left && second < right) {
                sum += integrate(f, g, left, second);
                sum += integrate(f, g, second, right);
            } else {
                sum += integrate(f, g, left, right);
            }
        }
    }
    return;
}


int main() {
//freopen("input.txt", "r", stdin);
//freopen("output.txt", "w", stdout);

    int n, m;
    cin >> n >> m;

    int f[n+1];
    int ffunc[n][3];
    int g[m+1];
    int gfunc[m][3];
    set <int> points;

    for (int i = 0; i < n+1; i++) {
        cin >> f[i];
        points.insert(f[i]);
    }

    for (int i = 0; i < n; i++) {
        for (int k = 0; k < 3; k++) {
            cin >> ffunc[i][k];
        }
    }

    for (int i = 0; i < m+1; i++) {
        cin >> g[i];
        points.insert(g[i]);
    }

    for (int i = 0; i < n; i++) {
        for (int k = 0; k < 3; k++) {
            cin >> gfunc[i][k];
        }
    }

    int fit = 0, git = 0;
    long double sum = 0.0;
    auto it1 = points.begin();
    auto it2 = points.begin();
    it2++;
    while (it2 != points.end()) {

        solve(ffunc[fit], gfunc[git], *it1, *it2, sum);

        if (f[fit+1] == *it2) {
            fit++;
        }
        if (g[git+1] == *it2) {
            git++;
        }
        it1++;
        it2++;
    }
    cout.precision(27);
    cout << sum;
    return 0;
}

The program doesn't pass some tests. It gets wrong answer.

What could be the problem?

Answer 1

In order to find the area between two curves you have to do the following:

1. Find the points at which the two curves intersect, these separate the x axis into segments.
2. For each segment, compute the integral for each of the two curves, this gives the area under each curve for this segment.
3. In each segment, you have one curve that is above the other and one that is below it. Subtract the area of the below function from the area of the above function. That gives the area between the curves.
4. Sum up the results you found for each segment.

This needs to slightly adapted in case the functions are not both non-negative. Luckily, all your functions are polynomials, so can write the anti-derivative in explicit form.

Answer 2

Area under curve (area between curve and zero abscissa) and two vertical lines is an integral of function on range (definite integral).

1. First you have to determine all unique ranges where only one function from g(x) set is paired with one from f(x) set. In general case that number can be greater than n and m.

   In each found range area is integral of function abs(f(x)-g(x))

2. Because you have functions given, you can find points of intersections between (f) and g(x) by finding all solutions of equation.

   (af - ag) * x^2 + (bf - bg) * x + cf - cg = 0

   If determinant is negative, they don't intersect. If solution is outside of particular range, those particular fragments do not intersect.

3. The solutions are used to divide ranges further. SO theoretically range can stay unchanged or be replaced by three ranges. It's probably prudent to have a list or ranges to be able modify it easily.

4. Calculate absolute value of integral for function h(x) = f(x) - g(x) in each range. You have functions in analytical form, so you can do that analytically.

   If H(x) is antiderivative of h(x) and h(x) is monotonic (we ensured that by finding ranges between intersections) , then integral in i-th range is S[i] = H(x[i]) - H(x[i-1]) where x[i-1] and x[i] are values of x that define bounds of range.

   Precision of calculation is an odd part there, perhaps the actual requirement (omitted by OP?) is to perform integration in discrete way, again, we would better to use monotonic h(x) than abs(f(x)-g(x)) to exclude imprecision caused by intersections. With analytical solution I expect precision to be 10^-9 if float is used, though real absolute precision would much depend on the values (very small or large values of x may decrease precision). Anyway there is multitude of discrete algorithms to evaluate definite integral and precision of it.

Answer 3

In solve, there is at least a corner case not accounted for:

if (a == 0) {
    if (b != 0) {
        long double x = -c/b;
        if (x > left && x < right) {
            sum += integrate(f, g, left, x);
            sum += integrate(f, g, x, right);
        }  
        // else {
        //     sum += integrate(f, g, left, right);
        // }
    } else {
        sum += integrate(f, g, left, right);
    }
    return;
}

In main, the local variables f, ffunc, g and gfunc are all variable length arrays (a non-standard extension provided by some compilers) and should rather all be standard containers like std::vector.

The posted code also uses a std::set to store all the points and to ensure their order. This may not be the most efficient way to accomplish that result, if the boundary points are already ordered in their own ranges.

Here, a possible different implementation is tested.