|
| 1 | +/** |
| 2 | + * @file |
| 3 | + * @brief Compute integral approximation of the function using [Riemann sum](https://en.wikipedia.org/wiki/Riemann_sum) |
| 4 | + * @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. |
| 5 | + * One very common application is approximating the area of functions or lines on a graph and the length of curves and other approximations. |
| 6 | + * 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. |
| 7 | + * 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. |
| 8 | + * 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. |
| 9 | + * 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. |
| 10 | + * \author [Benjamin Walton](https://github.com/bwalton24) |
| 11 | + * \author [Shiqi Sheng](https://github.com/shiqisheng00) |
| 12 | + */ |
| 13 | +#include <cassert> /// for assert |
| 14 | +#include <cmath> /// for mathematical functions |
| 15 | +#include <functional> /// for passing in functions |
| 16 | +#include <iostream> /// for IO operations |
| 17 | + |
| 18 | +/** |
| 19 | + * @namespace math |
| 20 | + * @brief Mathematical functions |
| 21 | + */ |
| 22 | +namespace math { |
| 23 | +/** |
| 24 | + * @brief Computes integral approximation |
| 25 | + * @param lb lower bound |
| 26 | + * @param ub upper bound |
| 27 | + * @param func function passed in |
| 28 | + * @param delta |
| 29 | + * @returns integral approximation of function from [lb, ub] |
| 30 | + */ |
| 31 | +double integral_approx(double lb, double ub, |
| 32 | + const std::function<double(double)>& func, |
| 33 | + double delta = .0001) { |
| 34 | + double result = 0; |
| 35 | + uint64_t numDeltas = static_cast<uint64_t>((ub - lb) / delta); |
| 36 | + for (int i = 0; i < numDeltas; i++) { |
| 37 | + double begin = lb + i * delta; |
| 38 | + double end = lb + (i + 1) * delta; |
| 39 | + result += delta * (func(begin) + func(end)) / 2; |
| 40 | + } |
| 41 | + return result; |
| 42 | +} |
| 43 | + |
| 44 | +/** |
| 45 | + * @brief Wrapper to evaluate if the approximated |
| 46 | + * value is within `.XX%` threshold of the exact value. |
| 47 | + * @param approx aprroximate value |
| 48 | + * @param exact expected value |
| 49 | + * @param threshold values from [0, 1) |
| 50 | + */ |
| 51 | +void test_eval(double approx, double expected, double threshold) { |
| 52 | + assert(approx >= expected * (1 - threshold)); |
| 53 | + assert(approx <= expected * (1 + threshold)); |
| 54 | +} |
| 55 | + |
| 56 | +/** |
| 57 | + * @brief Self-test implementations to |
| 58 | + * test the `integral_approx` function. |
| 59 | + * |
| 60 | + * @returns `void` |
| 61 | + */ |
| 62 | +} // namespace math |
| 63 | + |
| 64 | +static void test() { |
| 65 | + double test_1 = math::integral_approx( |
| 66 | + 3.24, 7.56, [](const double x) { return log(x) + exp(x) + x; }); |
| 67 | + std::cout << "Test Case 1" << std::endl; |
| 68 | + std::cout << "function: log(x) + e^x + x" << std::endl; |
| 69 | + std::cout << "range: [3.24, 7.56]" << std::endl; |
| 70 | + std::cout << "value: " << test_1 << std::endl; |
| 71 | + math::test_eval(test_1, 1924.80384023549, .001); |
| 72 | + std::cout << "Test 1 Passed!" << std::endl; |
| 73 | + std::cout << "=====================" << std::endl; |
| 74 | + |
| 75 | + double test_2 = math::integral_approx(0.023, 3.69, [](const double x) { |
| 76 | + return x * x + cos(x) + exp(x) + log(x) * log(x); |
| 77 | + }); |
| 78 | + std::cout << "Test Case 2" << std::endl; |
| 79 | + std::cout << "function: x^2 + cos(x) + e^x + log^2(x)" << std::endl; |
| 80 | + std::cout << "range: [.023, 3.69]" << std::endl; |
| 81 | + std::cout << "value: " << test_2 << std::endl; |
| 82 | + math::test_eval(test_2, 58.71291345202729, .001); |
| 83 | + std::cout << "Test 2 Passed!" << std::endl; |
| 84 | + std::cout << "=====================" << std::endl; |
| 85 | + |
| 86 | + double test_3 = math::integral_approx( |
| 87 | + 10.78, 24.899, [](const double x) { return x * x * x - x * x + 378; }); |
| 88 | + std::cout << "Test Case 3" << std::endl; |
| 89 | + std::cout << "function: x^3 - x^2 + 378" << std::endl; |
| 90 | + std::cout << "range: [10.78, 24.899]" << std::endl; |
| 91 | + std::cout << "value: " << test_3 << std::endl; |
| 92 | + math::test_eval(test_3, 93320.65915078377, .001); |
| 93 | + std::cout << "Test 3 Passed!" << std::endl; |
| 94 | + std::cout << "=====================" << std::endl; |
| 95 | + |
| 96 | + double test_4 = math::integral_approx( |
| 97 | + .101, .505, |
| 98 | + [](const double x) { return cos(x) * tan(x) * x * x + exp(x); }, |
| 99 | + .00001); |
| 100 | + std::cout << "Test Case 4" << std::endl; |
| 101 | + std::cout << "function: cos(x)*tan(x)*x^2 + e^x" << std::endl; |
| 102 | + std::cout << "range: [.101, .505]" << std::endl; |
| 103 | + std::cout << "value: " << test_4 << std::endl; |
| 104 | + math::test_eval(test_4, 0.566485986311631, .001); |
| 105 | + std::cout << "Test 4 Passed!" << std::endl; |
| 106 | + std::cout << "=====================" << std::endl; |
| 107 | + |
| 108 | + double test_5 = math::integral_approx( |
| 109 | + -1, 1, [](const double x) { return exp(-1 / (x * x)); }); |
| 110 | + std::cout << "Test Case 5" << std::endl; |
| 111 | + std::cout << "function: e^(-1/x^2)" << std::endl; |
| 112 | + std::cout << "range: [-1, 1]" << std::endl; |
| 113 | + std::cout << "value: " << test_5 << std::endl; |
| 114 | + math::test_eval(test_5, 0.1781477117815607, .001); |
| 115 | + std::cout << "Test 5 Passed!" << std::endl; |
| 116 | +} |
| 117 | + |
| 118 | +/** |
| 119 | + * @brief Main function |
| 120 | + * @returns 0 on exit |
| 121 | + */ |
| 122 | +int main() { |
| 123 | + test(); // run self-test implementations |
| 124 | + return 0; |
| 125 | +} |
0 commit comments