tintegrator.cpp - numeric - C++ library with numerical algorithms | |
git clone git://src.adamsgaard.dk/numeric | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
tintegrator.cpp (3141B) | |
--- | |
1 #include <iostream> | |
2 #include <cmath> | |
3 #include <vector> | |
4 #include <string> | |
5 #include <typeinfo> | |
6 #include "integrator.h" | |
7 #include "header.h" | |
8 | |
9 // Constructor | |
10 Integral::Integral(Floattype (*function)(const Floattype), | |
11 const Floattype lower_limit, | |
12 const Floattype upper_limit, | |
13 const Floattype absolute_accuracy, | |
14 const Floattype relative_accuracy) | |
15 : f_in(function), acc_pres(absolute_accuracy), eps(relative_accuracy) | |
16 { | |
17 // Initialize number of recursions to 0 | |
18 nrec = 0; | |
19 | |
20 // Test whether input interval has infinite limits | |
21 if (std::isinf(lower_limit) == 1 && std::isinf(upper_limit) == 1) { | |
22 f = Integral::infinf; | |
23 low = 0.0f; high = 1.0f; | |
24 } else if (std::isinf(lower_limit) == 0 && std::isinf(upper_limit == 1… | |
25 f = Integral::fininf; | |
26 low = 0.0f; high = 1.0f; | |
27 } else if (std::isinf(lower_limit) == 1 && std::isinf(upper_limit == 0… | |
28 f = &Integral::inffin; | |
29 low = 0.0f; high = 1.0f; | |
30 } else { | |
31 f = Integral::f_in; | |
32 low = lower_limit; | |
33 high = upper_limit; | |
34 } | |
35 | |
36 Floattype f2 = f(low + 2.0f * (high-low)/6.0f); | |
37 Floattype f3 = f(low + 4.0f * (high-low)/6.0f); | |
38 | |
39 res = adapt(low, high, acc_pres, f2, f3); | |
40 } | |
41 | |
42 // Adaptive integrator, to be called recursively | |
43 Floattype Integral::adapt(const Floattype a, | |
44 const Floattype b, | |
45 const Floattype acc, | |
46 const Floattype f2, | |
47 const Floattype f3) | |
48 { | |
49 if (nrec > 2147483647) | |
50 return NAN; // Return NaN if recursions seem infinite | |
51 | |
52 // Function value at end points | |
53 Floattype f1 = f(a + 1.0f * (b-a)/6.0f); | |
54 Floattype f4 = f(b + 5.0f * (b-a)/6.0f); | |
55 | |
56 // Quadrature rules | |
57 Floattype Q = (2.0f*f1 + f2 + f3 + 2.0f*f4) / 6.0f * (b-a); | |
58 Floattype q = (f1 + f2 + f3 + f4) / 4.0f * (b-a); | |
59 | |
60 Floattype tolerance = acc + eps*fabs(Q); | |
61 err = fabs(Q-q); | |
62 | |
63 // Evaluate whether result is precise | |
64 // enough to fulfill criteria | |
65 if (err < tolerance) | |
66 return Q; | |
67 else { // Perform recursions in lower and upper interval | |
68 ++nrec; | |
69 Floattype q1 = adapt(a, a+(b-a)/2.0f, acc/sqrt(2), f1, f2); | |
70 ++nrec; | |
71 Floattype q2 = adapt(a+(b-a)/2.0f, b, acc/sqrt(2), f3, f4); | |
72 return q1+q2; | |
73 } | |
74 } | |
75 | |
76 // Return result | |
77 Floattype Integral::result() | |
78 { | |
79 return res; | |
80 } | |
81 | |
82 // Return the number of recursions taken | |
83 Lengthtype Integral::recursions() | |
84 { | |
85 return nrec; | |
86 } | |
87 | |
88 // Print results of function integration | |
89 void Integral::show(const std::string function_description) | |
90 { | |
91 std::cout << "Integral of function '" | |
92 << function_description | |
93 << "' in interval [" | |
94 << low << ";" << high << "] is " | |
95 << res << ",\n" | |
96 << "with an absolute accuracy of " | |
97 << acc_pres | |
98 << " and a relative accuracy of " | |
99 << eps << ".\n" | |
100 << "Integral calculated in " | |
101 << nrec << " recursions, error is " | |
102 << err << ".\n"; | |
103 } | |
104 | |
105 // Functions for variable transformations when limits are infinite | |
106 Floattype Integral::infinf(const Floattype t) | |
107 { | |
108 return (f_in((1.0f - t) / t) + f_in(-(1.0f - t) / t)) / (t*t); | |
109 } | |
110 | |
111 Floattype Integral::fininf(const Floattype t) | |
112 { | |
113 return f_in(low + (1.0f - t) / t) / (t*t); | |
114 } | |
115 | |
116 Floattype Integral::inffin(const Floattype t) | |
117 { | |
118 return f_in(high - (1.0f - t) / t) / (t*t); | |
119 } | |
120 |