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

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.

• Please do not update the code in your question to incorporate feedback from answers, doing so goes against the Question + Answer style of Code Review. This is not a forum where you should keep the most updated version in your question. Please see what you may and may not do after receiving answers.
– Mast
May 30, 2020 at 11:46
• @Mast, you're right, it's been a while since I've posted here! May 30, 2020 at 11:48
• @Mast Should I move the extra information of the class below the questions and concerns? Having the class, the driver application and its results before towards the beginning of the post? May 30, 2020 at 11:49
• Now answers have arrived, the preferred choice is for you accept the current answer if it was useful to you, to wait 24h between posting questions and post a follow-up linking back to this question. Your new question can contain all the extra information (which probably should've been part of this question in the first revision). Do note that you raise a lot of questions in your question that are actually feature requests. You can note them as points of interest for a next iteration, a goal you're working towards, but they can not be the goal of the question itself.
– Mast
May 30, 2020 at 11:54
• Here's a follow-up question: codereview.stackexchange.com/questions/243339/… Jun 3, 2020 at 19:35

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 unitialized 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 whereever needed. According to the rule of zero, do not explicitely 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 that. I would suggest something like this, with a reasonable name of course (ommited 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.

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


or even

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


for the main.

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 }
{}

• I moved the structured information below the class and the questions and concerns as a reference to the thought in the design and implementation of my code. May 30, 2020 at 11:12
• @FrancisCugler We take answer invalidation quite seriously here, please do not add to or modify the code in the question after the first answer arrives.
– Mast
May 30, 2020 at 11:47
• @FrancisCugler That's unfortunate, but that risk is yours. To avoid this risk, make sure the first revision of your post is complete and correct.
– Mast
May 30, 2020 at 11:57
• @FrancisCugler And that's why I haven't rolled it back :-) But any edit carries a risk and that risk is yours the moment answers have come in. That's all I'm saying.
– Mast
May 30, 2020 at 12:01
• @FrancisCugler I made an edit about constexpr-ness of std::swap. About having other Instances of the class use dx and dy: Consider calculating those values on the fly since it isn’t expensive. Having getters to values which might be set to "meaningless" default values is bad style. Independent of how you choose to do this, it should still be private, as it is an implementation detail. Instances of the same class can access private members, see this simple example on godbolt. May 30, 2020 at 12:41

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:

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