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
94
95
96
97
98
|
#include "integration.h"
#include <tjutils/tjtest.h>
#ifdef HAVE_LIBGSL
#include <gsl/gsl_integration.h>
double Integrand::get_integral(double xmin, double xmax, unsigned int max_subintervals, double error_limit) const {
FunctionIntegral fi(*this,max_subintervals,error_limit);
return fi.get_integral(xmin,xmax);
}
//////////////////////////////////////////////////////////////
struct GslData4Integr {
gsl_integration_workspace* w;
};
//////////////////////////////////////////////////////////////
FunctionIntegral::FunctionIntegral(const Integrand& func, unsigned int max_subintervals, double error_limit)
: f(func), n_intervals(max_subintervals), errlimit(error_limit) {
gsldata=new GslData4Integr;
gsldata->w = gsl_integration_workspace_alloc(n_intervals);
}
FunctionIntegral::~FunctionIntegral() {
gsl_integration_workspace_free(gsldata->w);
delete gsldata;
}
double FunctionIntegral::get_integral(double xmin, double xmax) const {
double result, error;
gsl_function gsl_F;
gsl_F.function = &integrand;
gsl_F.params = (void*)&f;
gsl_integration_qags (&gsl_F, xmin, xmax, 0, errlimit, n_intervals, gsldata->w, &result, &error);
return result;
}
double FunctionIntegral::integrand(double x, void *params) {
const Integrand* integr=(const Integrand*)params;
return integr->evaluate(x);
}
//////////////////////////////////////////////////////////////
// Unit test
#ifndef NO_UNIT_TEST
class FunctionIntegralTest : public UnitTest {
struct QuadrFunction : public Integrand {
double evaluate(double x) const {return x*x;}
};
struct QuadrIntFunction : public Integrand {
double evaluate(double x) const {return qf.get_integral(0.0,x);}
QuadrFunction qf;
};
public:
FunctionIntegralTest() : UnitTest("FunctionIntegral") {}
private:
bool check() const {
Log<UnitTest> odinlog(this,"check");
// testing quadratic integral function
QuadrIntFunction func;
STD_string expected_str=ftos(1.0/12.0);
STD_string calculated_str=ftos(func.get_integral(0.0,1.0));
if(calculated_str!=expected_str) {
ODINLOG(odinlog,errorLog) << "integral=" << calculated_str << ", but expected integral=" << expected_str << STD_endl;
return false;
}
return true;
}
};
void alloc_FunctionIntegralTest() {new FunctionIntegralTest();} // create test instance
#endif
#else
#error "GNU Scientific library is missing!"
#endif
|
