1 files changed, 59 insertions, 0 deletions
diff --git a/src/boost/libs/math/example/double_exponential.cpp b/src/boost/libs/math/example/double_exponential.cpp
new file mode 100644
index 00000000..87e21b18
--- /dev/null
+++ b/src/boost/libs/math/example/double_exponential.cpp
@@ -0,0 +1,59 @@
+// Copyright Nick Thompson, 2017
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <iostream>
+#include <cmath>
+#include <boost/math/quadrature/tanh_sinh.hpp>
+#include <boost/math/quadrature/sinh_sinh.hpp>
+#include <boost/math/quadrature/exp_sinh.hpp>
+
+using boost::math::quadrature::tanh_sinh;
+using boost::math::quadrature::sinh_sinh;
+using boost::math::quadrature::exp_sinh;
+using boost::math::constants::pi;
+using boost::math::constants::half_pi;
+using boost::math::constants::half;
+using boost::math::constants::third;
+using boost::math::constants::root_pi;
+using std::log;
+using std::cos;
+using std::cosh;
+using std::exp;
+using std::sqrt;
+
+int main()
+{
+    std::cout << std::setprecision(std::numeric_limits<double>::digits10);
+    double tol = sqrt(std::numeric_limits<double>::epsilon());
+    // For an integral over a finite domain, use tanh_sinh:
+    tanh_sinh<double> tanh_integrator(tol, 10);
+    auto f1 = [](double x) { return log(x)*log(1-x); };
+    double Q = tanh_integrator.integrate(f1, (double) 0, (double) 1);
+    double Q_expected = 2 - pi<double>()*pi<double>()*half<double>()*third<double>();
+
+    std::cout << "tanh_sinh quadrature of log(x)log(1-x) gives " << Q << std::endl;
+    std::cout << "The exact integral is                        " << Q_expected << std::endl;
+
+    // For an integral over the entire real line, use sinh-sinh quadrature:
+    sinh_sinh<double> sinh_integrator(tol, 10);
+    auto f2 = [](double t) { return cos(t)/cosh(t);};
+    Q = sinh_integrator.integrate(f2);
+    Q_expected = pi<double>()/cosh(half_pi<double>());
+    std::cout << "sinh_sinh quadrature of cos(x)/cosh(x) gives " << Q << std::endl;
+    std::cout << "The exact integral is                        " << Q_expected << std::endl;
+
+    // For half-infinite intervals, use exp-sinh.
+    // Endpoint singularities are handled well:
+    exp_sinh<double> exp_integrator(tol, 10);
+    auto f3 = [](double t) { return exp(-t)/sqrt(t); };
+    Q = exp_integrator.integrate(f3, 0, std::numeric_limits<double>::infinity());
+    Q_expected = root_pi<double>();
+    std::cout << "exp_sinh quadrature of exp(-t)/sqrt(t) gives " << Q << std::endl;
+    std::cout << "The exact integral is                        " << Q_expected << std::endl;
+
+
+
+}