composite_simpson_rule.cpp File Reference

Implementation of the Composite Simpson Rule for the approximation. More...

#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <functional>
#include <iostream>
#include <map>

Include dependency graph for composite_simpson_rule.cpp:

Namespaces
namespace	numerical_methods
	for assert
namespace	simpson_method
	Contains the Simpson's method implementation.

Functions
double	numerical_methods::simpson_method::evaluate_by_simpson (std::int32_t N, double h, double a, const std::function< double(double)> &func)
double	numerical_methods::simpson_method::f (double x)
	A function f(x) that will be used to test the method. More...
double	numerical_methods::simpson_method::g (double x)
	Another test function. More...
double	numerical_methods::simpson_method::k (double x)
	Another test function. More...
double	numerical_methods::simpson_method::l (double x)
	Another test function. More...
static void	test (std::int32_t N, double h, double a, double b, bool used_argv_parameters)
	Self-test implementations. More...
int	main (int argc, char **argv)
	Main function. More...

Detailed Description

Implementation of the Composite Simpson Rule for the approximation.

The following is an implementation of the Composite Simpson Rule for the approximation of definite integrals. More info -> wiki: https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule

The idea is to split the interval in an EVEN number N of intervals and use as interpolation points the xi for which it applies that xi = x0 + i*h, where h is a step defined as h = (b-a)/N where a and b are the first and last points of the interval of the integration [a, b].

We create a table of the xi and their corresponding f(xi) values and we evaluate the integral by the formula: I = h/3 * {f(x0) + 4*f(x1) + 2*f(x2) + ... + 2*f(xN-2) + 4*f(xN-1) + f(xN)}

That means that the first and last indexed i f(xi) are multiplied by 1, the odd indexed f(xi) by 4 and the even by 2.

In this program there are 4 sample test functions f, g, k, l that are evaluated in the same interval.

Arguments can be passed as parameters from the command line argv[1] = N, argv[2] = a, argv[3] = b

N must be even number and a<b.

In the end of the main() i compare the program's result with the one from mathematical software with 2 decimal points margin.

Add sample function by replacing one of the f, g, k, l and the assert

Author

ggkogkou

Function Documentation

evaluate_by_simpson()

double numerical_methods::simpson_method::evaluate_by_simpson	(	std::int32_t	N,
		double	h,
		double	a,
		const std::function< double(double)> &	func
	)

                                                                    {
    std::map<std::int32_t, double>
        data_table;  // Contains the data points. key: i, value: f(xi)
    double xi = a;   // Initialize xi to the starting point x0 = a
 
    // Create the data table
    double temp = NAN;
    for (std::int32_t i = 0; i <= N; i++) {
        temp = func(xi);
        data_table.insert(
            std::pair<std::int32_t, double>(i, temp));  // add i and f(xi)
        xi += h;  // Get the next point xi for the next iteration
    }
 
    // Evaluate the integral.
    // Remember: f(x0) + 4*f(x1) + 2*f(x2) + ... + 2*f(xN-2) + 4*f(xN-1) + f(xN)
    double evaluate_integral = 0;
    for (std::int32_t i = 0; i <= N; i++) {
        if (i == 0 || i == N) {
            evaluate_integral += data_table.at(i);
        } else if (i % 2 == 1) {
            evaluate_integral += 4 * data_table.at(i);
        } else {
            evaluate_integral += 2 * data_table.at(i);
        }
    }
 
    // Multiply by the coefficient h/3
    evaluate_integral *= h / 3;
 
    // If the result calculated is nan, then the user has given wrong input
    // interval.
    assert(!std::isnan(evaluate_integral) &&
           "The definite integral can't be evaluated. Check the validity of "
           "your input.\n");
    // Else return
    return evaluate_integral;
}

f()

double numerical_methods::simpson_method::f

(

double

x

)

A function f(x) that will be used to test the method.

Parameters

x	The independent variable xi

Returns

the value of the dependent variable yi = f(xi)

{ return std::sqrt(x) + std::log(x); }

Here is the call graph for this function:

g()

double numerical_methods::simpson_method::g

(

double

x

)

Another test function.

{ return std::exp(-x) * (4 - std::pow(x, 2)); }

Here is the call graph for this function:

k()

double numerical_methods::simpson_method::k

(

double

x

)

Another test function.

Examples

/Users/runner/work/C-Plus-Plus/C-Plus-Plus/numerical_methods/rungekutta.cpp.

{ return std::sqrt(2 * std::pow(x, 3) + 3); }

Here is the call graph for this function:

l()

double numerical_methods::simpson_method::l

(

double

x

)

Another test function.

{ return x + std::log(2 * x + 1); }

Here is the call graph for this function:

main()

int main	(	int	argc,
		char **	argv
	)

Main function.

Parameters

argc	commandline argument count (ignored)
argv	commandline array of arguments (ignored)

Returns

0 on exit

Number of intervals to divide the integration interval. MUST BE EVEN

Starting and ending point of the integration in the real axis

Step, calculated by a, b and N

                                {
    std::int32_t N = 16;  /// Number of intervals to divide the integration
                          /// interval. MUST BE EVEN
    double a = 1, b = 3;  /// Starting and ending point of the integration in
                          /// the real axis
    double h = NAN;       /// Step, calculated by a, b and N
 
    bool used_argv_parameters =
        false;  // If argv parameters are used then the assert must be omitted
                // for the tst cases
 
    // Get user input (by the command line parameters or the console after
    // displaying messages)
    if (argc == 4) {
        N = std::atoi(argv[1]);
        a = std::atof(argv[2]);
        b = std::atof(argv[3]);
        // Check if a<b else abort
        assert(a < b && "a has to be less than b");
        assert(N > 0 && "N has to be > 0");
        if (N < 16 || a != 1 || b != 3) {
            used_argv_parameters = true;
        }
        std::cout << "You selected N=" << N << ", a=" << a << ", b=" << b
                  << std::endl;
    } else {
        std::cout << "Default N=" << N << ", a=" << a << ", b=" << b
                  << std::endl;
    }
 
    // Find the step
    h = (b - a) / N;
 
    test(N, h, a, b, used_argv_parameters);  // run self-test implementations
 
    return 0;
}

Here is the call graph for this function:

test()

static void test ( std::int32_t  N, double  h, double  a, double  b, bool  used_argv_parameters  )

static

Self-test implementations.

Parameters

N	is the number of intervals
h	is the step
a	is x0
b	is the end of the interval
used_argv_parameters	is 'true' if argv parameters are given and 'false' if not

                                            {
    // Call the functions and find the integral of each function
    double result_f = numerical_methods::simpson_method::evaluate_by_simpson(
        N, h, a, numerical_methods::simpson_method::f);
    assert((used_argv_parameters || (result_f >= 4.09 && result_f <= 4.10)) &&
           "The result of f(x) is wrong");
    std::cout << "The result of integral f(x) on interval [" << a << ", " << b
              << "] is equal to: " << result_f << std::endl;
 
    double result_g = numerical_methods::simpson_method::evaluate_by_simpson(
        N, h, a, numerical_methods::simpson_method::g);
    assert((used_argv_parameters || (result_g >= 0.27 && result_g <= 0.28)) &&
           "The result of g(x) is wrong");
    std::cout << "The result of integral g(x) on interval [" << a << ", " << b
              << "] is equal to: " << result_g << std::endl;
 
    double result_k = numerical_methods::simpson_method::evaluate_by_simpson(
        N, h, a, numerical_methods::simpson_method::k);
    assert((used_argv_parameters || (result_k >= 9.06 && result_k <= 9.07)) &&
           "The result of k(x) is wrong");
    std::cout << "The result of integral k(x) on interval [" << a << ", " << b
              << "] is equal to: " << result_k << std::endl;
 
    double result_l = numerical_methods::simpson_method::evaluate_by_simpson(
        N, h, a, numerical_methods::simpson_method::l);
    assert((used_argv_parameters || (result_l >= 7.16 && result_l <= 7.17)) &&
           "The result of l(x) is wrong");
    std::cout << "The result of integral l(x) on interval [" << a << ", " << b
              << "] is equal to: " << result_l << std::endl;
}