New answers tagged c
3
If you can use Boost, you'll get the benefit of exp-sinh quadrature, which appears to work well on this problem:
#include <iostream>
#include <cmath>
#include <boost/math/quadrature/exp_sinh.hpp>
using boost::math::quadrature::exp_sinh;
using std::exp;
using std::expm1;
using std::log;
int main() {
exp_sinh<double> integrator;
...
3
(I don't have any experience with Cubature, I don't know which algorithms they implement.)
If the integration rule you use has negative weights, like Newton-Cotes rules, then even if the integrand is strictly positive, the result you get may be negative. This was an issue we would commonly encounter when using MFEM to compute norms (probably other FEM ...
7
This might be an accuracy problem in computing the second term, because of those large exponentials when $x \gg 1$.
I would first work on that term: gather $e^x$ out from numerator and denominator and simplify it out (leaving a summand $e^{-x}$ in the numerator).
However, rather than guessing it might be necessary to understand better where the accuracy ...
Top 50 recent answers are included