Composite Simpson's Rule in C++

Question

I've been trying to write a function to approximate an the value of an integral using the Composite Simpson's Rule.

template <typename func_type>
double simp_rule(double a, double b, int n, func_type f){

    int i = 1; double area = 0;
    double n2 = n;
    double h = (b-a)/(n2-1), x=a;

    while(i <= n){

        area = area + f(x)*pow(2,i%2 + 1)*h/3;
        x+=h;
        i++;
    }
    area -= (f(a) * h/3);
    area -= (f(b) * h/3);

    return area;
    }

What I do is multiply each value of the function by either 2 or 4 (and h/3) with pow(2,i%2 + 1) and subtract off the edges as these should only have a weight of 1.

At first, I thought it worked just fine, however, when I compared it to my Trapezoidal Method function it was way more inaccurate which shouldn't be the case.

This is a simpler version of a code I previously wrote which had the same problem, I thought that if I cleaned it up a little the problem would go away, but alas. From another post, I get the idea that there's something going on with the types and the operations I'm doing on them which results in loss of precision, but I just don't see it.

Edit:

For completeness, I was running it for e^x from 1 to zero

\\function to be approximated
double f(double x){ double a = exp(x); return a; }

int main() {

    int n = 11; //this method works best for odd values of n
    double e = exp(1);
    double exact = e-1; //value of integral of e^x from 0 to 1

    cout << simp_rule(0,1,n,f) - exact;

use a debugger and a pocket calculator on the lines with h/3 and see what happens. is the runtime type of `f` return value also double? — Cee McSharpface, Jan 31 '20 at 14:15
@CeeMcSharpface I'm not sure I understood your question, but the value I get for n=11 is within 0.1 of the actual value. The difference should be about 10^-7 I belive. — AJ09, Jan 31 '20 at 14:21
as `func_type` is templated, can we be sure it is a `double` at runtime? consider including a minimal, complete example with a sample function so we could run and debug it? — Cee McSharpface, Jan 31 '20 at 14:25
@AlejandroMatos I had simplified your function and cleared somethings, then used it to integrate `x^2` on `1->5` using 100 steps and it returned `41.3463`. So, provide the code giving you undesired results to see. — asmmo, Jan 31 '20 at 14:36
You should have `i=0` to initialize and have the loop run as `while(i < n)`. The debugger could have showed you that your weights were going 4,2,4,2,4 instead of 2,4,2,4,2 — SirGuy, Jan 31 '20 at 14:37
@SirGuy It worked!! Thank you! If you write this as an answer I'll give you the check mark — AJ09, Jan 31 '20 at 16:21

Bob__ · Accepted Answer · 2020-02-01T16:13:03.077

The Simpson's Rule uses this approximation to estimate a definite integral:

Where

and

So that there are n + 1 equally spaced sample points x_i.

In the posted code, the parameter n passed to the function appears to be the number of points where the function is sampled (while in the previous formula n is the number of intervals, that's not a problem).

The (constant) distance between the points is calculated correctly

double h = (b - a) / (n - 1);

The while loop used to sum the weighted contributes of all the points iterates from x = a up to a point with an ascissa close to b, but probably not exactly b, due to rounding errors. This implies that the last calculated value of f, f(x_n), may be slightly different from the expected f(b).

This is nothing, though, compared to the error caused by the fact that those end points are summed inside the loop with the starting weight of 4 and then subtracted after the loop with weight 1, while all the inner points have their weight switched. As a matter of fact, this is what the code calculates:

$\frac{\Delta x}{3}\left ( 3f(x_0)+ 2f(x_1) + 4f(x_2) + ... + 2f(x_{n-1}) + 3f(x_{n}) \right )$

Also, using

pow(2, i%2 + 1)

To generate the sequence 4, 2, 4, 2, ..., 4 is a waste, in terms of efficency, and may add (depending on the implementation) other unnecessary rounding errors.

The following algorithm shows how to obtain the same (fixed) result, without a call to that library function.

template <typename func_type>
double simpson_rule(double a, double b,
                    int n, // Number of intervals
                    func_type f)
{
    double h = (b - a) / n;

    // Internal sample points, there should be n - 1 of them
    double sum_odds = 0.0;
    for (int i = 1; i < n; i += 2)
    {
        sum_odds += f(a + i * h);
    }
    double sum_evens = 0.0;
    for (int i = 2; i < n; i += 2)
    {
        sum_evens += f(a + i * h);
    }

    return (f(a) + f(b) + 2 * sum_evens + 4 * sum_odds) * h / 3;
}

Note that this function requires the number of intervals (e.g. use 10 instead of 11 to obtain the same results of OP's function) to be passed, not the number of points.

Testable here.

Thanks for your response. I have the end points in there from the start, they have a weight of 2 though which is why I subtract them at the end. Also, could you please elaborate on why using `pow(2,i%2+1)` is not good? — AJ09, Jan 31 '20 at 20:07
@AlejandroMatos Edited to clarify. The problem with `pow` is that it's just misused for tasks like this. It's potentially (it depends on the particular implementation) much more expansive than other simple workarounds. While it's often used in math textbooks as a great tool to generalize or simplify formulas, it may be the wrong tool when it comes to a numeric implementation in a computer program. — Bob__, Jan 31 '20 at 21:39

score 2 · Answer 2 · answered Apr 07 '20 at 17:44

The above excellent and accepted solution could benefit from liberal use of std::fma() and templatize on the floating point type. https://en.cppreference.com/w/cpp/numeric/math/fma

#include <cmath>
template <typename fptype, typename func_type>
double simpson_rule(fptype a, fptype b,
                    int n, // Number of intervals
                    func_type f)
{
    fptype h = (b - a) / n;

    // Internal sample points, there should be n - 1 of them
    fptype sum_odds = 0.0;
    for (int i = 1; i < n; i += 2)
    {
        sum_odds += f(std::fma(i,h,a));
    }
    fptype sum_evens = 0.0;
    for (int i = 2; i < n; i += 2)
    {
        sum_evens += f(std::fma(i,h,a);
    }

    return (std::fma(2,sum_evens,f(a)) + 
            std::fma(4,sum_odds,f(b))) * h / 3;
}

Composite Simpson's Rule in C++

2 Answers2