tintegrator.cpp - numeric - C++ library with numerical algorithms
(HTM) git clone git://src.adamsgaard.dk/numeric
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) 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