1 // Copyright Nick Thompson, 2017
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt
5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
6
7 #include <iostream>
8 #include <cmath>
9 #include <boost/math/quadrature/tanh_sinh.hpp>
10 #include <boost/math/quadrature/sinh_sinh.hpp>
11 #include <boost/math/quadrature/exp_sinh.hpp>
12
13 using boost::math::quadrature::tanh_sinh;
14 using boost::math::quadrature::sinh_sinh;
15 using boost::math::quadrature::exp_sinh;
16 using boost::math::constants::pi;
17 using boost::math::constants::half_pi;
18 using boost::math::constants::half;
19 using boost::math::constants::third;
20 using boost::math::constants::root_pi;
21 using std::log;
22 using std::cos;
23 using std::cosh;
24 using std::exp;
25 using std::sqrt;
26
main()27 int main()
28 {
29 std::cout << std::setprecision(std::numeric_limits<double>::digits10);
30 double tol = sqrt(std::numeric_limits<double>::epsilon());
31 // For an integral over a finite domain, use tanh_sinh:
32 tanh_sinh<double> tanh_integrator(tol, 10);
33 auto f1 = [](double x) { return log(x)*log(1-x); };
34 double Q = tanh_integrator.integrate(f1, (double) 0, (double) 1);
35 double Q_expected = 2 - pi<double>()*pi<double>()*half<double>()*third<double>();
36
37 std::cout << "tanh_sinh quadrature of log(x)log(1-x) gives " << Q << std::endl;
38 std::cout << "The exact integral is " << Q_expected << std::endl;
39
40 // For an integral over the entire real line, use sinh-sinh quadrature:
41 sinh_sinh<double> sinh_integrator(tol, 10);
42 auto f2 = [](double t) { return cos(t)/cosh(t);};
43 Q = sinh_integrator.integrate(f2);
44 Q_expected = pi<double>()/cosh(half_pi<double>());
45 std::cout << "sinh_sinh quadrature of cos(x)/cosh(x) gives " << Q << std::endl;
46 std::cout << "The exact integral is " << Q_expected << std::endl;
47
48 // For half-infinite intervals, use exp-sinh.
49 // Endpoint singularities are handled well:
50 exp_sinh<double> exp_integrator(tol, 10);
51 auto f3 = [](double t) { return exp(-t)/sqrt(t); };
52 Q = exp_integrator.integrate(f3, 0, std::numeric_limits<double>::infinity());
53 Q_expected = root_pi<double>();
54 std::cout << "exp_sinh quadrature of exp(-t)/sqrt(t) gives " << Q << std::endl;
55 std::cout << "The exact integral is " << Q_expected << std::endl;
56
57
58
59 }
60