Integration.h

Go to the documentation of this file.

#pragma once
#include "../../Matrix.h"
#include "../../sorting.h"
#include "../utils.h"
#include <functional>
 
namespace integration {
  std::pair<double, std::vector<double>> quad_adaptive_rec(
  const std::function<double(double)>& fun,
  const double& a,
  const double& b,
  const double& tol,
  const double& hMin,
  std::vector<double>& nodes,
  Matrix<double>& values) {
    double h = b - a;
    if(h < hMin) {
      // error
      std::cerr << "H too small! h = " << h << std::endl;
      exit(-1);
    }
 
    Matrix<double> x = { { a }, { a + (h / 2.0) }, { b } };
    if(values.rows() == 0) {
      values = x.Apply(fun);
    } else if(values.rows() < 3) {
      values = { { values(0, 0) }, { fun(x(1, 0)) }, { values(1, 0) } };
    }
 
    auto qSimpson = (h / 6.0) * (values(2, 0) + 4.0 * values(1, 0) + values(0, 0));
    auto qTrapez  = (h / 2.0) * (values(2, 0) + values(0, 0));
 
    auto error = fabs(qSimpson - qTrapez);
 
    nodes.push_back(x(1, 0));
 
    if(error <= tol) { return { qSimpson, nodes }; }
 
    auto hTol        = tol / 2.0;
    auto leftValues  = values.GetSlice(0, 1, 0, 0);
    auto resLeft     = quad_adaptive_rec(fun, a, x(1, 0), hTol, hMin, nodes, leftValues);
    auto rightValues = values.GetSlice(1, 2, 0, 0);
    auto resRight    = quad_adaptive_rec(fun, x(1, 0), b, hTol, hMin, nodes, rightValues);
 
    return { resLeft.first + resRight.first, nodes };
  }
} // namespace integration
 
std::pair<double, std::vector<double>> quadrature(
const std::function<double(double)>& fun,
const double& lower,
const double& upper,
const double& tol,
const double& hMin) {
  std::vector<double> nodes = { lower, upper };
  Matrix<double> q;
  auto res   = integration::quad_adaptive_rec(fun, lower, upper, tol, hMin, nodes, q);
  res.second = sort(res.second);
  return res;
}