C-Plus-Plus/math/integral_approximation.cpp

/**
 * @file
 * @brief Compute integral approximation of the function using [Riemann
 * sum](https://en.wikipedia.org/wiki/Riemann_sum)
 * @details In mathematics, a Riemann sum is a certain kind of approximation of
 * an integral by a finite sum. It is named after nineteenth-century German
 * mathematician Bernhard Riemann. One very common application is approximating
 * the area of functions or lines on a graph and the length of curves and other
 * approximations. The sum is calculated by partitioning the region into shapes
 * (rectangles, trapezoids, parabolas, or cubics) that form a region similar to
 * the region being measured, then calculating the area for each of these
 * shapes, and finally adding all of these small areas together. This approach
 * can be used to find a numerical approximation for a definite integral even if
 * the fundamental theorem of calculus does not make it easy to find a
 * closed-form solution. Because the region filled by the small shapes is
 * usually not the same shape as the region being measured, the Riemann sum will
 * differ from the area being measured. This error can be reduced by dividing up
 * the region more finely, using smaller and smaller shapes. As the shapes get
 * smaller and smaller, the sum approaches the Riemann integral. \author
 * [Benjamin Walton](https://github.com/bwalton24) \author [Shiqi
 * Sheng](https://github.com/shiqisheng00)
 */
#include <cassert>     /// for assert
#include <cmath>       /// for mathematical functions
#include <cstdint>
#include <functional>  /// for passing in functions
#include <iostream>    /// for IO operations

/**
 * @namespace math
 * @brief Mathematical functions
 */
namespace math {
/**
 * @brief Computes integral approximation
 * @param lb lower bound
 * @param ub upper bound
 * @param func function passed in
 * @param delta
 * @returns integral approximation of function from [lb, ub]
 */
double integral_approx(double lb, double ub,
                       const std::function<double(double)>& func,
                       double delta = .0001) {
    double result = 0;
    uint64_t numDeltas = static_cast<uint64_t>((ub - lb) / delta);
    for (int i = 0; i < numDeltas; i++) {
        double begin = lb + i * delta;
        double end = lb + (i + 1) * delta;
        result += delta * (func(begin) + func(end)) / 2;
    }
    return result;
}

/**
 * @brief Wrapper to evaluate if the approximated
 * value is within `.XX%` threshold of the exact value.
 * @param approx aprroximate value
 * @param exact expected value
 * @param threshold values from [0, 1)
 */
void test_eval(double approx, double expected, double threshold) {
    assert(approx >= expected * (1 - threshold));
    assert(approx <= expected * (1 + threshold));
}

/**
 * @brief Self-test implementations to
 * test the `integral_approx` function.
 *
 * @returns `void`
 */
}  // namespace math

static void test() {
    double test_1 = math::integral_approx(
        3.24, 7.56, [](const double x) { return log(x) + exp(x) + x; });
    std::cout << "Test Case 1" << std::endl;
    std::cout << "function: log(x) + e^x + x" << std::endl;
    std::cout << "range: [3.24, 7.56]" << std::endl;
    std::cout << "value: " << test_1 << std::endl;
    math::test_eval(test_1, 1924.80384023549, .001);
    std::cout << "Test 1 Passed!" << std::endl;
    std::cout << "=====================" << std::endl;

    double test_2 = math::integral_approx(0.023, 3.69, [](const double x) {
        return x * x + cos(x) + exp(x) + log(x) * log(x);
    });
    std::cout << "Test Case 2" << std::endl;
    std::cout << "function: x^2 + cos(x) + e^x + log^2(x)" << std::endl;
    std::cout << "range: [.023, 3.69]" << std::endl;
    std::cout << "value: " << test_2 << std::endl;
    math::test_eval(test_2, 58.71291345202729, .001);
    std::cout << "Test 2 Passed!" << std::endl;
    std::cout << "=====================" << std::endl;

    double test_3 = math::integral_approx(
        10.78, 24.899, [](const double x) { return x * x * x - x * x + 378; });
    std::cout << "Test Case 3" << std::endl;
    std::cout << "function: x^3 - x^2 + 378" << std::endl;
    std::cout << "range: [10.78, 24.899]" << std::endl;
    std::cout << "value: " << test_3 << std::endl;
    math::test_eval(test_3, 93320.65915078377, .001);
    std::cout << "Test 3 Passed!" << std::endl;
    std::cout << "=====================" << std::endl;

    double test_4 = math::integral_approx(
        .101, .505,
        [](const double x) { return cos(x) * tan(x) * x * x + exp(x); },
        .00001);
    std::cout << "Test Case 4" << std::endl;
    std::cout << "function: cos(x)*tan(x)*x^2 + e^x" << std::endl;
    std::cout << "range: [.101, .505]" << std::endl;
    std::cout << "value: " << test_4 << std::endl;
    math::test_eval(test_4, 0.566485986311631, .001);
    std::cout << "Test 4 Passed!" << std::endl;
    std::cout << "=====================" << std::endl;

    double test_5 = math::integral_approx(
        -1, 1, [](const double x) { return exp(-1 / (x * x)); });
    std::cout << "Test Case 5" << std::endl;
    std::cout << "function: e^(-1/x^2)" << std::endl;
    std::cout << "range: [-1, 1]" << std::endl;
    std::cout << "value: " << test_5 << std::endl;
    math::test_eval(test_5, 0.1781477117815607, .001);
    std::cout << "Test 5 Passed!" << std::endl;
}

/**
 * @brief Main function
 * @returns 0 on exit
 */
int main() {
    test();  // run self-test implementations
    return 0;
}