A simple definite integrator class of a single variable in C++

Question

---

I have written a simple Integrator class in C++17 that can perform either a definite single integration of a single variable or a definite double integration of a single variable using the same integrand.

---

Here's my Integrator class:

Integrator.h

#pragma once

#include <algorithm>
#include <utility>
#include <functional>

struct Limits {
    double lower;
    double upper;

    Limits(double a = 0, double b = 0) : lower{ a }, upper{ b } {
        if (a > b) std::swap(lower, upper);
    }

    void applyLimits(double a, double b) {
        lower = a;
        upper = b;
        if (a > b) std::swap(lower, upper);
    }
};

class Integrator {
private:
    Limits limits_;
    std::function<double(double)> integrand_;

    double dx_;
    double dy_;  
    double integral_; 
    int step_size_;

public:
    Integrator(Limits limits, int stepSize, std::function<double(double)> integrand) 
        : limits_{ limits }, 
        step_size_{ stepSize }, 
        integrand_{ integrand }, 
        dx_{ 0 }, dy_{ 0 } 
    {}
    ~Integrator() = default;

    constexpr double dx() const { return this->dx_; }
    constexpr double dy() const { return this->dy_; }
    constexpr double integral() const { return this->integral_; }

    Limits limits() const { return limits_; }    
    std::function<double(double)>* integrand() { return &this->integrand_; }

    // This is always a 1st order of integration!
    constexpr double evaluate() {
        double distance = limits_.upper - limits_.lower;      // Distance is defined as X0 to XN. (upperLimit - lowerLimit) 
        dx_ = distance / step_size_;                          // Calculate the amount of iterations by dividing 
                                                              // the x-distance by the dx stepsize
        integral_ = 0;                                        // Initialize area to zero
        for (auto i = 0; i < step_size_; i++) {               // For each dx step or iteration calculate the area at Xi
            dy_ = integrand_(limits_.lower + i * dx_);
            double area = dy_ * dx_;                          // Where the width along x is defines as dxStepSize*i 
            integral_ += area;                                // and height(dy) is f(x) at Xi. Sum all of the results
        }

        return integral_;
    }

    // This will perform a second order of integration where the inner limits are defined
    // by [lower, y] where "upper" is not used directly. This may be expanded in the future...
    double integrate(double lower = 0.0, double upper = 0.0) {
        // Since we are not using the inner upper limit directly
        // make sure that it is still greater than the lower limit
        if (upper <= lower) {
            upper = lower + 1;
        }

        // As the code currently stands this temporary is not necessary as I could have
        // used the values from the arguments directly, but I wanted to keep it
        // for consistency reasons as this might be expanded in the future where the use
        // of the upper bound inner limit will be taken into context. 
        Limits limits(lower, upper);

        double outerSum = 0;
        dy_ = static_cast<double>(limits_.upper - limits_.lower) / step_size_;

        for (int i = 0; i < step_size_; i++) {
            double yi = limits_.lower+i*dy_;
            double dx_ = static_cast<double>(yi - limits.lower) / step_size_;
            double innerSum = 0;

            for (int j = 0; j < step_size_; j++) {
                double xi = limits.lower + dx_ * j;
                double fx = integrand_(xi);                
                double innerArea = fx*dx_;
                innerSum += innerArea;
            }
            double outerArea = innerSum * dy_;
            outerSum += outerArea;
        }

        integral_ = outerSum;
        return integral_;
    }
};

This is my driver application:

main.cpp

#include <iostream>
#include <exception>
#include <cmath>

#include "Integrator.h"

constexpr double PI = 3.14159265358979;

constexpr double funcA(double x) {
    return x;
}

constexpr double funcB(double x) {
    return (x*x);
}

constexpr double funcC(double x) {
    return ((0.5*(x*x)) + (3*x) - (1/x));
}

double funcD(double x) {
    return sin(x);
}

int main() {
    //using namespace util; 
    try {

        std::cout << "Integration of f(x) = x from a=3.0 to b=5.0\nwith an expected output of 8\n";
        Integrator integratorA(Limits(3.0, 5.0), 10000, &funcA);
        std::cout << integratorA.evaluate() << '\n';        

        std::cout << "\n\nIntegration of f(x) = x^2 from a=2.0 to b=20.0\nwith an expected output of 2664\n";
        Integrator integratorB(Limits(2.0, 20.0), 10000, &funcB);
        std::cout << integratorB.evaluate() << '\n';

        std::cout << "\n\nIntegration of f(x) = (1\\2)x^2 + 3x - (1\\x) from a=1.0 to b=10.0\nwith an expected output of 312.6974\n";
        Integrator integratorC(Limits(1.0, 10.0), 10000, &funcC);
        std::cout << integratorC.evaluate() << '\n';

        std::cout << "\n\nIntegration of f(x) = sin(x) from a=0.0 to b=" <<PI<< "\nwith an expected output of 2\n";
        Integrator integratorD(Limits(0.0, PI), 10000, &funcD);
        std::cout << integratorD.evaluate() << '\n';


        std::cout << "\n\nTesting Double Integration of f(x) = (1\\2)x^2 + 3x - (1\\x) from [3,5] and [1,y]\nwith an expected output of 65.582\n";
        Integrator integratorE(Limits(3, 5), 500, &funcC);
        //double dy = integratorE.limits().upper - integratorE.limits().lower;
        integratorE.integrate(1);
        std::cout << integratorE.integral() << '\n';    

    } catch (const std::exception& e) {
        std::cerr << e.what() << std::endl;
        return EXIT_FAILURE;
    }

    return EXIT_SUCCESS;
}

And this is my output to the console when I run the program:

Integration of f(x) = x from a=3.0 to b=5.0
with an expected output of 8
7.9998


Integration of f(x) = x^2 from a=2.0 to b=20.0
with an expected output of 2664
2663.64


Integration of f(x) = (1\2)x^2 + 3x - (1\x) from a=1.0 to b=10.0
with an expected output of 312.6974
312.663


Integration of f(x) = sin(x) from a=0.0 to b=3.14159
with an expected output of 2
2


Testing Double Integration of f(x) = (1\2)x^2 + 3x - (1\x) from [3,5] and [1,y]
with an expected output of 65.582
65.3933

---

Here are my questions and concerns about the above code:

• What kind of improvements can be made to this code? I'm referring to "readability", "generically", and "portability".
  • I know that this isn't within a namespace as that is not primary concern within the context of this question. I can always put this in some defined namespace!
• Are there any apparent "code smells"?
• I have comments in my integrate function on not using the inner upper bounds...
  • How would I be able to incorporate the use of a defined inner upper bound?
• How can I extend my integrate function to perform even high orders of integration?
• Considering that the current implementation of performing a double integration has an O(n^2) complexity is there a way to reduce to this O(n) or O(log N)? If so, how?
• Are there any other optimizations that can be incorporated?
• Would the use of threads, multithreading, and parallel-programming be applicable here?
• Should I template this class?
• I'm also interested in any and all suggestions, tips, and feedback!

---

---

Extra useful information in regards to the design and implementation of my class

Its user-defined constructor requires three parameters/arguments in order to create an instance of an Integrator object.

• Its first requirement is the limits of integration that is defined by a simple Limits struct.
• Its second requirement is the step_size, normally the width of dx or the number of divisions in calculating the area of integration by parts.
• The third and final requirement is an std::function<double(double)> object.

About the Limits struct:

• It contains the lower and upper bounds of integration from a to b.
• It has a basic user-defined default constructor that takes the lower and upper bounds of integration as arguments. Both arguments can default to 0.
• It also contains an applyLimits(lower,upper) function. This simply acts as its constructor does with respect to its members by setting them or updating them.
• Access is purely public as these limits can be changed by the user at any given time. There is no restriction on the changing of the limits of integration.
• Both its constructor and its applyLimits() function will check if lower is greater than upper and if so it will swap them.

About the function-objects:

• They can be any of the following:
  • function object
  • function pointer
  • functor
  • lambda expression.
• These function objects can be defined as either constexpr or non-constexpr.
• Any are valid as long as they have the signature double(double) and can be stored in an std::function<> object.

About the construction and the use of the Integrator class object:

• What it can do
  • It can perform a definite integral of a single variable through the use of its evaluate() function.
  • It can also perform a second integral of the same integrand of a single variable through its function integrate(lower, upper).
  • It can also give you both the current dy and dx values, the integrand, and the current integral as well as the limits of integration.
• Construction
  • The limits or outer limits are defined when the object is instantiated through its user-defined constructor.
    • This is the default behavior for both single and double integrations.
  • The higher the step_size the more accurate the approximation.
    • Trade-offs: accuracy versus decrease in performance, time of execution taken.
  • The function object is stored as its integrand.
  • Versatility in being able to retrieve it back from the Integrator object and being able to use it at any time.
  • The inner limits are defined when calling its integrate() function.
    • The inner limits of integration are from [lower,y].
    • lower is passed into the function as an argument and
    • y is calculated on each iteration.
    • Currently within this context, upper is ignored for the inner limits and will default to 1>lower so that the Limits struct doesn't swap them.
  • See the note below in regards to expanding this class...

---

I'm considering expanding this to also allow for the user input of the inner upper limit of integration to be defined by the user and apply it within the algorithm to generate the appropriate values of integration. This has yet to be implemented and is one of my concerns. I would also like to incorporate an easy way to perform triple, quad, and quint integrations if performance bottlenecks can be reduced to a minimum while still being able to give an accurate approximation without a major decrease in performance. I would like to have the capabilities of an Integrator object to possibly accept another Integrator object as an argument.

Answer 1

You made multiple non-trivial edits while I wrote my answer, so there might be some divergence. (Personal annotation: Code should be (mostly) self explanatory. Don’t add a wall of text beforehand that words out what the code says anyway.)

Due to lack of expertise, I will not comment on possible mathematical improvements or multithreading.

Clear interface

I am bit confused by the Integrator class. The usage as shown in your main is as expected, but why are dx_, dy_ and integral_ member variables, which can be accessed, but do not contain any meaningful content (Or are even uninitialised for integral_!) until evaluate() or integrate() was called?

If this is meant to be some kind of result caching, then it should happen completely internally, maybe with an std::optional<double> integral_, which is set the first time something is calculated and then returned the next time. Also, both functions should not share the cached result. Since this is nothing but a wild guess, I’ll assume the smallest sensible interface as depicted by the main in the following.

struct Limits

In my opinion, applyLimits is completely redundant to the non-default constructor and introduces code duplication. It should be completely removed, since it can be replaced as follows:

some_limits.applyLimits(3., 4.);  //your function call
some_limits = Limits{3., 4.};     //shorter and clearer replacement

lower and upper should not be public (although you mention that this is intended) as lower <= upper is an invariant which cannot be guaranteed if the user meddles with the variables directly.

class Integrator

In the name of RAII, never have a constructor not initialize a member variable, in this case integral_!

As mentioned above, I will argue for a simplified interface here: Remove the member variables dx_, dy_ and integral_ as well as their respective getters completely and initialize them locally wherever needed. According to the rule of zero, do not explicitly default the destructor, as it is redundant and even deletes the move constructors!

Since your algorithm breaks for negative step_size_, use size_t instead of int as its type.

The loop over i in evaluate and the one over j in integrate are again code duplication; consider refactoring those. I would suggest something like this, with a reasonable name of course (omitted comments for brevity):

constexpr double evaluate() {
    double distance = limits_.upper - limits_.lower;
    auto dx = distance / step_size_;
                                                          
    return mathematically_descriptive_name(dx, limits_);
}

private:
double mathematically_descriptive_name(double dx, const Limits& limits) {
    double result = 0.;
    for (size_t i = 0; i < step_size_; ++i) {
        auto dy = integrand_(limits.lower + i * dx);
        auto area = dy * dx;
        result += area;
    }
    return result;
}

The loop in integrate can then be replaced with:

auto innerSum = mathematically_descriptive_name(dx, limits);

Whilst implementing this, I tripped over the fact that in integrate both the member variable limits_ as well as the local variable limits are used; you should make the names more distinguishable from each other to avoid confusion.

General style

Since you are using C++17, I would suggest a widespread use of [[nodiscard]]. Additionally, now that those additional member variables disappeared, all your functions can be const! With my interpretation of your interface, you could even make everything constexpr* and calculate everything at compile time – you would need to replace std::function though, maybe by templating the class over the function used.

Nowadays curly initialization, as already used by your constructors, is the way to go, e.g. use

Integrator integratorA{Limits{3.0, 5.0}, 10000, &funcA};

or even

auto integratorA = Integrator{Limits{3.0, 5.0}, 10000, &funcA};

for the main.

About templates

I would template both the struct and the class over a template<typename Field> instead of using double to increase flexibility of usage. Additionally, as mentioned earlier, for a constexpr evaluation you could consider using template<typename Func> and throwing lambdas in as parameters.

---

*std::swap is not constexpr before C++20. Until then, one could do a small trick to work around this like

constexpr Limits(double a = 0., double b = 0.) :
    lower{ a < b ? a : b },
    upper{ a < b ? b : a }
{}

Answer 2

You've implemented Riemann sums to numerically integrate functions. That's a good method if you may have very ugly/discontinuous functions and you don't care how long the integrals take. Plus it's simple and generally well understood. If the simplest choice is good enough for your application, then stick with it by all means.

However, there are other algorithms that will evaluate the integrand at fewer points and can handle definite integrals with infinite bounds.

I'm not going to dive into the alternative methods here, but I'll point you to two resources that explain the methods better than I can:

• https://en.wikipedia.org/wiki/Numerical_integration#Methods_for_one-dimensional_integrals. This is a really good article. I think the pictures show how you can get a more accurate integral with fewer evaluations.

• https://www.boost.org/doc/libs/1_73_0/libs/math/doc/html/math_toolkit/gauss.html. boost::math::quadrature is Boost's version of your code. You may enjoy reading the docs and/or the source code to see how they implement more performant algorithms. As a general rule, whenever you implement something general in C++ it's worth checking if one of the major C++ general purpose libraries has a version.

---

Integrator(..., int stepSize, ...)

stepSize is only useful in some integration algorithms. IMO that implies this argument is a leaky abstraction. Also, why should this be an int?

I think what you really want is a way to control the precision of your answer. Maybe a double maximumError argument could achieve that?

---

Why is Integrator a class rather than a function?

---

Typically, integral(from: a, to:b) == -integral(from:b, to:a) (https://en.wikipedia.org/wiki/Integral#Conventions). In your implementation, they are equivalent.

---

definite double integration of a single variable

This confused me because you actually introduce a second variable in the limits of integration of the inner integral. Also you have some little bugs in the integrate function which I think you would have easily caught if you added more test cases.

Imagine your single definite integral functions had the signature template<typename F> double integrate1x(double lower, double upper, F f). Then you could implement your outer integral with the same method:

// \int_(L)^(U) \int_(g(y))^(h(y)) f(x) dx dy
template <typename F, G, H>
double integrate2x(double L, double U, G g, H h, F f) {
    return integrate1x(L, U, [&](double y) {
        return integrate1x(g(y), h(y), f);
    });
}