riemann.cc

#include <algorithm>
#include <cmath>
#include <iostream>
#include <vector>

template <typename Func, typename T>
T riemann(const Func& f, T x_min, T x_max, size_t N) {
  T dx = (x_max - x_min) / N;
  T sum = 0;

  for (T x = x_min; x < x_max; x += dx) {
    sum += f(x);
  }

  return sum * dx;
}

int main(int argc, char** argv) {
  std::vector<std::string> arguments(argv + 1, argv + argc);
  std::vector<double> coefficients(arguments.size());

  std::transform(arguments.begin(), arguments.end(), coefficients.begin(),
      [](std::string s) -> double { return std::stod(s); });

  auto sin = [](double x) { return std::sin(x); };

  std::cout << "integral of sin(x) from 0 to pi:" << std::endl;
  std::cout << "N=10: " << riemann(sin, 0., M_PI, 10) << std::endl;
  std::cout << "N=100: " << riemann(sin, 0., M_PI, 100) << std::endl;
  std::cout << "N=1000: " << riemann(sin, 0., M_PI, 1000) << std::endl;
  std::cout << std::endl;

  auto polynomial = [coefficients](double x) {
    double sum = 0;

    for (size_t i = 0; i < coefficients.size(); i++) {
      sum += coefficients[i] * std::pow(x, i);
    }

    return sum;
  };

  std::cout << "integral of polynomial from 0 to 1:" << std::endl;
  std::cout << "N=10: " << riemann(polynomial, 0., 1., 10) << std::endl;
  std::cout << "N=100: " << riemann(polynomial, 0., 1., 10) << std::endl;
  std::cout << "N=1000: " << riemann(polynomial, 0., 1., 10) << std::endl;
}