1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
|
// Copyright (C) 2013 Steve Taylor (steve98654@gmail.com)
// License: Boost Software License See LICENSE.txt for full license
#ifndef DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSONh_
#define DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSONh_
#include "integrate_function_adapt_simpson_abstract.h"
#include "../assert.h"
// ----------------------------------------------------------------------------------------
namespace dlib
{
template <typename T, typename funct>
T impl_adapt_simp_stop(const funct& f, T a, T b, T fa, T fm, T fb, T is, int cnt)
{
const int maxint = 500;
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 ((std::abs(i1-i2) <= std::abs(is)) || (m <= a) || (b <= m))
{
Q = i1;
}
else
{
if(cnt < maxint)
{
cnt = cnt + 1;
Q = impl_adapt_simp_stop(f,a,m,fa,fml,fm,is,cnt)
+ impl_adapt_simp_stop(f,m,b,fm,fmr,fb,is,cnt);
}
}
return Q;
}
// ----------------------------------------------------------------------------------------
template <typename T, typename funct>
T integrate_function_adapt_simp(
const funct& f,
T a,
T b,
T tol = 1e-10
)
{
// make sure requires clause is not broken
DLIB_ASSERT(b > a && tol > 0,
"\t T integrate_function_adapt_simp()"
<< "\n\t Invalid arguments were given to this function."
<< "\n\t a: " << a
<< "\n\t b: " << b
<< "\n\t tol: " << tol
);
T eps = std::numeric_limits<T>::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;
int cnt = 0;
return impl_adapt_simp_stop(f, a, b, fa, fm, fb, is, cnt);
}
}
// ----------------------------------------------------------------------------------------
#endif // DLIB_INTEGRATE_FUNCTION_ADAPT_SIMPSONh_
|