Introduction
Introduction Statistics Contact Development Disclaimer Help
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
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.