38#define G_ROMBERG "GIntegral::romberg(double&, double&, int&)"
39#define G_TRAPZD "GIntegral::trapzd(double&, double&, int&, double)"
40#define G_POLINT "GIntegral::polint(double*, double*, int, double, double*)"
53 0.973906528517171720077964012084452,
54 0.865063366688984510732096688423493,
55 0.679409568299024406234327365114874,
56 0.433395394129247190799265943165784,
57 0.148874338981631210884826001129720
62 0.066671344308688137593568809893332,
63 0.149451349150580593145776339657697,
64 0.219086362515982043995534934228163,
65 0.269266719309996355091226921569469,
66 0.295524224714752870173892994651338
71 0.995657163025808080735527280689003,
72 0.930157491355708226001207180059508,
73 0.780817726586416897063717578345042,
74 0.562757134668604683339000099272694,
75 0.294392862701460198131126603103866
80 0.032558162307964727478818972459390,
81 0.075039674810919952767043140916190,
82 0.109387158802297641899210590325805,
83 0.134709217311473325928054001771707,
84 0.147739104901338491374841515972068
89 0.011694638867371874278064396062192,
90 0.054755896574351996031381300244580,
91 0.093125454583697605535065465083366,
92 0.123491976262065851077958109831074,
93 0.142775938577060080797094273138717,
94 0.149445554002916905664936468389821
99 0.999333360901932081394099323919911,
100 0.987433402908088869795961478381209,
101 0.954807934814266299257919200290473,
102 0.900148695748328293625099494069092,
103 0.825198314983114150847066732588520,
104 0.732148388989304982612354848755461,
105 0.622847970537725238641159120344323,
106 0.499479574071056499952214885499755,
107 0.364901661346580768043989548502644,
108 0.222254919776601296498260928066212,
109 0.074650617461383322043914435796506
114 0.016296734289666564924281974617663,
115 0.037522876120869501461613795898115,
116 0.054694902058255442147212685465005,
117 0.067355414609478086075553166302174,
118 0.073870199632393953432140695251367,
119 0.005768556059769796184184327908655,
120 0.027371890593248842081276069289151,
121 0.046560826910428830743339154433824,
122 0.061744995201442564496240336030883,
123 0.071387267268693397768559114425516
128 0.001844477640212414100389106552965,
129 0.010798689585891651740465406741293,
130 0.021895363867795428102523123075149,
131 0.032597463975345689443882222526137,
132 0.042163137935191811847627924327955,
133 0.050741939600184577780189020092084,
134 0.058379395542619248375475369330206,
135 0.064746404951445885544689259517511,
136 0.069566197912356484528633315038405,
137 0.072824441471833208150939535192842,
138 0.074507751014175118273571813842889,
139 0.074722147517403005594425168280423
144 0.999902977262729234490529830591582,
145 0.997989895986678745427496322365960,
146 0.992175497860687222808523352251425,
147 0.981358163572712773571916941623894,
148 0.965057623858384619128284110607926,
149 0.943167613133670596816416634507426,
150 0.915806414685507209591826430720050,
151 0.883221657771316501372117548744163,
152 0.845710748462415666605902011504855,
153 0.803557658035230982788739474980964,
154 0.757005730685495558328942793432020,
155 0.706273209787321819824094274740840,
156 0.651589466501177922534422205016736,
157 0.593223374057961088875273770349144,
158 0.531493605970831932285268948562671,
159 0.466763623042022844871966781659270,
160 0.399424847859218804732101665817923,
161 0.329874877106188288265053371824597,
162 0.258503559202161551802280975429025,
163 0.185695396568346652015917141167606,
164 0.111842213179907468172398359241362,
165 0.037352123394619870814998165437704
170 0.008148377384149172900002878448190,
171 0.018761438201562822243935059003794,
172 0.027347451050052286161582829741283,
173 0.033677707311637930046581056957588,
174 0.036935099820427907614589586742499,
175 0.002884872430211530501334156248695,
176 0.013685946022712701888950035273128,
177 0.023280413502888311123409291030404,
178 0.030872497611713358675466394126442,
179 0.035693633639418770719351355457044,
180 0.000915283345202241360843392549948,
181 0.005399280219300471367738743391053,
182 0.010947679601118931134327826856808,
183 0.016298731696787335262665703223280,
184 0.021081568889203835112433060188190,
185 0.025370969769253827243467999831710,
186 0.029189697756475752501446154084920,
187 0.032373202467202789685788194889595,
188 0.034783098950365142750781997949596,
189 0.036412220731351787562801163687577,
190 0.037253875503047708539592001191226
195 0.000274145563762072350016527092881,
196 0.001807124155057942948341311753254,
197 0.004096869282759164864458070683480,
198 0.006758290051847378699816577897424,
199 0.009549957672201646536053581325377,
200 0.012329447652244853694626639963780,
201 0.015010447346388952376697286041943,
202 0.017548967986243191099665352925900,
203 0.019938037786440888202278192730714,
204 0.022194935961012286796332102959499,
205 0.024339147126000805470360647041454,
206 0.026374505414839207241503786552615,
207 0.028286910788771200659968002987960,
208 0.030052581128092695322521110347341,
209 0.031646751371439929404586051078883,
210 0.033050413419978503290785944862689,
211 0.034255099704226061787082821046821,
212 0.035262412660156681033782717998428,
213 0.036076989622888701185500318003895,
214 0.036698604498456094498018047441094,
215 0.037120549269832576114119958413599,
216 0.037334228751935040321235449094698,
217 0.037361073762679023410321241766599
309 if (
this != &integral) {
385 std::sort(bounds.begin(), bounds.end());
394 for (
int i = 0; i < bounds.size()-1; ++i) {
395 value +=
romberg(bounds[i], bounds[i+1], order);
442 bool converged =
false;
451 std::string msg =
"Requested integration order "+
453 "maximum number of iterations "+
455 "integration order or increase the (maximum) "
456 "number of iterations.";
475 #if defined(G_NAN_CHECK)
476 if (is_notanumber(s[
m_iter]) || is_infinite(s[
m_iter])) {
477 m_message =
"*** ERROR: GIntegral::romberg"
502 else if (std::abs(dss) <=
m_eps * std::abs(result)) {
508 if (std::abs(result) > 0) {
531 " exceeds absolute tolerance of "+
534 " iterations. Result "+
538 std::string origin =
"GIntegral::romberg("+
591 std::string msg =
"Function kernel not set. Please set a function "
592 "kernel before calling the method.";
613 result = 0.5*(b-a)*(y_a + y_b);
622 for (
int j = 1; j < n-1; ++j) {
635 ". Looks like n is too large.";
642 double tnm = double(it);
643 double del = (b-a)/tnm;
654 ". Step is too small to make sense.";
661 double x = a + 0.5*del;
663 for (
int j = 0; j < it; ++j, x+=del) {
675 result = 0.5*(result + (b-a)*
sum/tnm);
719 double c = 0.5*(a + b);
729 double S = (h/6.0) * (fa + 4.0*fc + fb);
732 double epsilon = (std::abs(S) > 0) ?
m_eps * S :
m_eps;
742 m_message =
"Integration uncertainty exceeds relative tolerance "
748 std::string origin =
"GIntegral::adaptive_simpson("+
773 const double alpha = std::sqrt(2.0/3.0);
774 const double beta = 1.0/std::sqrt(5.0);
775 const double x1 = 0.942882415695480;
776 const double x2 = 0.641853342345781;
777 const double x3 = 0.236383199662150;
778 const double x[12] = {0,-x1,-alpha,-x2,-beta,-x3,0.0,x3,beta,x2,alpha,x1};
779 const double eps = std::numeric_limits<double>::epsilon();
797 double m = 0.5 * (a+b);
798 double h = 0.5 * (b-a);
810 for (
int i = 1; i < 12; ++i) {
816 double i2 = (h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));
819 double i1 = (h/1470.0) * ( 77.0 * (y[0]+y[12])+
820 432.0 * (y[2]+y[10])+
825 double is = h*(0.0158271919734802 * (y[0]+y[12])+
826 0.0942738402188500 * (y[1]+y[11])+
827 0.155071987336585 * (y[2]+y[10])+
828 0.188821573960182 * (y[3]+y[9])+
829 0.199773405226859 * (y[4]+y[8])+
830 0.224926465333340 * (y[5]+y[7])+
831 0.242611071901408 * y[6]);
834 double erri1 = std::abs(i1-is);
835 double erri2 = std::abs(i2-is);
838 double r = (erri2 != 0.0) ? erri1/erri2 : 1.0;
839 double toler = (r > 0.0 && r < 1.0) ? tol/r : tol;
883 if (
m_eps < 1.12e-14) {
886 " cannot be acheived. Please relax the integration "
889 std::string origin =
"GIntegral::gauss_kronrod("+
897 double h = 0.5 * (b - a);
898 double abs_h = std::abs(h);
899 double c = 0.5 * (b + a);
908 for (
int k = 0; k < 5; ++k) {
912 double fval = fval1 + fval2;
921 for (
int k = 0; k < 5; ++k) {
925 double fval = fval1 + fval2;
934 double mean = 0.5 * res21;
936 for (
int k = 0; k < 5; ++k) {
938 (std::abs(fv1[k] - mean) + std::abs(fv2[k] - mean)) +
940 (std::abs(fv3[k] - mean) + std::abs(fv4[k] - mean)));
944 double error =
rescale_error((res21 - res10) * h, resabs, resasc);
947 if (error <
m_eps * std::abs(result)) {
950 if (std::abs(result) > 0) {
952 m_relerr = error / std::abs(result);
960 for (
int k = 0; k < 10; ++k) {
963 for (
int k = 0; k < 11; ++k) {
975 if (error <
m_eps * std::abs(result)) {
978 if (std::abs(result) > 0) {
980 m_relerr = error / std::abs(result);
988 for (
int k = 0; k < 21; ++k) {
991 for (
int k = 0; k < 22; ++k) {
1001 error =
rescale_error ((res87 - res43) * h, resabs, resasc);
1002 if (error <
m_eps * std::abs(result)) {
1005 if (std::abs(result) > 0) {
1007 m_relerr = error / std::abs(result);
1015 "absolute tolerance of "+
1021 std::string origin =
"GIntegral::gauss_kronrod("+
1049 result.append(
"=== GIntegral ===");
1067 result.append(
" (fixed: ");
1072 result.append(
" (maximum: ");
1080 result.append(
"Result accurate.");
1090 result.append(
"in standard output");
1194 std::vector<double> c(n, 0.0);
1195 std::vector<double> d(n, 0.0);
1198 double dif = std::abs(x-xa[1]);
1202 for (
int i = 0; i < n; ++i) {
1203 double dift = std::abs(x-xa[i+1]);
1217 for (
int m = 1; m < n; ++m) {
1220 for (
int i = 0; i < n-m; ++i) {
1221 double ho = xa[i+1] - x;
1222 double hp = xa[i+m+1] - x;
1223 double w = c[i+1] - d[i];
1224 double den = ho - hp;
1228 " encountered. Two values in xa array are"
1240 *dy = (2*(ns+1) < (n-m)) ? c[ns+1] : d[ns--];
1268 const double& eps,
const double& S,
1269 const double& fa,
const double& fb,
1271 const int& bottom)
const
1279 double c = 0.5*(a + b);
1281 double d = 0.5*(a + c);
1282 double e = 0.5*(c + b);
1290 double h12 = h / 12.0;
1291 double Sleft = h12 * (fa + 4.0*fd + fc);
1292 double Sright = h12 * (fc + 4.0*fe + fb);
1293 double S2 = Sleft + Sright;
1299 if (std::abs(S2 - S) <= 15.0 *
eps) {
1301 value = S2 + (S2 - S)/15.0;
1306 else if (bottom <= 0) {
1307 value = S2 + (S2 - S)/15.0;
1341 const double& toler)
const
1344 const double alpha = std::sqrt(2.0/3.0);
1345 const double beta = 1.0/std::sqrt(5.0);
1351 double m = 0.5*(a+b);
1352 double h = 0.5*(b-a);
1355 double mll = m - alpha*h;
1356 double ml = m - beta*h;
1357 double mr = m + beta*h;
1358 double mrr = m + alpha*h;
1369 double i2 = h/6.0*(fa+fb+5.0*(fml+fmr));
1372 double i1 = h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);
1375 if (std::abs(i1-i2) <= toler*is || mll <= a || b <= mrr) {
1407 const double& result_abs,
1408 const double& result_asc)
const
1411 err = std::abs(err);
1414 if (result_asc != 0.0 && err != 0.0) {
1415 double scale = std::pow((200.0 * err / result_asc), 1.5);
1417 err = result_asc * scale ;
1423 if (result_abs > 2.2250738585072014e-308 / (50.0 * 2.2204460492503131e-16)) {
1424 double min_err = 50.0 * 2.2204460492503131e-16 * result_abs;
1425 if (min_err > err) {
Exception handler interface definition.
Single parameter function abstract base class definition.
Integration class interface definition.
double sum(const GVector &vector)
Computes vector sum.
Single parameter function abstract base class.
virtual double eval(const double &x)=0
GIntegral class interface definition.
const bool & silent(void) const
Get silence flag.
const std::string & message(void) const
Return integration status message.
std::string m_message
Status message (if result is invalid)
int m_calls
Number of function calls used.
const double & eps(void) const
Get relative precision.
double adaptive_simpson_aux(const double &a, const double &b, const double &eps, const double &S, const double &fa, const double &fb, const double &fc, const int &bottom) const
Auxiliary function for adaptive Simpson's method.
const int & iter(void) const
Return number of iterations.
double adaptive_gauss_kronrod_aux(const double &a, const double &b, const double &fa, const double &fb, const double &is, const double &toler) const
Adaptive Gauss-Lobatto-Kronrod integration helper.
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
double m_relerr
Absolute integration error.
void copy_members(const GIntegral &integral)
Copy class members.
double m_eps
Requested relative integration precision.
const int & fixed_iter(void) const
Return fixed number of iterations.
const int & calls(void) const
Get number of function calls.
bool m_has_abserr
Has absolute integration error.
bool m_has_relerr
Has relative integration error.
GIntegral * clone(void) const
Clone integral.
double trapzd(const double &a, const double &b, const int &n=1, double result=0.0)
Perform Trapezoidal integration.
void clear(void)
Clear integral.
bool m_terminate
Signals termination of subdivision.
bool m_silent
Suppress integration warnings in console.
const GFunction * kernel(void) const
Get function kernel.
int m_iter
Number of iterations used.
double gauss_kronrod(const double &a, const double &b) const
Gauss-Kronrod integration.
virtual ~GIntegral(void)
Destructor.
int m_max_iter
Maximum number of iterations.
void init_members(void)
Initialise class members.
double adaptive_simpson(const double &a, const double &b) const
Adaptive Simpson's integration.
GFunction * m_kernel
Pointer to function kernel.
std::string print(const GChatter &chatter=NORMAL) const
Print integral information.
GIntegral(void)
Void constructor.
void free_members(void)
Delete class members.
GIntegral & operator=(const GIntegral &integral)
Assignment operator.
double adaptive_gauss_kronrod(const double &a, const double &b) const
Adaptive Gauss-Lobatto-Kronrod integration.
const int & max_iter(void) const
Return maximum number of iterations.
double rescale_error(double err, const double &result_abs, const double &result_asc) const
Rescale errors for Gauss-Kronrod integration.
const bool & is_valid(void) const
Signal if integration result is valid.
double m_abserr
Absolute integration error.
bool m_isvalid
Integration result valid (true=yes)
double polint(double *xa, double *ya, int n, double x, double *dy)
Perform Polynomial interpolation.
int m_fix_iter
Fixed number of iterations.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
void warning(const std::string &origin, const std::string &message)
Emits warning.