8
\$\begingroup\$

This is a follow up of A simple definite integrator class of a single variable in C++.

I took most of the advice from Emma X and some from sudo rm -rf slash.


Here is my fully revised code:

Integrator.h

#pragma once

#include <type_traits>

template <typename Field>
struct Limits {
    Field lower;
    Field upper;

    constexpr Limits(Field a = 0, Field b = 0) : 
        lower{ a }, 
        upper{ b }
    {}
};

template <typename LimitType, typename Func>
class Integrator {
private:
    Limits<LimitType> limits_;
    size_t step_size_;
    Func integrand_;

public:
    Integrator(Limits<LimitType> limits, size_t stepSize, Func integrand) :
        limits_{ limits },
        step_size_{ stepSize },
        integrand_{ integrand }
    {}

    Limits<LimitType> limits() const { return limits_; }
    Func* integrand() { return &integrand_; }

    // This is always a 1st order integration!   
    constexpr auto evaluate() {
        if (limits_.lower == limits_.upper) {
            return 0.0;
        } else if (limits_.lower > limits_.upper) {
            limits_ = Limits{ limits_.upper, limits_.lower };
            auto distance = limits_.upper - limits_.lower;
            auto dx = distance / step_size_;
            return -calculate(dx);
        } else {
            auto distance = limits_.upper - limits_.lower;
            auto dx = distance / step_size_;
            return calculate(dx);
        }
    }

    // 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...
    template<typename ValueType>
    constexpr ValueType doubleIntegral(Limits<LimitType> innerLimits) {
        // Currently the inner upper bound is not being used, if it is <= lower bound,
        // it will default to lower bound + 1, just to ensure that an integration is
        // is calculated. Eventually this may be expanded to return true values
        // such as 0 when the lower and upper bounds are equal or the (-)integration if a > b.
        if (innerLimits.upper <= innerLimits.lower) innerLimits.upper = innerLimits.lower + static_cast<LimitType>(1);

        ValueType sum = 0;
        ValueType dy = (ValueType)(limits_.upper - limits_.lower) / step_size_;                   
        for (size_t i = 0; i < step_size_; i++) {
            ValueType yi = (ValueType)(limits_.lower + i * dy);
            ValueType dx = (yi - innerLimits.lower) / step_size_;                    
            Integrator inner{ innerLimits, step_size_, integrand_ };
            sum += inner.calculate(dx) * dy;    
        }    
        return sum;
    }

private:
    template<typename ValueType>
    constexpr auto calculate(ValueType dx) ->
    std::enable_if_t<std::is_invocable_v<Func&, ValueType&>, ValueType> {
        ValueType result = 0.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;
    }
};  

main.cpp

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

#include "Integrator.h"

constexpr double PI = 3.14159265358979;

double funcA0(double x) {
    return x;
}

template <typename T>
constexpr T funcA(T x) {
    return x;
}

template <typename T>
constexpr T funcB(T x) {
    return x * x;
}

template <typename T>
T funcC(T x) {
    return sin(x);
}

template <typename T>
constexpr T funcD(T x) {
    return (0.5*(x*x) + (3 * x) - (1 / x));
}

double funcE(double x) {
    return cos(x);
}

int main() {
    try {

        std::cout << "Integration of f(x) = x from a=3.0 to b=3.0 with an expected output of 0\n";
        std::cout << "Testing non template function object\n";
        Integrator integratorA0A{ Limits{3.0, 3.0}, 10000, funcA0 };
        std::cout << integratorA0A.evaluate() << '\n';

        std::cout << "Testing non template function object by reference\n";
        Integrator integratorA0B{ Limits{3.0, 3.0}, 10000, &funcA0 };
        std::cout << integratorA0B.evaluate() << '\n';

        std::cout << "Testing function template object\n";
        Integrator integratorA0C{ Limits{3.0, 3.0}, 10000, funcA<double> };
        std::cout << integratorA0C.evaluate() << '\n';

        std::cout << "Testing function template object by reference\n";
        Integrator integratorA0D{ Limits{3.0, 3.0}, 10000, &funcA<double> };
        std::cout << integratorA0D.evaluate() << '\n';


        std::cout << "\nIntegration of f(x) = x from a=3.0 to b=5.0 with an expected output of 8\n";
        std::cout << "Testing non template function object\n";
        Integrator integratorA1A{ Limits{3.0, 5.0}, 10000, funcA0 };
        std::cout << integratorA1A.evaluate() << '\n';

        std::cout << "Testing non template function object by reference\n";
        Integrator integratorA1B{ Limits{3.0, 5.0}, 10000, &funcA0 };
        std::cout << integratorA1B.evaluate() << '\n';

        std::cout << "Testing function template object\n";
        Integrator integratorA1C{ Limits{3.0, 5.0}, 10000, funcA<double> };
        std::cout << integratorA1C.evaluate() << '\n';

        std::cout << "Testing function template object by reference\n";
        Integrator integratorA1D{ Limits{3.0, 5.0}, 10000, &funcA<double> };
        std::cout << integratorA1D.evaluate() << '\n';


        std::cout << "\nIntegration of f(x) = x from a=5.0 to b=3.0 with an expected output of -8\n";
        std::cout << "Testing non template function object\n";
        Integrator integratorA2A{ Limits{5.0, 3.0}, 10000, funcA0 };
        std::cout << integratorA2A.evaluate() << '\n';

        std::cout << "Testing non template function object by reference\n";
        Integrator integratorA2B{ Limits{5.0, 3.0}, 10000, &funcA0 };
        std::cout << integratorA2B.evaluate() << '\n';

        std::cout << "Testing function template object\n";
        Integrator integratorA2C{ Limits{5.0, 3.0}, 10000, funcA<double> };
        std::cout << integratorA2C.evaluate() << '\n';

        std::cout << "Testing function template object by reference\n";
        Integrator integratorA2D{ Limits{5.0, 3.0}, 10000, &funcA<double> };
        std::cout << integratorA2D.evaluate() << '\n';


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

        std::cout << "\n\nIntegration of f(x) = sin(x) from a=0.0 to b=" << PI << " with an expected output of 2\n";
        Integrator integratorC{ Limits{0.0, PI}, 10000, &funcC<double> };
        std::cout << integratorC.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 integratorD{ Limits{1.0, 10.0 }, 10000, &funcD<double> };
        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 integrator2D{ Limits{3, 5}, 10000, &funcD<double> };
        std::cout << integrator2D.doubleIntegral<double>(1.0) << '\n';

        std::cout << "\n\nTesting Double Integration of f(x) = cos(x) from [0," << PI / 2 << "] and [0.25,y]\nwith an expected output of 0.61137\n";
        Integrator integratorE{ Limits{0.0, PI / 2}, 10000, &funcE };
        std::cout << integratorE.doubleIntegral<double>(0.25) << '\n';

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

    return EXIT_SUCCESS;
}

Output

Integration of f(x) = x from a=3.0 to b=3.0 with an expected output of 0
Testing non template function object
0
Testing non template function object by reference
0
Testing function template object
0
Testing function template object by reference
0

Integration of f(x) = x from a=3.0 to b=5.0 with an expected output of 8
Testing non template function object
7.9998
Testing non template function object by reference
7.9998
Testing function template object
7.9998
Testing function template object by reference
7.9998

Integration of f(x) = x from a=5.0 to b=3.0 with an expected output of -8
Testing non template function object
-7.9998
Testing non template function object by reference
-7.9998
Testing function template object
-7.9998
Testing function template object by reference
-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) = sin(x) from a=0.0 to b=3.14159 with an expected output of 2
2


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


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


Testing Double Integration of f(x) = cos(x) from [0,1.5708] and [0.25,y]
with an expected output of 0.61137
0.611325

List of Fixed Issues:

  • I removed the default destructor.
  • I removed the unnecessary member variables and their equivalent functions.
    • I'm now caching them on the stack.
  • I removed the member function from the Limits struct.
  • I also removed the swap from the limits.
    • if a==b then the integral is 0 and
    • if a>b then the integral is (-) the integral from b to a.
    • These checks for the limits of integration I moved into the Integrator class.
  • I tried to make sure everything follows RAII.
  • The class now works in a constexpr context provided the function objects are constexpr.
  • I moved the loop calculations for the approximation from the evaluate function into its own private member function.
  • I'm using {} initialization for all object instantiations.
  • I templated both classes.
  • I can also pass in a function object that is also templated.

All of the changes to my classes are all based on the advice of Emma X. I did take some of the advice from sudo rm -rf slash in that I read up on the links that they had provided. As for changing the structure of the code from a class to be just a set of functions, I did not. I have my own reasons why I want this to be a class object. I plan on expanding this later for doing different types of integrations and I will need this class to keep track of its integrand and its limits of integration.

Yes, I do understand that the current implementation is based on Riemann Sums. This is the intentional base structure of my class since it is the simplest form in approximate integrals. Eventually, I plan on expanding this to incorporate different methods and algorithms. These will be in the form of member functions naming them according to their algorithm.


Other changes that were made:

  • I renamed integrate to doubleIntegral.
  • I also changed its signature to take a Limits struct instead of 2 individual parameters.
  • I also removed the local innerLimits variable since it is now being passed in, and just update it accordingly if needed.
  • In both the evaluate and doubleIntegral function, I have early returns if the limits are equal.
  • I only have a check to see if a > b for the 1st order integration to return the (-)integral. I did not apply this behavior to the doubleIntegral function yet.
  • I removed a lot of unnecessary local variables in my functions and inlined many of the calculations.
  • Finally I added a sanity check to make sure that the template argument that is passed into the construct is an invokable object.
  • The last of the changes, I modified my driver program slightly to show different cases of the same function being used in different contexts and to show the results when the limits are equal, less than and greater than... and I added a second example of a double integration using a trig function.

Where do I go from here?

  • I'd like to incorporate the innerLimit's upper into the calculations when calculating doubleIntegral.
  • I'd like to know if there is any room for improvements, efficiency, readability, portability, etc.
  • Are there still any existing code smells?
    • -Note-, mentioned in the previous post: I'm not concerned with namespace; that can be done at any time!
  • Emma X had mentioned the use of [[nodiscard]]; however, I'm not familiar with its use... I've never really used that feature of C++ before. It is new to me.
  • I have not yet tried to use lambdas within this code yet; that is on my to-do list.
  • I would also like to be able to add in the ability or the concept to incorporate infinity as the limits of integration and since this is a numerical type integrator, only the upper or lower bound can go to ±∞, but not both. I'm not at the stage of returning back a family of functions with a constant of integration; that's a little beyond my skills at the moment since that would involved a function lookup table and an extensive parser.

Let me know what you think. I look forward to your advice, constructive criticism, opinions, positions, and all feedback is welcomed!

\$\endgroup\$
3
  • 1
    \$\begingroup\$ How do you argue that initializing the loop variable before the for-loop helps with optimization? It compiles to the same assembly, but your way is worse in terms of readability and scoping. \$\endgroup\$
    – Emma X
    Commented Jun 3, 2020 at 21:26
  • \$\begingroup\$ @EmmaX Well, it may not make a difference these days, compilers have improved quite a bit, but I only have to declare and initialize it once outside the loop, then just iterate through the loop. This is an old trick going back to when I first started learning the C/C++ languages... These days, the compilers will probably optimize it anyway... Then again this trick is commonly used for extreme efficiency when working with "nested loops" and it usually pertains to the inner loop's counter. \$\endgroup\$ Commented Jun 3, 2020 at 23:17
  • 2
    \$\begingroup\$ @EmmaX After thinking about it for a little while, yeah it's not necessary to have the counters outside of the for loops... If there were to be any kind of gain it would be minimal at best and the tradeoff for readability simply outweighs it. So I reverted my code to put them back in their respective for loops and I removed the highlighted bullet that mentions them. \$\endgroup\$ Commented Jun 4, 2020 at 1:20

3 Answers 3

6
\$\begingroup\$

#pragma once is non-standard.


size_t is never defined. If it's a misspelling of std::size_t, we should be including one of the headers that defines it (rather than relying on an artefact of our compiler both to transitively include it and to declare a standard-namespace equivalent, neither of which is portable).


Implicit conversion is generally a bad idea, and especially so for template types. Make the Limits constructor explicit:

    explicit constexpr Limits(Field a = 0, Field b = 0) :
        lower{ std::move(a) },
        upper{ std::move(b) }
    {}

We'll need to explicitly create a Limits when we call doubleIntegral() but I consider that a Good Thing. (In fact, why not just accept a LimitType scalar there, since we never use the upper member?)


If we could use C++20, then we could constrain Func template argument:

#include <concepts>
template <typename LimitType, std::regular_invocable<LimitType> Func>

private: is redundant at the start of a class definition.


step_size is a misleading name: it's a step count. We should ensure it's not zero if we want to avoid infinities.

Using std::size_t is probably wrong here, given that it's commonly too large to convert to double without loss of precision for larger values.


limits() and integrand() should both be const, and return reference-to-const. Don't return by value when a copy isn't needed, and don't return a pointer that can't be null.


No need for else in evaluate(), where the if statement returns. The last two branches can be combined if we simply normalise limits:

        if (limits_.lower == limits_.upper) {
            return 0.0;
        }

        if (limits_.lower > limits_.upper) {
            std::swap(limits_.lower, limits_.upper);
        }
        auto distance = limits_.upper - limits_.lower;
        auto dx = distance / step_size_;
        return calculate(dx);

I'd actually move the normalisation into the constructor, since nothing changes limits subsequently.

It's not clear that 0.0 is an empty result of the correct type. We ought to be returning a default-constructed std::invoke_result_t<Func, LimitType>.

The division producing dx may be inexact - it's better to use the exact difference between adjacent values.


I don't think doubleIntegral() is helpful. It's too inflexible, compared to the obvious approach of simply composing integration. The latter lets us integrate functions of two (or more) variables, for example.


calculate() can be const.

Instead of evaluating the function at one side of each strip, it's better to do so at the midpoint. Or evaluate at both sides and compute the area of a trapezium rather than a rectangle.

Further possibility - adaptive ẟx so that it's larger when f″(x) is larger, and smaller when f(x) is close to linear.


I see no need for a class here at all - a simple function would be easier to use and more flexible:

#include <concepts>
#include <type_traits>

template<typename Func, std::floating_point Arithmetic,
         std::unsigned_integral Count = unsigned short>
requires std::regular_invocable<Func, Arithmetic>
constexpr auto integrate(Func&& func,
                         Arithmetic from, Arithmetic to,
                         Count step_count = 10000u)
{
    std::invoke_result_t<Func, Arithmetic> result = {};
    if (to == from) {
        return result;
    }
    if (step_count == 0) {
        step_count = 1;
    }
    // sum trapezia bounded by x0, x1
    auto x0 = from;
    auto y0 = func(from);
    for (auto i = 1u;  i <= step_count;  ++i) {
        auto x1 = ((step_count - i) * from + i * to) / step_count;
        auto y1 = func(x1);
        result += (y0 + y1) * (x1 - x0) / 2;
        x0 = x1; y0 = y1;       // this iteration's end is next one's start
    }
    return result;
}
  • For thoroughness, we might include <functional> and replace func(x1) with std::invoke(func, x1)).
  • Instead of the counted loop, in C++20 we might iterate over a transform of std::views::iota to get our x1 values.

Usage changes just slightly:

static constexpr double π = 3.14159265358979;

std::cout << "Integration of f(x) = x from a=3.0 to b=3.0 with an expected output of 0\n";
std::cout << integrate(funcA0, 3.0, 3.0) << '\n';

std::cout << "\nIntegration of f(x) = x from a=3.0 to b=5.0 with an expected output of 8\n";
std::cout << integrate(funcA0, 3.0, 5.0) << '\n';

std::cout << "\nIntegration of f(x) = x from a=5.0 to b=3.0 with an expected output of -8\n";
std::cout << integrate(funcA0, 5.0, 3.0) << '\n';


std::cout << "\n\nIntegration of f(x) = x² from a=2.0 to b=20.0 with an expected output of 2664\n";
std::cout << integrate(funcB<double>, 2.0, 20.0) << '\n';

        std::cout << "\n\nIntegration of f(x) = sin(x) from a=0.0 to b=π with an expected output of 2\n";
std::cout << integrate<double(*)(double)>(std::sin, 0.0, π) << '\n';

std::cout << "\n\nIntegration of f(x) = ½x² + 3x - x⁻¹ from a=1.0 to b=10.0 \nwith an expected output of 312.6974\n";
std::cout << integrate(funcD<double>, 1.0, 10.0) << '\n';

std::cout << "\n\nTesting Double Integration of f(x) = ½x² + 3x - x⁻¹ from [3,5] and [1,y]\nwith an expected output of 65.582\n";
std::cout << integrate([](double y){ return integrate(funcD<double>, 1.0, y, 10000u); },
                       3.0, 5.0) << '\n';

std::cout << "\n\nTesting Double Integration of f(x) = cos(x) from [0,π/2] and [0.25,y]\nwith an expected output of 0.61137\n";
std::cout << integrate([](double y){ return integrate(funcE, 0.25, y, 10000u); },
                       0.0, π/2) << '\n';

In the test code, we've misspelt std::sin and std::cos. Also, prefer not to use all-caps for C++ identifiers - we normally use that convention for preprocessor macros (which don't obey C++ rules of scope, for instance).


Finally, note that there can be numerical stability problems with particular functions that are hard to eliminate because the order of summation can produce catastrophic subtraction or dissimilar-magnitude addition. I offer no immediate solution to this problem.

\$\endgroup\$
1
  • 3
    \$\begingroup\$ Great answer. I'd also like to stress that this should not have been a class at all; most of the problems come from that fact alone. About the double integral, it was very weird to begin with. From the name I would have expected it to be able to integrate a function of two variables, i.e. calculating the volume under its surface. \$\endgroup\$
    – G. Sliepen
    Commented Feb 3 at 15:15
5
\$\begingroup\$

Use a better PI

Typical double can require 17 significant decimal digits in source code to encode the best double. Even this "17" may be more on select systems.

OP used 15 and 17 would have resulted in a better value of machine_pi. OP's PI is 7 ULP away from the best machine_pi.

//                       1 23456789 12345 
// constexpr double PI = 3.14159265358979;
constexpr double PI = 3.1415926535897932;

There is nearly no downside to using even more digits when expressing a constant in source code.

  constexpr double PI_15 = 3.14159265358979;
  constexpr double PI_17 = 3.1415926535897932;
  constexpr double PI_Many = 3.1415926535897932384626433832795;
  
  3.14159265358979_000737349451810587197542190551757812500000000
  3.14159265358979_311599796346854418516159057617187500000000000
  3.14159265358979_311599796346854418516159057617187500000000000
  3.14159265358979_323846264338327950288419716939937510582097494... π

  0x3.243f6a8885a22
  0x3.243f6a8885a30
  0x3.243f6a8885a30
  0x3.243f6a8885a308d... π

As I posted another answer concerning an entirely different review aspect, I'll make this one wiki.

\$\endgroup\$
2
  • 1
    \$\begingroup\$ Or since C++20: std::numbers::pi \$\endgroup\$
    – G. Sliepen
    Commented Feb 22 at 9:47
  • 1
    \$\begingroup\$ @G.Sliepen Yes, std::numbers::pi works for this specific case. Still, in general, FP constants should be coded with sufficient precision. \$\endgroup\$ Commented Feb 22 at 16:45
4
\$\begingroup\$

Correct algorithm??

Below looks wrong.

    for (size_t i = 0; i < step_size_; ++i) {
        auto dy = integrand_(limits_.lower + i * dx);
        auto area = dy * dx;
        result += area;
    }

For simplicity, assume we are integrating from 0.0.

That calculates the sum of dx*(y(0) + y(dx) + y(dx*2) ... y(dx*(step_size-1)))

I'd expect integration to follow the Trapezoidal rule. That includes both ends (weighed half of the mid-samples).

Sum of dx*(y(limits_.lower)/2 + y(dx) + y(dx*2) ... y(dx*(step_size-1.0)) + y(limits_.upper)/2).


Minor: efficiency and reduced errors

We could perform the * dx once, after the loop.

This improves efficiency and slight reduced error generation. It may negatively (or positively) affect usable computational range.


Step size

Using a power-of-2 for step_size has computational advantages as it reduces division errors.


I'd expect a better computational result with these changes.

Maybe something like

        ValueType result = 
            (integrand_(limits_.lower) + integrand_(limits_.upper))/2;
        for (size_t i = 1; i < step_size_; ++i) {
            result += integrand_(limits_.lower + i * dx);
        }
        result *= dx;
\$\endgroup\$
4
  • \$\begingroup\$ Yes. You're saying that 10000 is Too Big, it covers up tiny sins. Better to magnify them by supplying a value closer to 10. \$\endgroup\$
    – J_H
    Commented Feb 23 at 17:12
  • \$\begingroup\$ @J_H Not quite. Yet using power-of-2 for step_size has computational advantages. \$\endgroup\$ Commented Feb 23 at 17:29
  • \$\begingroup\$ Oh, sorry, I meant "use a small value during testing" so we notice issues. A production call would of course use a suitably large value. // Also, mathematic code is usually improved if OP cites their references. \$\endgroup\$
    – J_H
    Commented Feb 23 at 17:31
  • \$\begingroup\$ @J_H Yes, I see your point about small values now. \$\endgroup\$ Commented Feb 23 at 17:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.