Implemented a numerical quadrature method based on an adaptive

Simpson rule.  Added unit tests and supporting examples for this
function.

dlib/integrate_function_adapt_simpson.h

+// 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

dlib/numerical_integration/integrate_function_adapt_simpson.h

+// 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 <typename T, typename funct>
+T integrate_function_adapt_simp(const funct& f, T a, T b, T tol)
+{
+    T eps = std::numeric_limits<double>::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 <typename T, typename funct>
+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__
+

dlib/numerical_integration/integrate_function_adapt_simpson_abstract.h

+// 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 <typename T, typename funct>
+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__
+

dlib/test/numerical_integration.cpp

+// 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 <math.h>
+#include <dlib/matrix.h>
+#include <dlib/numeric_constants.h>
+#include <dlib/integrate_function_adapt_simpson.h>
+
+#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 <<dlib::LINFO << "Testing integrate_function_adapt_simpson";
+
+            matrix<double,23,1> 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;
+}
+

examples/integrate_function_adapt_simp_ex.cpp

+// 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 <iostream>
+#include <stdint.h>
+#include <dlib/matrix.h>
+#include <dlib/numeric_constants.h>
+#include <dlib/integrate_function_adapt_simpson.h>
+
+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<double,5,1> 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;
+}
+