Numerical integration in C++: Newton-Cotes formulas

Question

I have tried my hand with implementing simple quadrature formulas in C++.

Definite integral: $$\int_a^b f(x) dx$$ Domain of integration \$[a, b]\$ divided into \$n\$ intervals of equal length \$h = (b - a) / n\$.

Midpoint of interval \$h_i = b_i - a_i\$ is $$x_i = a + h / 2 + (i - 1) h$$ with \$i = 1, \dots, n \$.

Newton-Cotes quadrature formulas implemented in the code:

• Midpoint rule (constant approx) $$F_M = h \sum_{i=1}^n f(x_i)$$
• Trapezoidal rule (linear approx) $$F_T = \frac{h}{2} \left(f(a) + f(b) \right) + h \sum_{i=1}^{n-1} f(a_i)$$
• Simpson's rule (quadratic approx) $$F_S = \frac{h}{6} \left(f(a) + f(b) \right) + \frac{h}{3} \sum_{i=1}^{n-1} f(a_i) + 2 \frac{h}{3} \sum_{i=1}^{n} f(x_i) $$

I decided to go with a parametrized factory pattern that has no subclasses involved.

I employ \$ f(x) = \sqrt{x} \exp{(-x)} \$ as test function, to be integrated between \$a = 1 \$ and \$b = 3\$.

I wonder whether the #include "NewtonCotesFormulas.cpp" at the end of the NewtonCotesFormulas.h header file could be somehow avoided. I used it because without I'd get linking errors.

Here is the code.

main.cpp:

#include <iostream>
#include <cmath>
#include "NewtonCotesFormulas.h"
#include "NewtonCotesFactory.h"

using namespace std;

double func(double x) {
    return sqrt(x) * exp(-x);
}

int main()
{
    cout << "Choose method:\n"
         << "Midpoint numeric quadrature -----> (1)\n"
         << "Trapezoidal numeric quadrature --> (2)\n"
         << "Simpson's numeric quadrature ----> (3)\n" << endl;
// //user input code:
//    int choice;
//    cin >> choice;
//    while (choice < 1 || choice > 3) {
//        cout << "Pick again a number from 1 to 3\n";
//        cin >> choice;
//    }
    NewtonCotesFactory creator;
// alternative to user input, print all results at once:
    for (int choice : {1,2,3}) {
        NewtonCotesFormulas* rule = creator.createQuadrature(choice, 1, 3, 0.000001);
        double result = rule->printConvergenceValues(4, &func);
        cout << "Final result within tolerance: " << result << endl << endl;
        delete rule;
    }

    return 0;
}

NewtonCotesFormulas.h:

#ifndef NEWTONCOTESFORMULAS_H_INCLUDED
#define NEWTONCOTESFORMULAS_H_INCLUDED

#include <string>

class NewtonCotesFormulas {
public:
    NewtonCotesFormulas(double, double, double, std::string);
    virtual ~NewtonCotesFormulas();
    virtual double computeIntegral(int, double (*func)(double)) = 0;
    double convergeToTol(int, double (*func)(double));
    double printConvergenceValues(int, double (*func)(double));

protected:
    double inf, sup;
    double tol;
    std::string method_name;
};

class Midpoint: public NewtonCotesFormulas {
public:
    Midpoint(double, double, double);
    ~Midpoint(){};
    double computeIntegral(int, double (*func)(double)) override;
};

class Trapezoidal: public NewtonCotesFormulas {
public:
    Trapezoidal(double, double, double);
    ~Trapezoidal(){};
    double computeIntegral(int, double (*func)(double)) override;
};

class Simpsons: public NewtonCotesFormulas {
public:
    Simpsons(double, double, double);
    ~Simpsons(){};
    double computeIntegral(int, double (*func)(double)) override;
};

#include "NewtonCotesFormulas.cpp"
#endif // NEWTONCOTESFORMULAS_H_INCLUDED

NewtonCotesFormulas.cpp:

#include <iomanip>
#include "NewtonCotesFormulas.h"

using namespace std;

NewtonCotesFormulas::NewtonCotesFormulas(double a, double b, double inp_tol, string name):
    inf(a), sup(b), tol(inp_tol), method_name(name) {}

NewtonCotesFormulas::~NewtonCotesFormulas() {}

double NewtonCotesFormulas::convergeToTol(int intervals, double (*func)(double)) {
    double old_integral = this->computeIntegral(intervals, func);
    intervals *= 2;
    double new_integral = this->computeIntegral(intervals, func);
    double difference = abs(new_integral - old_integral);

    while(difference > tol) {
        cout << setprecision(8) << string(9 - to_string(intervals).size(), ' ') <<  // blank padding used to format output table
                intervals << string(3,' ') << fixed << new_integral << string(3,' ') << scientific << difference << endl;
        old_integral = new_integral;
        intervals *= 2;
        new_integral = this->computeIntegral(intervals, func);
        difference = abs(new_integral - old_integral);
    }
    cout << setprecision(8) << fixed << string(9 - to_string(intervals).size(), ' ') <<
            intervals << string(3,' ') << new_integral << string(3,' ') << scientific << difference << endl;

    return new_integral;
}

double NewtonCotesFormulas::printConvergenceValues(int intervals, double (*func)(double)) {
    cout << method_name << " method, with tolerance " << tol << "\n";
    cout << "intervals" << string(5,' ') << "integral " << string(9,' ') << "tol\n";
    cout << setprecision(8) << string(9 - to_string(intervals).size(), ' ') <<
            intervals << string(3,' ') << this->computeIntegral(intervals, func) <<
            string(3,' ') << fixed << 0.0 << endl;
    double result = this->convergeToTol(intervals, func);

    return result;
}

Midpoint::Midpoint(double a, double b, double inp_tol):
    NewtonCotesFormulas(a, b, inp_tol, "Midpoint"){}

Trapezoidal::Trapezoidal(double a, double b, double inp_tol):
    NewtonCotesFormulas(a, b, inp_tol, "Trapezoidal"){}

Simpsons::Simpsons(double a, double b, double inp_tol):
    NewtonCotesFormulas(a, b, inp_tol, "Simpson's"){}

double Midpoint::computeIntegral(int intervals, double (*func)(double)) {
    double interval_width = (sup - inf) / intervals;
    double result = 0;
    for (int i = 1; i <= intervals; ++i)
        result += func(inf + (i - 0.5) * interval_width);

    return result * interval_width;
}

double Trapezoidal::computeIntegral(int intervals, double (*func)(double)) {
    double interval_width = (sup - inf) / intervals;
    double result = (func(inf) + func(sup)) / 2.;
    for (int i = 1; i <= intervals - 1; ++i)
        result += func(inf + i * interval_width);

    return result * interval_width;
}

double Simpsons::computeIntegral(int intervals, double (*func)(double)) {
    double interval_width = (sup - inf) / intervals;
    double result = (func(inf) + func(sup)) / 6;
    for (int i = 1; i <= intervals - 1; ++i)
        result += func(inf + i * interval_width) / 3 + 2 * func(inf + (i - 0.5) * interval_width) / 3;
    result += 2 * func(inf + (intervals - 0.5) * interval_width) / 3;

    return result * interval_width;
}

NewtonCotesFactory.h:

#ifndef NEWTONCOTESFACTORY_H_INCLUDED
#define NEWTONCOTESFACTORY_H_INCLUDED

#include <string>
#include <stdexcept>
#include "NewtonCotesFormulas.h"

class NewtonCotesFactory {      // Parametrized factory, no children
public:
    virtual NewtonCotesFormulas* createQuadrature(int choice, double a, double b, double inp_tol) const {
        if (choice == 1)
            return new Midpoint(a, b, inp_tol);
        else if (choice == 2)
            return new Trapezoidal(a, b, inp_tol);
        else if (choice == 3)
            return new Simpsons(a, b, inp_tol);
        else
            throw std::runtime_error("Pick either choice = 1, 2, or 3");

        return 0;  // never gets here
    }
    virtual ~NewtonCotesFactory(){}
};

#endif // NEWTONCOTESFACTORY_H_INCLUDED

Results look fine:

Choose method:
Midpoint numeric quadrature -----> (1)
Trapezoidal numeric quadrature --> (2)
Simpson's numeric quadrature ----> (3)

Midpoint method, with tolerance 1e-06
intervals     integral          tol
        4   0.40715731   0.00000000
        8   0.40807542   9.18106750e-04
       16   0.40829709   2.21674991e-04
       32   0.40835199   5.49009778e-05
       64   0.40836569   1.36924160e-05
      128   0.40836911   3.42104469e-06
      256   0.40836996   8.55132348e-07
Final result within tolerance: 4.08369962e-01

Trapezoidal method, with tolerance 1.00000000e-06
intervals     integral          tol
        4   4.10757439e-01   0.00000000
        8   0.40895737   1.80006398e-03
       16   0.40851640   4.40978614e-04
       32   0.40840674   1.09651812e-04
       64   0.40837937   2.73754169e-05
      128   0.40837253   6.84150045e-06
      256   0.40837082   1.71022788e-06
      512   0.40837039   4.27547766e-07
Final result within tolerance: 4.08370390e-01

Simpson's method, with tolerance 1.00000000e-06
intervals     integral          tol
        4   4.08357353e-01   0.00000000
        8   0.40836940   1.20498403e-05
       16   0.40837019   7.90455982e-07
Final result within tolerance: 4.08370194e-01

Process returned 0 (0x0)   execution time : 0.063 s
Press any key to continue.

Answer 1

I wonder whether the #include "NewtonCotesFormulas.cpp" at the end of the NewtonCotesFormulas.h header file could be somehow avoided. I used it because without I'd get linking errors.

Yea, that's simply wrong. You need to give both that CPP file and the main.cpp file to the compiler, so it knows that both are part of the program to build.

---

Don’t write using namespace std;.

You can, however, in a CPP file (not H file) or inside a function put individual using std::string; etc. (See SF.7.)

---

⧺SL.io.50 Don't use endl.

---

NewtonCotesFormulas* rule = creator.createQuadrature(choice, 1, 3, 0.000001);
    ⋮
delete rule;

⧺C.149 — no naked new or delete.

You can simply have createQuadrature return a unique_ptr<NewtonCotesFormulas> instead. Use auto when you call it.

---

class NewtonCotesFactory {      // Parametrized factory, no children
public:
    virtual NewtonCotesFormulas* createQuadrature(int choice, double a, double b, double inp_tol) const {

Not only does it not have any derived classes, the comment even says "no children". So what is the purpose of using virtual on the function? There is no member data at all, so it can be a static member. Actually, why is it even a class? You just need a single stand-alone free function.

As for the class as-written, never write a destructor or other special member function with an empty body. That will stop the compiler from realizing that it is trivial or has the built-in meaning with no changes, and causes it to stop auto-generating some of the other members because now it's "custom". Use =default instead.

---

NewtonCotesFormulas::~NewtonCotesFormulas() {}
Just don't define this. The only reason for having it is to make it virtual, but you could have put the definition directly in-line in the header. But as explained above, don't do that for an empty destructor; just write it as =default.

NewtonCotesFormulas(double, double, double, std::string); Don't pass the string by value ! (Though, for a constructor with a "sink" parameter, you could use the "sink" idiom but I don't think you were intending to use this. Generally, pass strings as const ... & but it's better to use a string_view for parameters.)

---

double old_integral = this->computeIntegral(intervals, func);
Don't write this-> generally to access members. They are in scope. That is not the C++ way.

Use const whenever you can.

---

good!

I see you used override on the derived classes.

---

not "wrong" but...

You're limiting the supplied function to being a plain non-member (or static member) function, rather than any kind of callable thing like a lambda.

And, you are preventing the implementation from inlining the call to the supplied func (in two different ways), which can have a serious effect on performance especially when the function is so simple.

You might want to look into doing this part better, like how the standard library handles the function passed to sort and everywhere that comparison or callback functions are passed, in a future study session.

I think the most important thing to master next would be how to work with projects that have more than one CPP file.

Keep it up!

Answer 2

Use an enum to give names to choices

Consider using an enum, or even better an enum class, to enumerate the possible choices for the Newton-Cotes formulas:

enum class Formula {
    MIDPOINT = 1,
    TRAPEZOIDAL = 2,
    SIMPSONS = 3,
};

See below for how to use it.

Don't use (virtual) classes unnecessarily

There is no reason for class NewtonCotesFactory to have virtual functions. Nothing is inheriting from it. In fact, it has no member variables and only a single member function that does anything useful, so it should not be a class to begin with. Just create a stand-alone function:

std::unique_ptr<NewtonCotesFormula> createQuadrature(Formula choice, ...) {
    switch (choice) {
    case Formulas::MIDPOINT:
        return new Midpoint(a, b, inp_tol);
    case Formulas::TRAPEZOIDAL:
        return new Trapezoidal(a, b, inp_tol);
    case Formulas::SIMPSONS:
        return new Simpsons(a, b, inp_tol);
    }
}

It's also possible to restructure the other classes to not need virtual functions. The actual computeIntegral() functions don't depend on much, except for the the sup and inf values. Why not pass those as parameters as well? And again you can write them as stand-alone functions:

template<typename Function>
double Midpoint(double inf, double sup, int intervals, Function func) {
    double interval_width = (sup - inf) / intervals;
    double result = 0;

    for (int i = 0; i < intervals; ++i)
        result += func(inf + (i + 0.5) * interval_width);

    return result * interval_width;
}

Since you are passing functions to other functions already: you can pass both the function you wish to integrate over and the quadrature function to convergeToTol():

template<typename Function, function Quadrature>
double convergeToTol(double inf, double sup, int intervals, double tol, Function func, Quadrature computeIntegral) {
    ...
    while (difference > tol) {
        ...
        new_integral = computeIntegral(inf, sup, intervals, func);
        ...
    }
    ...
}

Now that functions are no longer in classes, the question is how createQuadrature() should be rewritten. One way is to use std::function to be able to return a function:

std::function<double(double, double, int, std::function<double(double)>)>
createQuadrature(Formula choice) {
    switch (choice) {
    case Formulas::MIDPOINT:
        return Midpoint<std::function<double(double)>>;
    case Formulas::TRAPEZOIDAL:
        return Trapezoidal<std::function<double(double)>>;
    case Formulas::SIMPSONS:
        return Simpsons<std::function<double(double)>>;
    }
}

The above is now more in a functional programming style, which I think is more appropriate for the problem.

Answer 3

Prefer Composition to Inheritance

Currently, you have objects that store and hide the interval you’re integrating over, the tolerance, and the method of integration. Each class has a different virtual method that you pass a number of intervals and a C-style pointer to the function to be integrated. The basic idea of a object defining a type of integration method, which you can re-use or store, is sound. (Although you could overload a non-member integration function to take its data members as arguments.)

This won’t work well for an application where you want to try integrating the same function using different methods and tolerances. You really have two somewhat orthogonal parts to the integral: the function you want to integrate plus the interval over which you want to integrate it, and the method you want to use to integrate it, along with its parameters.

Client code should create or receive these two things separately, and compose them together with a function call equivalent to, “Integrate f from a to b, by whatever method this object represents.” You don’t want your API to have to specify and pass along a number-of-intervals parameter separately. Among other problems, that limits the API to only ever handle regular partitions.

You also shouldn’t encapsulate elements of what you’re integrating, like the limits of integration, along with elements of the numeric method for integrating it, such as inp_tol, if you want to be able to either integrate a function over different intervals with the same method, or a function over the same interval with several different methods.

The basic idea of having different Simpsons, Trapezoidal and Midpoint methods seems sound, but in that case, the instances should encapsulate n and inp_tol.

You might also think it worthwhile to have a DefiniteIntegral class that encapsulates the function to be integrated along with the limits of integration. At that point, your API could support any or all of the following ways to compose your two orthogonal objects:

• integrationMethod.computeIntegral(f, a, b);
• f_a_to_b.computeIntegral(integrationMethod);
• integrationMethod.computeIntegral(f_a_to_b);

If you then want an object, let’s say int_f, that evaluates the integral of f from a to b according to a given method, you can then create one with

const auto int_f = std::bind_from( &Midpoint::computeIntegral,
                                   Midpoint( intervals, inp_tol ),
                                   f );
area = int_f( a, b );

Declare Constants Where You Can

Numeric integration is a great use case for a pure computation. A function like NewtonCotesFormulas::computeIntegral should be declared constexpr if some integrals can be computed purely at compile time, and const if it does not modify any per-instance data.

Similarly, you do not want to create a shallow copy of a std::string object that says “midpoint” for every instance, or dynamically allocate a new copy of the string on the heap every time you copy the object. You should actually use typeid to determine the type of these objects at runtime. This will work because they contain at least one virtual function, so a virtual function table will be created for them.

An alternative to storing any per-instance data at all is to write a virtual function that returns a const reference to a private: static const std::string object. For example,

// Declared as static in the definition of Midpoint:
const std::string Midpoint::method_name = "midpoint";

// Declared as virtual in the base class:    
constexpr const std::string& Midpoint::get_methodname() noexcept
{
  return Midpoint::method_name; // Good practice to specify the namespace when re-using identifiers.
}

This uses no additional memory (beyond one more pointer in the virtual function table) and runs with zero overhead when the type of the object is known at compile time.

However, if you do want to have a data member in each instance that contains the name of the method, it should be a pointer to a static string constant. You will never want to customize or modify constants such as "Simpson's" or "trapezoidal". All copies of the string are then shallow copies, which not only is much more efficient, it allows constexpr instances to exist, potentially allowing the compiler to pre-calculate integrals whose parameters are all compile-time constants, and greatly optimizing your program.

Pass the Right Types

You should pass the function to be integrated as a std::function object instead of a C-style function pointer. This enables you to use many more types of functions, such as a member call to an object, a lambda expression, or an object such as std::bind( normalDistribution, mean, variance ).

You pass the number of intervals as an int, which might be zero or negative, then never check that the domain is valid. Since an int is not necessarily 32 bits wide, you might want to declare the type as unsigned long or perhaps std::uint32_t from <cstdint>. Then if you declare what requesting zero intervals means, there would be no possible domain errors Or, if the number of partitions is encapsulated in the integration method, it can be checked and throw an exception when set to an invalid value, and then the critical path that you want to optimize can be guaranteed safe, and declared noexcept.

Shut Down the Factory

Just use the constructor. The factory pattern would be useful if you want to, for example, create a singleton as needed and return a reference to it. You’re not doing anything like that here, but creating individual instances with their own data.

Instead of client code writing

constexpr int QUAD_MIDPOINT = 1;
NewtonCotesFormulas* const quadrature =
  NewtonCotesFactory::createQuadrature(QUAD_MIDPOINT, a, b, inp_tol );

You can instead write

const auto& quadrature = Midpoint( a, b, inp_tol );

or, if you need virtual function dispatch,

const NewtonCotesFormulas& quadrature =
  Midpoint( a, b, inp_tol );

This works fine, because an object of a derived class implicitly converts to a reference to the base class. Access through the reference to an abstract base class will properly dispatch all virtual functions.

If you decide you do need a factory, it should return either a reference to a static object, such as a singleton—which should never return NULL—or a std::unique_ptr<NewtonCotesFormulas>. (In many use cases, there is a good third option, returning an object that will be created with guaranteed copy elision, but that is not possible here because you want a reference to a pure virtual class.) The smart pointer will automatically free the object it owns, once and only once, when its lifetime expires, and a naked new will not.

If you do want a factory method that creates an object of a derived class dynamically and returns it as a smart pointer to a base class, which would be useful if you need to store references to the abstract interface in a data structure, a possible implementation would be:

template<class T, class... Args>
  std::unique_ptr<NewtonCotesFormulas> make_NewtonCotes( Args&&... args )
  {
    return std::unique_ptr<NewtonCotesFormulas>(
      static_cast<NewtonCotesFormulas*>(new T(args...)));
  }

const auto mid = make_NewtonCotes<Midpoint>( a, b, inp_tol );
const double area = mid->computeIntegral( n, f );

While that sample gives you a std::unique_ptr<NewtonCotesFormulas> object, you can also initialize different kinds of smart pointers from it:

const std::shared_ptr<NewtonCotesFormulas> shared = make_NewtonCotes<Midpoint>( a, b, inp_tol );
const std::unique_ptr<const NewtonCotesFormulas> cmid = make_NewtonCotes<Midpoint>( a, b, inp_tol );

These work for any type of object that can be constructed or implicitly converted from a std::unique_ptr<NewtonCotesMethod>. It would not work for a std::weak_ptr, for example.

This function template constructs a new object, dynamically, of the class you specify (Midpoint), forwarding the arguments (here, a, b and inp_tol, but any arguments will work) to the constructor of the derived class. If any of them can be moved instead of copied, it will do that. It wraps the created object in a smart pointer to the base class, which will manage its memory automatically. Dereferencing the smart pointer will use virtual function dispatch.

Since this removes the need for all derived classes to have the same type of constructor—you could call make_NewtonCoates<MethodA>(foo) or make_NewtonCoates<MethodB>(foo, bar, baz)—this probably removes any need for the abstract base class to specify any data members. Each daughter class can define whatever ones it needs.

Consider Wrapping Your Library in a Namespace

If you declare all identifiers in namespace NumericIntegration and then using namespace NumericIntegration, you’ll get the best of both worlds: the namespace is completely transparent to your clients, but if there is, unlikely as it might seem, another library that declares an identifier with the same name, you can still link your program, and write NumericIntegration::foo to disambiguate as needed.

Listen to the Other Answers Too

Don’t #include a .cpp file, avoid using namespace std; use the default constructors and destructors instead of empty braes, and so on. I won’t repeat all of that here.

Answer 4

I have adapted my code in the original question to accommodate for the more functional approach suggested by @G-Sliepen in his answer.

Also a timing subroutine was added here and to the old OOP code, in order to compare performances.

Functional version of code implementing Newton-Cotes quadrature formulas:

main.cpp:

#include <iostream>
#include <cmath>
#include <cctype>
#include <functional>
#include <chrono>
#include "NewtonCotesFormulas.h"

using std::cin;     using std::cout;
using std::pow;     using std::function;
using std::exp;
namespace chrono = std::chrono;

// helps with input
void flushCin() {
    cin.clear();
    cin.ignore();
}

// timer
unsigned int stopwatch()
{
    static auto start_time = chrono::steady_clock::now();

    auto end_time = chrono::steady_clock::now();
    auto delta = chrono::duration_cast<chrono::microseconds>(end_time - start_time);

    start_time = end_time;

    return delta.count();
}

// use a functor this time
struct func {
    double operator()(double x) {
        return sqrt(x) * exp(-x);
    }
};
//double func(double x) {
//    return sqrt(x) * exp(-x);
//}

int main() {
    cout << "Choose method:\n"
         << "Midpoint numeric quadrature -----> (1)\n"
         << "Trapezoidal numeric quadrature --> (2)\n"
         << "Simpson's numeric quadrature ----> (3)\n\n";

    // user picks a quadrature method
//    int choice_int;
//    cin >> choice_int;
//    while (choice_int < 1 || choice_int > 3) {
//        flushCin();
//        cout << "Pick again a number from 1 to 3\n";
//        cin >> choice_int;
//    }
//    Formula choice = static_cast<Formula>(choice_int);

    // cycle through all quadrature formulas
    stopwatch();
    for (Formula choice : {Formula::Midpoint, Formula::Trapezoidal, Formula::Simpsons}) {
        double result = printConvergenceValues(choice, 1, 3, 4, 0.000001, func() );
                           // above parameters:choice, inf, sup, intervals, tol, func
        cout << "Final result within tolerance: " << result << "\n\n";
    }
    cout << "Performance time: " << stopwatch() << " microseconds.\n\n";

    return 0;
}

NewtonCotesFormulas.h:

#ifndef NEWTONCOTESFORMULAS_H_INCLUDED
#define NEWTONCOTESFORMULAS_H_INCLUDED

#include <iostream>
#include <functional>

enum class Formula {
    Midpoint = 1,
    Trapezoidal,
    Simpsons
};

std::ostream& operator<< (std::ostream&, const Formula&);

template<typename Function, typename Quadrature>
double convergeToTol(const double, const double, int, const double, const Function&, const Quadrature&);

template<typename Function>
double printConvergenceValues(const Formula&, const double, const double, const int, const double, const Function&);

template<typename Function> double Midpoint(const double, const double, const int, const Function&);

template<typename Function> double Trapezoidal(const double, const double, const int, const Function&);

template<typename Function> double Simpsons(const double, const double, const int, const Function&);

std::function<double (double, double, int, std::function<double (double) >) >
createQuadrature(const Formula&);

#include "NewtonCotesFormulas.cpp"
#endif // NEWTONCOTESFORMULAS_H_INCLUDED

NewtonCotesFormulas.cpp:

#include <iomanip>
#include <string>
#include <cmath>
#include <stdexcept>
#include "NewtonCotesFormulas.h"

using std::string;          using std::cout;
using std::setprecision;    using std::fixed;
using std::scientific;      using std::to_string;
using std::abs;             using std::function;
using std::runtime_error;   using std::ostream;

// this prints the enum on cout
ostream& operator<< (ostream& os, const Formula& method) {
    if (method == Formula::Midpoint)
        os << "Midpoint";
    else if (method == Formula::Trapezoidal)
        os << "Trapezoidal";
    else if (method == Formula::Simpsons)
        os << "Simpsons";

    return os;
}

template<typename Function, typename Quadrature>
double convergeToTol(const double inf, const double sup, int intervals, const double tol,
                     const Function& func, const Quadrature& my_method) {
    double old_integral = my_method(inf, sup, intervals, func);
    intervals *= 2;
    double new_integral = my_method(inf, sup, intervals, func);
    double difference = abs(new_integral - old_integral);

    cout.setf(ios_base::scientific);
    cout.precision(8);
    while(difference > tol) {
        cout << string(9 - to_string(intervals).size(), ' ') << intervals <<
                string(3,' ') << fixed << new_integral << string(3,' ') <<
                scientific << difference << '\n';
        old_integral = new_integral;
        intervals *= 2;
        new_integral = my_method(inf, sup, intervals, func);
        difference = abs(new_integral - old_integral);
    }
    cout << fixed << string(9 - to_string(intervals).size(), ' ') <<
            intervals << string(3,' ') << new_integral << string(3,' ') <<
            scientific << difference << '\n';

    return new_integral;
}

template<typename Function>
double printConvergenceValues(const Formula& choice, const double inf, const double sup, 
                              const int intervals, const double tol, const Function& func)
{
    function<double (double, double, int, function<double (double) >) >
        my_method = createQuadrature(choice);

    cout << choice << " method, with tolerance " << tol << "\n";
    cout << "intervals" << string(5,' ') << "integral " << string(9,' ') << "tol\n";
    cout << setprecision(8) << fixed << string(9 - to_string(intervals).size(), ' ') <<
            intervals << string(3,' ') << my_method(inf, sup, intervals, func) <<
            string(3,' ') << 0.0 << '\n';

    double result = convergeToTol(inf, sup, intervals, tol, func, my_method);

    return result;
}

template <typename Function>
double Midpoint(const double inf, const double sup, const int intervals, const Function& func) {
    double interval_width = (sup - inf) / intervals;
    double result = 0;
    for (int i = 1; i <= intervals; ++i)
        result += func(inf + (i - 0.5) * interval_width);

    return result * interval_width;
}

template <typename Function>
double Trapezoidal(const double inf, const double sup, const int intervals, const Function& func) {
    double interval_width = (sup - inf) / intervals;
    double result = (func(inf) + func(sup)) / 2.;
    for (int i = 1; i <= intervals - 1; ++i)
        result += func(inf + i * interval_width);

    return result * interval_width;
}

template <typename Function>
double Simpsons(const double inf, const double sup, const int intervals, const Function& func) {
    double interval_width = (sup - inf) / intervals;
    double result = (func(inf) + func(sup)) / 6;
    for (int i = 1; i <= intervals - 1; ++i)
        result += func(inf + i * interval_width) / 3 + 2 * func(inf + (i - 0.5) * interval_width) / 3;
    result += 2 * func(inf + (intervals - 0.5) * interval_width) / 3;

    return result * interval_width;
}

function<double (double, double, int, function<double (double) >) >
createQuadrature(const Formula& choice) {
    if (choice == Formula::Midpoint)
        return Midpoint<function<double (double)> >;
    else if (choice == Formula::Trapezoidal)
        return Trapezoidal<function<double (double)>>;
    else if (choice == Formula::Simpsons)
        return Simpsons<function<double (double)>>;
    else
        throw runtime_error("Pick between choice = 1, 2, or 3");

    return 0;  // never gets here
}

Numerical results are obviously identical among the two versions. Performances are as following:

• Old OO code:

    Performance time: 23375 microseconds.
  

• Functional approach, code above:

    Performance time: 22211 microseconds.
  

I'd say performances are similar. I suppose the introduction of templates balances out the removal of virtual calls, in terms of slowing down computational times?

The functional approach is surely more elegant, but also a bit more complicate to write, and read as well. Or perhaps it's just me not being overly familiar with functional programming.