diff --git a/dlib/integrate_function_adapt_simpson.h b/dlib/integrate_function_adapt_simpson.h new file mode 100644 index 000000000..5110b279e --- /dev/null +++ b/dlib/integrate_function_adapt_simpson.h @@ -0,0 +1,9 @@ +// Copyright (C) 2013 Steve Taylor (steve98654@gmail.com) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_HEADER +#define DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_HEADER + +#include "matrix.h" +#include "numerical_integration/integrate_function_adapt_simpson.h" + +#endif // DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_HEADER diff --git a/dlib/numerical_integration/integrate_function_adapt_simpson.h b/dlib/numerical_integration/integrate_function_adapt_simpson.h new file mode 100644 index 000000000..c8602fbbb --- /dev/null +++ b/dlib/numerical_integration/integrate_function_adapt_simpson.h @@ -0,0 +1,88 @@ +// Copyright (C) 2013 Steve Taylor (steve98654@gmail.com) +// License: Boost Software License See LICENSE.txt for full license + +#ifndef DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON__ +#define DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON__ + +#include "../matrix.h" + +class adapt_simp +{ +public: + +template +T integrate_function_adapt_simp(const funct& f, T a, T b, T tol) +{ + T eps = std::numeric_limits::epsilon(); + + if(tol < eps) + { + tol = eps; + } + + const T ba = b-a; + const T fa = f(a); + const T fb = f(b); + const T fm = f((a+b)/2); + + T is =ba/8*(fa+fb+fm+ f(a + 0.9501*ba) + f(a + 0.2311*ba) + f(a + 0.6068*ba) + + f(a + 0.4860*ba) + f(a + 0.8913*ba)); + + if(is == 0) + { + is = b-a; + } + + is = is*tol/eps; + + int cnt = 0; + + T tstvl = adapt_simp_stop(f, a, b, fa, fm, fb, is, cnt); + + return tstvl; + +} + +private: + +template +T adapt_simp_stop(const funct& f, T a, T b, T fa, T fm, T fb, T is, int &cnt) +{ + int MAXINT = 1000; + + T m = (a + b)/2.0; + T h = (b - a)/4.0; + T fml = f(a + h); + T fmr = f(b - h); + T i1 = h/1.5*(fa+4.0*fm+fb); + T i2 = h/3.0*(fa+4.0*(fml+fmr)+2.0*fm+fb); + i1 = (16.0*i2 - i1)/15.0; + T Q = 0; + + if((is+(i1-i2) == is) || (m <= a) || (b <= m)) + { + if((m <= a) || (b <= m)) + { + DLIB_ASSERT(1, "\tintegrate_function_adapt_simpson::adapt_simp_stop" + << "\n\tmidpoint evaluation occurred at endpoint"); + } + + Q = i1; + } + else + { + if(cnt < MAXINT) + {cnt = cnt + 1; + + Q = adapt_simp_stop(f,a,m,fa,fml,fm,is,cnt) + + adapt_simp_stop(f,m,b,fm,fmr,fb,is,cnt); + } + } + + return Q; +} + +}; + +#endif //DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON.h__ + diff --git a/dlib/numerical_integration/integrate_function_adapt_simpson_abstract.h b/dlib/numerical_integration/integrate_function_adapt_simpson_abstract.h new file mode 100644 index 000000000..48e539dfc --- /dev/null +++ b/dlib/numerical_integration/integrate_function_adapt_simpson_abstract.h @@ -0,0 +1,24 @@ +// Copyright (C) 2013 Steve Taylor (steve98654@gmail.com) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_ABSTRACT__ +#define DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON_ABSTRACT__ + +#include "matrix.h" + +templte +T integrate_function_adapt_simp(const funct& f, T a, T b, T tol); +/*! + requires + - b > a + - tol > 0 + - f to be real valued single variable function + + ensures + - returns an approximation of the integral of f over the domain [a,b] + using the adaptive Simpson method outlined in + ander, W. and W. Gautshi, "Adaptive Quadrature -- Revisited" + BIT, Vol. 40, (2000), pp.84-101 +!*/ + +#endif // DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSON__ + diff --git a/dlib/test/numerical_integration.cpp b/dlib/test/numerical_integration.cpp new file mode 100644 index 000000000..06a169b04 --- /dev/null +++ b/dlib/test/numerical_integration.cpp @@ -0,0 +1,231 @@ +// Copyright (C) 2013 Steve Taylor (steve98654@gmail.com) +// License: Boost Software License See LICENSE.txt for the full license. + +// This function test battery is given in: +// +// Test functions taken from Pedro Gonnet's dissertation at ETH: +// Adaptive Quadrature Re-Revisited +// http://e-collection.library.ethz.ch/eserv/eth:65/eth-65-02.pdf + +#include +#include +#include +#include + +#include "tester.h" + +namespace +{ + using namespace test; + using namespace dlib; + using namespace std; + + logger dlog("test.numerical_integration"); + + class numerical_integration_tester : public tester + { + public: + numerical_integration_tester ( + ) : + tester ("test_numerical_integration", + "Runs tests on the numerical integration function.", + 0 + ) + {} + + void perform_test() + { + + dlog < m; + double tol = 1e-10; + double eps = 1e-4; + + adapt_simp ad; + + m(0) = ad.integrate_function_adapt_simp(&gg1, 0.0, 1.0, tol); + m(1) = ad.integrate_function_adapt_simp(&gg2, 0.0, 1.0, tol); + m(2) = ad.integrate_function_adapt_simp(&gg3, 0.0, 1.0, tol); + m(3) = ad.integrate_function_adapt_simp(&gg4, 0.0, 1.0, tol); + m(4) = ad.integrate_function_adapt_simp(&gg5, -1.0, 1.0, tol); + m(5) = ad.integrate_function_adapt_simp(&gg6, 0.0, 1.0, tol); + m(6) = ad.integrate_function_adapt_simp(&gg7, 0.0, 1.0, tol); + m(7) = ad.integrate_function_adapt_simp(&gg8, 0.0, 1.0, tol); + m(8) = ad.integrate_function_adapt_simp(&gg9, 0.0, 1.0, tol); + m(9) = ad.integrate_function_adapt_simp(&gg10, 0.0, 1.0, tol); + m(10) = ad.integrate_function_adapt_simp(&gg11, 0.0, 1.0, tol); + m(11) = ad.integrate_function_adapt_simp(&gg12, 1e-6, 1.0, tol); + m(12) = ad.integrate_function_adapt_simp(&gg13, 0.0, 10.0, tol); + m(13) = ad.integrate_function_adapt_simp(&gg14, 0.0, 10.0, tol); + m(14) = ad.integrate_function_adapt_simp(&gg15, 0.0, 10.0, tol); + m(15) = ad.integrate_function_adapt_simp(&gg16, 0.01, 1.0, tol); + m(16) = ad.integrate_function_adapt_simp(&gg17, 0.0, pi, tol); + m(17) = ad.integrate_function_adapt_simp(&gg18, 0.0, 1.0, tol); + m(18) = ad.integrate_function_adapt_simp(&gg19, -1.0, 1.0, tol); + m(19) = ad.integrate_function_adapt_simp(&gg20, 0.0, 1.0, tol); + m(20) = ad.integrate_function_adapt_simp(&gg21, 0.0, 1.0, tol); + m(21) = ad.integrate_function_adapt_simp(&gg22, 0.0, 5.0, tol); + + // Here we compare the approximated integrals against + // highly accurate approximations generated either from + // the exact integral values or Mathematica's NIntegrate + // function using a working precision of 20. + + DLIB_TEST(abs(m(0) - 1.7182818284590452354) < eps); + DLIB_TEST(abs(m(1) - 0.7000000000000000000) < eps); + DLIB_TEST(abs(m(2) - 0.6666666666666666667) < eps); + DLIB_TEST(abs(m(3) - 0.2397141133444008336) < eps); + DLIB_TEST(abs(m(4) - 1.5822329637296729331) < eps); + DLIB_TEST(abs(m(5) - 0.4000000000000000000) < eps); + DLIB_TEST(abs(m(6) - 2.0000000000000000000) < eps); + DLIB_TEST(abs(m(7) - 0.8669729873399110375) < eps); + DLIB_TEST(abs(m(8) - 1.1547005383792515290) < eps); + DLIB_TEST(abs(m(9) - 0.6931471805599453094) < eps); + DLIB_TEST(abs(m(10) - 0.3798854930417224753) < eps); + DLIB_TEST(abs(m(11) - 0.7775036341124982763) < eps); + DLIB_TEST(abs(m(12) - 0.5000000000000000000) < eps); + DLIB_TEST(abs(m(13) - 1.0000000000000000000) < eps); + DLIB_TEST(abs(m(14) - 0.4993633810764567446) < eps); + DLIB_TEST(abs(m(15) - 0.1121393035410217 ) < eps); + DLIB_TEST(abs(m(16) - 0.2910187828600526985) < eps); + DLIB_TEST(abs(m(17) + 0.4342944819032518276) < eps); + DLIB_TEST(abs(m(18) - 1.56439644406905 ) < eps); + DLIB_TEST(abs(m(19) - 0.1634949430186372261) < eps); + DLIB_TEST(abs(m(20) - 0.0134924856494677726) < eps); + } + + static double gg1(double x) + { + return pow(e,x); + } + + static double gg2(double x) + { + if(x > 0.3) + { + return 1.0; + } + else + { + return 0; + } + } + + static double gg3(double x) + { + return pow(x,0.5); + } + + static double gg4(double x) + { + return 23.0/25.0*cosh(x)-cos(x); + } + + static double gg5(double x) + { + return 1/(pow(x,4) + pow(x,2) + 0.9); + } + + static double gg6(double x) + { + return pow(x,1.5); + } + + static double gg7(double x) + { + return pow(x,-0.5); + } + + static double gg8(double x) + { + return 1/(1 + pow(x,4)); + } + + static double gg9(double x) + { + return 2/(2 + sin(10*pi*x)); + } + + static double gg10(double x) + { + return 1/(1+x); + } + + static double gg11(double x) + { + return 1.0/(1 + pow(e,x)); + } + + static double gg12(double x) + { + return x/(pow(e,x)-1.0); + } + + static double gg13(double x) + { + return sqrt(50)*pow(e,-50.0*pi*x*x); + } + + static double gg14(double x) + { + return 25.0*pow(e,-25.0*x); + } + + static double gg15(double x) + { + return 50.0/(pi*(2500.0*x*x+1)); + } + + static double gg16(double x) + { + return 50.0*pow((sin(50.0*pi*x)/(50.0*pi*x)),2); + } + + static double gg17(double x) + { + return cos(cos(x)+3*sin(x)+2*cos(2*x)+3*cos(3*x)); + } + + static double gg18(double x) + { + return log10(x); + } + + static double gg19(double x) + { + return 1/(1.005+x*x); + } + + static double gg20(double x) + { + return 1/cosh(20.0*(x-1.0/5.0)) + 1/cosh(400.0*(x-2.0/5.0)) + + 1/cosh(8000.0*(x-3.0/5.0)); + } + + static double gg21(double x) + { + return 1.0/(1.0+(230.0*x-30.0)*(230.0*x-30.0)); + } + + static double gg22(double x) + { + if(x < 1) + { + return (x + 1.0); + } + else if(x >= 1 && x <= 3) + { + return (3.0 - x); + } + else + { + return 2.0; + } + } + + }; + + numerical_integration_tester a; +} + diff --git a/examples/integrate_function_adapt_simp_ex.cpp b/examples/integrate_function_adapt_simp_ex.cpp new file mode 100644 index 000000000..3de64b937 --- /dev/null +++ b/examples/integrate_function_adapt_simp_ex.cpp @@ -0,0 +1,106 @@ +// Copyright (C) 2013 Steve Taylor (steve98654@gmai.com) +// The contents of this file are in the public domain. See +// LICENSE_FOR_EXAMPLE_PROGRAMS.txt + +/* + + This example demonstrates the usage of the numerical quadrature function + integrate_function_adapt_simpson. This function takes as input a single variable + function, the endpoints of a domain over which the function will be integrated, and a + tolerance parameter. It outputs an approximation of the integral of this function + over the specified domain. The algorithm is based on the adaptive Simpson method outlined in: + + Numerical Integration method based on the adaptive Simpson method in + Gander, W. and W. Gautschi, "Adaptive Quadrature – Revisited," + BIT, Vol. 40, 2000, pp. 84-101 + +*/ + +#include +#include +#include +#include +#include + +using namespace std; +using namespace dlib; + +// Here we define a class that consists of the set of functions that we +// wish to integrate and comment in the domain of integration. + +class int_fcts { + +public: + + // x in [0,1] + static double gg1(double x) + { + return pow(e,x); + } + + // x in [0,1] + static double gg2(double x) + { + return x*x; + } + + // x in [0, pi] + static double gg3(double x) + { + return 1/(x*x + cos(x)*cos(x)); + } + + // x in [-pi, pi] + static double gg4(double x) + { + return sin(x); + } + + // x in [0,2] + static double gg5(double x) + { + return 1/(1 + x*x); + } +}; + +// Examples +int main() +{ + // We first define a matrix m to which we will store the approximated values of the + // integrals of our int_fcts functions. + + matrix m; + + // Next we define a tolerance parameter. Roughly speaking, a lower tolerance will + // result in a more accurate approximation of the true integral. However, there are + // instances where too small of a tolerance may yield a less accurate approximation + // than a larger tolerance. We recommend taking the tolerance to be in the + // [1e-10, 1e-8] region. + + double tol = 1e-10; + + // Here we instantiate a class which contains the numerical quadrature method. + adapt_simp ad; + + // Here we compute the integrals of the five functions defined above using the same + // tolerance level for each. + + m(0) = ad.integrate_function_adapt_simp(&int_fcts::gg1, 0.0, 1.0, tol); + m(1) = ad.integrate_function_adapt_simp(&int_fcts::gg2, 0.0, 1.0, tol); + m(2) = ad.integrate_function_adapt_simp(&int_fcts::gg3, 0.0, pi, tol); + m(3) = ad.integrate_function_adapt_simp(&int_fcts::gg4, -pi, pi, tol); + m(4) = ad.integrate_function_adapt_simp(&int_fcts::gg5, 0.0, 2.0, tol); + + // We finally print out the values of each of the approximated integrals to ten + // significant digits. + + cout << "\nThe integral of exp(x) for x in [0,1] is " << std::setprecision(10) << m(0) << endl; + cout << "The integral of x^2 for in [0,1] is " << std::setprecision(10) << m(1) << endl; + cout << "The integral of 1/(x^2 + cos(x)^2) for in [0,pi] is " << std::setprecision(10) << m(2) << endl; + cout << "The integral of sin(x) for in [-pi,pi] is " << std::setprecision(10) << m(3) << endl; + cout << "The integral of 1/(1+x^2) for in [0,2] is " << std::setprecision(10) << m(4) << endl; + cout << "" << endl; + + return 0; +} +