composite_simpson_rule.cpp

Go to the documentation of this file.

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

namespace numerical_methods {
namespace simpson_method {
double 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;
}

double f(double x) { return std::sqrt(x) + std::log(x); }
double g(double x) { return std::exp(-x) * (4 - std::pow(x, 2)); }
double k(double x) { return std::sqrt(2 * std::pow(x, 3) + 3); }
double l(double x) { return x + std::log(2 * x + 1); }
}  // namespace simpson_method
}  // namespace numerical_methods

static void test(std::int32_t N, double h, double a, double b,
                 bool used_argv_parameters) {
    // 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;
}

int main(int argc, char** argv) {
    std::int32_t N = 16;  
    double a = 1, b = 3;  
    double h = NAN;       
 
    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;
}