1 // Copyright (C) 2013 Steve Taylor (steve98654@gmail.com) 2 // License: Boost Software License See LICENSE.txt for the full license. 3 4 // This function test battery is given in: 5 // 6 // Test functions taken from Pedro Gonnet's dissertation at ETH: 7 // Adaptive Quadrature Re-Revisited 8 // http://e-collection.library.ethz.ch/eserv/eth:65/eth-65-02.pdf 9 10 #include <math.h> 11 #include <dlib/matrix.h> 12 #include <dlib/numeric_constants.h> 13 #include <dlib/numerical_integration.h> 14 #include "tester.h" 15 16 namespace 17 { 18 using namespace test; 19 using namespace dlib; 20 using namespace std; 21 22 logger dlog("test.numerical_integration"); 23 24 class numerical_integration_tester : public tester 25 { 26 public: numerical_integration_tester()27 numerical_integration_tester ( 28 ) : 29 tester ("test_numerical_integration", 30 "Runs tests on the numerical integration function.", 31 0 32 ) 33 {} 34 perform_test()35 void perform_test() 36 { 37 38 dlog <<dlib::LINFO << "Testing integrate_function_adapt_simpson"; 39 40 matrix<double,23,1> m; 41 double tol = 1e-10; 42 double eps = 1e-8; 43 44 m(0) = integrate_function_adapt_simp(&gg1, 0.0, 1.0, tol); 45 m(1) = integrate_function_adapt_simp(&gg2, 0.0, 1.0, tol); 46 m(2) = integrate_function_adapt_simp(&gg3, 0.0, 1.0, tol); 47 m(3) = integrate_function_adapt_simp(&gg4, 0.0, 1.0, tol); 48 m(4) = integrate_function_adapt_simp(&gg5, -1.0, 1.0, tol); 49 m(5) = integrate_function_adapt_simp(&gg6, 0.0, 1.0, tol); 50 m(6) = integrate_function_adapt_simp(&gg7, 0.0, 1.0, tol); 51 m(7) = integrate_function_adapt_simp(&gg8, 0.0, 1.0, tol); 52 m(8) = integrate_function_adapt_simp(&gg9, 0.0, 1.0, tol); 53 m(9) = integrate_function_adapt_simp(&gg10, 0.0, 1.0, tol); 54 m(10) = integrate_function_adapt_simp(&gg11, 0.0, 1.0, tol); 55 m(11) = integrate_function_adapt_simp(&gg12, 1e-6, 1.0, tol); 56 m(12) = integrate_function_adapt_simp(&gg13, 0.0, 10.0, tol); 57 m(13) = integrate_function_adapt_simp(&gg14, 0.0, 10.0, tol); 58 m(14) = integrate_function_adapt_simp(&gg15, 0.0, 10.0, tol); 59 m(15) = integrate_function_adapt_simp(&gg16, 0.01, 1.0, tol); 60 m(16) = integrate_function_adapt_simp(&gg17, 0.0, pi, tol); 61 m(17) = integrate_function_adapt_simp(&gg18, 0.0, 1.0, tol); 62 m(18) = integrate_function_adapt_simp(&gg19, -1.0, 1.0, tol); 63 m(19) = integrate_function_adapt_simp(&gg20, 0.0, 1.0, tol); 64 m(20) = integrate_function_adapt_simp(&gg21, 0.0, 1.0, tol); 65 m(21) = integrate_function_adapt_simp(&gg22, 0.0, 5.0, tol); 66 67 // Here we compare the approximated integrals against 68 // highly accurate approximations generated either from 69 // the exact integral values or Mathematica's NIntegrate 70 // function using a working precision of 20. 71 72 DLIB_TEST(abs(m(0) - 1.7182818284590452354) < 1e-11); 73 DLIB_TEST(abs(m(1) - 0.7000000000000000000) < eps); 74 DLIB_TEST(abs(m(2) - 0.6666666666666666667) < eps); 75 DLIB_TEST(abs(m(3) - 0.2397141133444008336) < eps); 76 DLIB_TEST(abs(m(4) - 1.5822329637296729331) < 1e-11); 77 DLIB_TEST(abs(m(5) - 0.4000000000000000000) < eps); 78 DLIB_TEST(abs(m(6) - 2.0000000000000000000) < 1e-4); 79 DLIB_TEST(abs(m(7) - 0.8669729873399110375) < eps); 80 DLIB_TEST(abs(m(8) - 1.1547005383792515290) < eps); 81 DLIB_TEST(abs(m(9) - 0.6931471805599453094) < eps); 82 DLIB_TEST(abs(m(10) - 0.3798854930417224753) < eps); 83 DLIB_TEST(abs(m(11) - 0.7775036341124982763) < eps); 84 DLIB_TEST(abs(m(12) - 0.5000000000000000000) < eps); 85 DLIB_TEST(abs(m(13) - 1.0000000000000000000) < eps); 86 DLIB_TEST(abs(m(14) - 0.4993633810764567446) < eps); 87 DLIB_TEST(abs(m(15) - 0.1121393035410217 ) < eps); 88 DLIB_TEST(abs(m(16) - 0.2910187828600526985) < eps); 89 DLIB_TEST(abs(m(17) + 0.4342944819032518276) < 1e-5); 90 DLIB_TEST(abs(m(18) - 1.56439644406905 ) < eps); 91 DLIB_TEST(abs(m(19) - 0.1634949430186372261) < eps); 92 DLIB_TEST(abs(m(20) - 0.0134924856494677726) < eps); 93 } 94 gg1(double x)95 static double gg1(double x) 96 { 97 return pow(e,x); 98 } 99 gg2(double x)100 static double gg2(double x) 101 { 102 if(x > 0.3) 103 { 104 return 1.0; 105 } 106 else 107 { 108 return 0; 109 } 110 } 111 gg3(double x)112 static double gg3(double x) 113 { 114 return pow(x,0.5); 115 } 116 gg4(double x)117 static double gg4(double x) 118 { 119 return 23.0/25.0*cosh(x)-cos(x); 120 } 121 gg5(double x)122 static double gg5(double x) 123 { 124 return 1/(pow(x,4) + pow(x,2) + 0.9); 125 } 126 gg6(double x)127 static double gg6(double x) 128 { 129 return pow(x,1.5); 130 } 131 gg7(double x)132 static double gg7(double x) 133 { 134 return pow(x,-0.5); 135 } 136 gg8(double x)137 static double gg8(double x) 138 { 139 return 1/(1 + pow(x,4)); 140 } 141 gg9(double x)142 static double gg9(double x) 143 { 144 return 2/(2 + sin(10*pi*x)); 145 } 146 gg10(double x)147 static double gg10(double x) 148 { 149 return 1/(1+x); 150 } 151 gg11(double x)152 static double gg11(double x) 153 { 154 return 1.0/(1 + pow(e,x)); 155 } 156 gg12(double x)157 static double gg12(double x) 158 { 159 return x/(pow(e,x)-1.0); 160 } 161 gg13(double x)162 static double gg13(double x) 163 { 164 return sqrt(50.0)*pow(e,-50.0*pi*x*x); 165 } 166 gg14(double x)167 static double gg14(double x) 168 { 169 return 25.0*pow(e,-25.0*x); 170 } 171 gg15(double x)172 static double gg15(double x) 173 { 174 return 50.0/(pi*(2500.0*x*x+1)); 175 } 176 gg16(double x)177 static double gg16(double x) 178 { 179 return 50.0*pow((sin(50.0*pi*x)/(50.0*pi*x)),2); 180 } 181 gg17(double x)182 static double gg17(double x) 183 { 184 return cos(cos(x)+3*sin(x)+2*cos(2*x)+3*cos(3*x)); 185 } 186 gg18(double x)187 static double gg18(double x) 188 { 189 return log10(x); 190 } 191 gg19(double x)192 static double gg19(double x) 193 { 194 return 1/(1.005+x*x); 195 } 196 gg20(double x)197 static double gg20(double x) 198 { 199 return 1/cosh(20.0*(x-1.0/5.0)) + 1/cosh(400.0*(x-2.0/5.0)) 200 + 1/cosh(8000.0*(x-3.0/5.0)); 201 } 202 gg21(double x)203 static double gg21(double x) 204 { 205 return 1.0/(1.0+(230.0*x-30.0)*(230.0*x-30.0)); 206 } 207 gg22(double x)208 static double gg22(double x) 209 { 210 if(x < 1) 211 { 212 return (x + 1.0); 213 } 214 else if(x >= 1 && x <= 3) 215 { 216 return (3.0 - x); 217 } 218 else 219 { 220 return 2.0; 221 } 222 } 223 224 }; 225 226 numerical_integration_tester a; 227 } 228 229