inte_kronrod_gsl.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, Jerry Gagelman and Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_INTE_GSL_KRONROD_H
24 #define O2SCL_INTE_GSL_KRONROD_H
25 
26 /** \file inte_kronrod_gsl.h
27  \brief File defining GSL-based integration coefficients, workspace,
28  and \ref o2scl::inte_kronrod_gsl
29 */
30 
31 #include <cmath>
32 #include <gsl/gsl_integration.h>
33 #include <o2scl/inte.h>
34 #include <o2scl/inte_gsl.h>
35 #include <o2scl/string_conv.h>
36 
37 /** \brief A namespace for the GSL adaptive Gauss-Kronrod integration
38  coefficients
39 
40  The storage scheme follows that of QUADPACK: symmetry about the
41  origin permits storing only those abscissae in the interval
42  \f$ [0, 1). \f$ Similiarly for the weights. The odd-indexed
43  abscissae \c xgk[1], \c xgk[3],... belong to the low-order Gauss
44  rule. For the full Gauss-Kronrod rule, the array \c wgk[] holds
45  the weights corresponding to the abscissae \c xgk[]. For the
46  low-order rule, \c wg[i] is the weight correpsonding to
47  \c xgk[2*i+1].
48 
49  <b>Documentation from GSL</b>: \n
50  Gauss quadrature weights and kronrod quadrature abscissae and
51  weights as evaluated with 80 decimal digit arithmetic by
52  L. W. Fullerton, Bell Labs, Nov. 1981.
53 */
55 
56  /** \brief Abscissae of the 15-point Kronrod rule */
57  static const double qk15_xgk[8]={
58  0.991455371120812639206854697526329,
59  0.949107912342758524526189684047851,
60  0.864864423359769072789712788640926,
61  0.741531185599394439863864773280788,
62  0.586087235467691130294144838258730,
63  0.405845151377397166906606412076961,
64  0.207784955007898467600689403773245,
65  0.000000000000000000000000000000000
66  };
67 
68  /** \brief Weights of the 7-point Gauss rule */
69  static const double qk15_wg[4] =
70  {
71  0.129484966168869693270611432679082,
72  0.279705391489276667901467771423780,
73  0.381830050505118944950369775488975,
74  0.417959183673469387755102040816327
75  };
76 
77  /** \brief Weights of the 15-point Kronrod rule */
78  static const double qk15_wgk[8] =
79  {
80  0.022935322010529224963732008058970,
81  0.063092092629978553290700663189204,
82  0.104790010322250183839876322541518,
83  0.140653259715525918745189590510238,
84  0.169004726639267902826583426598550,
85  0.190350578064785409913256402421014,
86  0.204432940075298892414161999234649,
87  0.209482141084727828012999174891714
88  };
89 
90  /** \brief Abscissae of the 21-point Kronrod rule */
91  static const double qk21_xgk[11] =
92  {
93  0.995657163025808080735527280689003,
94  0.973906528517171720077964012084452,
95  0.930157491355708226001207180059508,
96  0.865063366688984510732096688423493,
97  0.780817726586416897063717578345042,
98  0.679409568299024406234327365114874,
99  0.562757134668604683339000099272694,
100  0.433395394129247190799265943165784,
101  0.294392862701460198131126603103866,
102  0.148874338981631210884826001129720,
103  0.000000000000000000000000000000000
104  };
105 
106  /** \brief Weights of the 10-point Gauss rule */
107  static const double qk21_wg[5] =
108  {
109  0.066671344308688137593568809893332,
110  0.149451349150580593145776339657697,
111  0.219086362515982043995534934228163,
112  0.269266719309996355091226921569469,
113  0.295524224714752870173892994651338
114  };
115 
116  /** \brief Weights of the 21-point Kronrod rule */
117  static const double qk21_wgk[11] =
118  {
119  0.011694638867371874278064396062192,
120  0.032558162307964727478818972459390,
121  0.054755896574351996031381300244580,
122  0.075039674810919952767043140916190,
123  0.093125454583697605535065465083366,
124  0.109387158802297641899210590325805,
125  0.123491976262065851077958109831074,
126  0.134709217311473325928054001771707,
127  0.142775938577060080797094273138717,
128  0.147739104901338491374841515972068,
129  0.149445554002916905664936468389821
130  };
131 
132  /** \brief Abscissae of the 31-point Kronrod rule */
133  static const double qk31_xgk[16] =
134  {
135  0.998002298693397060285172840152271,
136  0.987992518020485428489565718586613,
137  0.967739075679139134257347978784337,
138  0.937273392400705904307758947710209,
139  0.897264532344081900882509656454496,
140  0.848206583410427216200648320774217,
141  0.790418501442465932967649294817947,
142  0.724417731360170047416186054613938,
143  0.650996741297416970533735895313275,
144  0.570972172608538847537226737253911,
145  0.485081863640239680693655740232351,
146  0.394151347077563369897207370981045,
147  0.299180007153168812166780024266389,
148  0.201194093997434522300628303394596,
149  0.101142066918717499027074231447392,
150  0.000000000000000000000000000000000
151  };
152 
153  /** \brief Weights of the 15-point Gauss rule */
154  static const double qk31_wg[8] =
155  {
156  0.030753241996117268354628393577204,
157  0.070366047488108124709267416450667,
158  0.107159220467171935011869546685869,
159  0.139570677926154314447804794511028,
160  0.166269205816993933553200860481209,
161  0.186161000015562211026800561866423,
162  0.198431485327111576456118326443839,
163  0.202578241925561272880620199967519
164  };
165 
166  /** \brief Weights of the 31-point Kronrod rule */
167  static const double qk31_wgk[16] =
168  {
169  0.005377479872923348987792051430128,
170  0.015007947329316122538374763075807,
171  0.025460847326715320186874001019653,
172  0.035346360791375846222037948478360,
173  0.044589751324764876608227299373280,
174  0.053481524690928087265343147239430,
175  0.062009567800670640285139230960803,
176  0.069854121318728258709520077099147,
177  0.076849680757720378894432777482659,
178  0.083080502823133021038289247286104,
179  0.088564443056211770647275443693774,
180  0.093126598170825321225486872747346,
181  0.096642726983623678505179907627589,
182  0.099173598721791959332393173484603,
183  0.100769845523875595044946662617570,
184  0.101330007014791549017374792767493
185  };
186 
187  /** \brief Abscissae of the 41-point Kronrod rule */
188  static const double qk41_xgk[21] =
189  {
190  0.998859031588277663838315576545863,
191  0.993128599185094924786122388471320,
192  0.981507877450250259193342994720217,
193  0.963971927277913791267666131197277,
194  0.940822633831754753519982722212443,
195  0.912234428251325905867752441203298,
196  0.878276811252281976077442995113078,
197  0.839116971822218823394529061701521,
198  0.795041428837551198350638833272788,
199  0.746331906460150792614305070355642,
200  0.693237656334751384805490711845932,
201  0.636053680726515025452836696226286,
202  0.575140446819710315342946036586425,
203  0.510867001950827098004364050955251,
204  0.443593175238725103199992213492640,
205  0.373706088715419560672548177024927,
206  0.301627868114913004320555356858592,
207  0.227785851141645078080496195368575,
208  0.152605465240922675505220241022678,
209  0.076526521133497333754640409398838,
210  0.000000000000000000000000000000000
211  };
212 
213  /** \brief Weights of the 20-point Gauss rule */
214  static const double qk41_wg[11] =
215  {
216  0.017614007139152118311861962351853,
217  0.040601429800386941331039952274932,
218  0.062672048334109063569506535187042,
219  0.083276741576704748724758143222046,
220  0.101930119817240435036750135480350,
221  0.118194531961518417312377377711382,
222  0.131688638449176626898494499748163,
223  0.142096109318382051329298325067165,
224  0.149172986472603746787828737001969,
225  0.152753387130725850698084331955098
226  };
227 
228  /** \brief Weights of the 41-point Kronrod rule */
229  static const double qk41_wgk[21] =
230  {
231  0.003073583718520531501218293246031,
232  0.008600269855642942198661787950102,
233  0.014626169256971252983787960308868,
234  0.020388373461266523598010231432755,
235  0.025882133604951158834505067096153,
236  0.031287306777032798958543119323801,
237  0.036600169758200798030557240707211,
238  0.041668873327973686263788305936895,
239  0.046434821867497674720231880926108,
240  0.050944573923728691932707670050345,
241  0.055195105348285994744832372419777,
242  0.059111400880639572374967220648594,
243  0.062653237554781168025870122174255,
244  0.065834597133618422111563556969398,
245  0.068648672928521619345623411885368,
246  0.071054423553444068305790361723210,
247  0.073030690332786667495189417658913,
248  0.074582875400499188986581418362488,
249  0.075704497684556674659542775376617,
250  0.076377867672080736705502835038061,
251  0.076600711917999656445049901530102
252  };
253 
254  /** \brief Abscissae of the 51-point Kronrod rule */
255  static const double qk51_xgk[26] =
256  {
257  0.999262104992609834193457486540341,
258  0.995556969790498097908784946893902,
259  0.988035794534077247637331014577406,
260  0.976663921459517511498315386479594,
261  0.961614986425842512418130033660167,
262  0.942974571228974339414011169658471,
263  0.920747115281701561746346084546331,
264  0.894991997878275368851042006782805,
265  0.865847065293275595448996969588340,
266  0.833442628760834001421021108693570,
267  0.797873797998500059410410904994307,
268  0.759259263037357630577282865204361,
269  0.717766406813084388186654079773298,
270  0.673566368473468364485120633247622,
271  0.626810099010317412788122681624518,
272  0.577662930241222967723689841612654,
273  0.526325284334719182599623778158010,
274  0.473002731445714960522182115009192,
275  0.417885382193037748851814394594572,
276  0.361172305809387837735821730127641,
277  0.303089538931107830167478909980339,
278  0.243866883720988432045190362797452,
279  0.183718939421048892015969888759528,
280  0.122864692610710396387359818808037,
281  0.061544483005685078886546392366797,
282  0.000000000000000000000000000000000
283  };
284 
285  /** \brief Weights of the 25-point Gauss rule */
286  static const double qk51_wg[13] =
287  {
288  0.011393798501026287947902964113235,
289  0.026354986615032137261901815295299,
290  0.040939156701306312655623487711646,
291  0.054904695975835191925936891540473,
292  0.068038333812356917207187185656708,
293  0.080140700335001018013234959669111,
294  0.091028261982963649811497220702892,
295  0.100535949067050644202206890392686,
296  0.108519624474263653116093957050117,
297  0.114858259145711648339325545869556,
298  0.119455763535784772228178126512901,
299  0.122242442990310041688959518945852,
300  0.123176053726715451203902873079050
301  };
302 
303  /** \brief Weights of the 51-point Kronrod rule */
304  static const double qk51_wgk[26] =
305  {
306  0.001987383892330315926507851882843,
307  0.005561932135356713758040236901066,
308  0.009473973386174151607207710523655,
309  0.013236229195571674813656405846976,
310  0.016847817709128298231516667536336,
311  0.020435371145882835456568292235939,
312  0.024009945606953216220092489164881,
313  0.027475317587851737802948455517811,
314  0.030792300167387488891109020215229,
315  0.034002130274329337836748795229551,
316  0.037116271483415543560330625367620,
317  0.040083825504032382074839284467076,
318  0.042872845020170049476895792439495,
319  0.045502913049921788909870584752660,
320  0.047982537138836713906392255756915,
321  0.050277679080715671963325259433440,
322  0.052362885806407475864366712137873,
323  0.054251129888545490144543370459876,
324  0.055950811220412317308240686382747,
325  0.057437116361567832853582693939506,
326  0.058689680022394207961974175856788,
327  0.059720340324174059979099291932562,
328  0.060539455376045862945360267517565,
329  0.061128509717053048305859030416293,
330  0.061471189871425316661544131965264,
331  0.061580818067832935078759824240066
332  };
333 
334  /** \brief Abscissae of the 61-point Kronrod rule */
335  static const double qk61_xgk[31] =
336  {
337  0.999484410050490637571325895705811,
338  0.996893484074649540271630050918695,
339  0.991630996870404594858628366109486,
340  0.983668123279747209970032581605663,
341  0.973116322501126268374693868423707,
342  0.960021864968307512216871025581798,
343  0.944374444748559979415831324037439,
344  0.926200047429274325879324277080474,
345  0.905573307699907798546522558925958,
346  0.882560535792052681543116462530226,
347  0.857205233546061098958658510658944,
348  0.829565762382768397442898119732502,
349  0.799727835821839083013668942322683,
350  0.767777432104826194917977340974503,
351  0.733790062453226804726171131369528,
352  0.697850494793315796932292388026640,
353  0.660061064126626961370053668149271,
354  0.620526182989242861140477556431189,
355  0.579345235826361691756024932172540,
356  0.536624148142019899264169793311073,
357  0.492480467861778574993693061207709,
358  0.447033769538089176780609900322854,
359  0.400401254830394392535476211542661,
360  0.352704725530878113471037207089374,
361  0.304073202273625077372677107199257,
362  0.254636926167889846439805129817805,
363  0.204525116682309891438957671002025,
364  0.153869913608583546963794672743256,
365  0.102806937966737030147096751318001,
366  0.051471842555317695833025213166723,
367  0.000000000000000000000000000000000
368  };
369 
370  /** \brief Weights of the 30-point Gauss rule */
371  static const double qk61_wg[15] =
372  {
373  0.007968192496166605615465883474674,
374  0.018466468311090959142302131912047,
375  0.028784707883323369349719179611292,
376  0.038799192569627049596801936446348,
377  0.048402672830594052902938140422808,
378  0.057493156217619066481721689402056,
379  0.065974229882180495128128515115962,
380  0.073755974737705206268243850022191,
381  0.080755895229420215354694938460530,
382  0.086899787201082979802387530715126,
383  0.092122522237786128717632707087619,
384  0.096368737174644259639468626351810,
385  0.099593420586795267062780282103569,
386  0.101762389748405504596428952168554,
387  0.102852652893558840341285636705415
388  };
389 
390  /** \brief Weights of the 61-point Kronrod rule */
391  static const double qk61_wgk[31] =
392  {
393  0.001389013698677007624551591226760,
394  0.003890461127099884051267201844516,
395  0.006630703915931292173319826369750,
396  0.009273279659517763428441146892024,
397  0.011823015253496341742232898853251,
398  0.014369729507045804812451432443580,
399  0.016920889189053272627572289420322,
400  0.019414141193942381173408951050128,
401  0.021828035821609192297167485738339,
402  0.024191162078080601365686370725232,
403  0.026509954882333101610601709335075,
404  0.028754048765041292843978785354334,
405  0.030907257562387762472884252943092,
406  0.032981447057483726031814191016854,
407  0.034979338028060024137499670731468,
408  0.036882364651821229223911065617136,
409  0.038678945624727592950348651532281,
410  0.040374538951535959111995279752468,
411  0.041969810215164246147147541285970,
412  0.043452539701356069316831728117073,
413  0.044814800133162663192355551616723,
414  0.046059238271006988116271735559374,
415  0.047185546569299153945261478181099,
416  0.048185861757087129140779492298305,
417  0.049055434555029778887528165367238,
418  0.049795683427074206357811569379942,
419  0.050405921402782346840893085653585,
420  0.050881795898749606492297473049805,
421  0.051221547849258772170656282604944,
422  0.051426128537459025933862879215781,
423  0.051494729429451567558340433647099
424  };
425 
426 }
427 
428 #ifndef DOXYGEN_NO_O2NS
429 namespace o2scl {
430 #endif
431 
432  /** \brief Integration workspace for the GSL integrators
433 
434  This is a simple rewrite of \c inte_gslgration_workspace
435  into a class.
436 
437  QUADPACK workspace documentation:
438  \verbatim
439  c parameters (meaning at output)
440  c limit - integer
441  c maximum number of error estimates the list
442  c can contain
443  c last - integer
444  c number of error estimates currently in the list
445  c maxerr - integer
446  c maxerr points to the nrmax-th largest error
447  c estimate currently in the list
448  c ermax - double precision
449  c nrmax-th largest error estimate
450  c ermax = elist(maxerr)
451  c elist - double precision
452  c vector of dimension last containing
453  c the error estimates
454  c iord - integer
455  c vector of dimension last, the first k elements
456  c of which contain pointers to the error
457  c estimates, such that
458  c elist(iord(1)),..., elist(iord(k))
459  c form a decreasing sequence, with
460  c k = last if last.le.(limit/2+2), and
461  c k = limit+1-last otherwise
462  c nrmax - integer
463  c maxerr = iord(nrmax)
464  \endverbatim
465  \verbatim
466  c alist - real
467  c vector of dimension at least limit, the first
468  c last elements of which are the left
469  c end points of the subintervals in the partition
470  c of the given integration range (a,b)
471  c blist - real
472  c vector of dimension at least limit, the first
473  c last elements of which are the right
474  c end points of the subintervals in the partition
475  c of the given integration range (a,b)
476  c rlist - real
477  c vector of dimension at least limit, the first
478  c last elements of which are the
479  c integral approximations on the subintervals
480  c elist - real
481  c vector of dimension at least limit, the first
482  c last elements of which are the moduli of the
483  c absolute error estimates on the subintervals
484  c iord - integer
485  c vector of dimension at least limit, the first k
486  c elements of which are pointers to the
487  c error estimates over the subintervals,
488  c such that elist(iord(1)), ...,
489  c elist(iord(k)) form a decreasing sequence,
490  c with k = last if last.le.(limit/2+2), and
491  c k = limit+1-last otherwise
492  c last - integer
493  c number of subintervals actually produced in the
494  c subdivision process
495  \endverbatim
496 
497  */
499 
500  public:
501 
503 
505 
506  /// Maximum number of subintervals allocated
507  size_t limit;
508  /// Current number of subintervals being used
509  size_t size;
510  /// Counter for extrapolation routine
511  size_t nrmax;
512  /// Index of current subinterval
513  size_t i;
514  /// Depth of subdivisions reached
516  /// Left endpoints of subintervals
517  double *alist;
518  /// Right endpoints of subintervals
519  double *blist;
520  /// Integral approximations on subintervals
521  double *rlist;
522  /// Integral error estimates
523  double *elist;
524  /// Linear ordering vector for sort routine
525  size_t *order;
526  /// Numbers of subdivisions made
527  size_t *level;
528 
529  /// Allocate a workspace
530  int allocate(size_t sz);
531 
532  /// Free allocated workspace memory
533  int free();
534 
535  /** \brief Initialize the workspace for an integration with limits
536  \c a and \c b.
537  */
538  int initialise(double a, double b);
539 
540  /** \brief Update the workspace with the result and error
541  from the first integration
542  */
543  int set_initial_result(double result, double error);
544 
545  /** \brief Retrieve the ith result from the
546  workspace stack
547 
548  The workspace variable \c i is used to specify which
549  interval is requested.
550  */
551  int retrieve(double *a, double *b, double *r, double *e) const;
552 
553  /** \brief Sort the workspace stack
554 
555  This routine maintains the descending ordering in the list of
556  the local error estimated resulting from the interval
557  subdivision process. at each call two error estimates are
558  inserted using the sequential search method, top-down for the
559  largest error estimate and bottom-up for the smallest error
560  estimate.
561 
562  Originally written in QUADPACK by R. Piessens and E. de
563  Doncker, translated into C for GSL by Brian Gough, and then
564  rewritten for \o2.
565  */
566  int qpsrt();
567 
568  /** \brief Determine which new subinterval to add to the workspace
569  stack and perform update
570  */
571  int update(double a1, double b1, double area1, double error1,
572  double a2, double b2, double area2, double error2);
573 
574  /// Add up all of the contributions to construct the final result
575  double sum_results();
576 
577  /** \brief Test whether the proposed subdivision falls
578  before floating-point precision
579  */
580  int subinterval_too_small(double a1, double a2, double b2);
581 
582  /// Push a new interval to the workspace stack
583  void append_interval(double a1, double b1, double area1, double error1);
584 
585  };
586 
587  /** \brief Basic Gauss-Kronrod integration class (GSL)
588 
589  This class provides the basic Gauss-Kronrod integration function
590  and integration workspace for some of the GSL-based integration
591  classes.
592 
593  The main function of interest is \ref set_rule(), which
594  sets the integration rule for the GSL integration classes
595  which inherit from this base class. The argument to
596  set rule should be selected from the following list:
597  - 1: 15-point Gauss-Kronrod integration rule (default)
598  - 2: 21-point rule
599  - 3: 31-point rule
600  - 4: 41-point rule
601  - 5: 51-point rule
602  - 6: 61-point rule
603  Any value other than 1-6 forces the error handler to be
604  called. All classes default to the 15-point rule,
605  except for \ref o2scl::inte_qags_gsl which defaults to
606  the 21-point rule for singularities.
607 
608  The integration coefficients for use with this class and
609  associated children are stored in the \ref o2scl_inte_gk_coeffs
610  namespace.
611 
612  \future Convert 'double *' objects to ubvector
613  */
614  template<class func_t> class inte_kronrod_gsl : public inte_gsl,
615  public inte<func_t> {
616 
617  protected:
618 
619  /// The integration workspace
621 
622  /// Size of Gauss-Kronrod arrays
623  int n_gk;
624 
625  /// Gauss-Kronrod abscissae pointer
626  const double *x_gk;
627 
628  /// Gauss weight pointer
629  const double *w_g;
630 
631  /// Gauss-Kronrod weight pointer
632  const double *w_gk;
633 
634  /// Scratch space
635  double *f_v1;
636 
637  /// Scratch space
638  double *f_v2;
639 
640  public:
641 
642  inte_kronrod_gsl() {
643  w=new inte_workspace_gsl;
644  w->allocate(1000);
645  n_gk=0;
646  set_rule(1);
647  }
648 
649  ~inte_kronrod_gsl() {
650  w->free();
651  delete w;
652  if (n_gk>0) {
653  delete[] f_v1;
654  delete[] f_v2;
655  }
656  }
657 
658  /** \brief Get the Gauss-Kronrod integration rule
659 
660  This returns the index of the GSL integration rule
661  a number between 1 and 6 (inclusive)
662  */
663  int get_rule() {
664  switch (n_gk) {
665  case 8:
666  return 1;
667  case 11:
668  return 2;
669  case 16:
670  return 3;
671  case 21:
672  return 4;
673  case 26:
674  return 5;
675  case 31:
676  return 6;
677  }
678  return 0;
679  }
680 
681  /// Set the Gauss-Kronrod integration rule to be used
682  void set_rule(int rule) {
683 
684  using namespace o2scl_inte_gk_coeffs;
685 
686  if (n_gk > 0) {
687  delete[] f_v1;
688  delete[] f_v2;
689  }
690 
691  switch (rule) {
692  case 1:
693  n_gk=8;
694  x_gk=qk15_xgk;
695  w_g=qk15_wg;
696  w_gk=qk15_wgk;
697  break;
698  case 2:
699  n_gk=11;
700  x_gk=qk21_xgk;
701  w_g=qk21_wg;
702  w_gk=qk21_wgk;
703  break;
704  case 3:
705  n_gk=16;
706  x_gk=qk31_xgk;
707  w_g=qk31_wg;
708  w_gk=qk31_wgk;
709  break;
710  case 4:
711  n_gk=21;
712  x_gk=qk41_xgk;
713  w_g=qk41_wg;
714  w_gk=qk41_wgk;
715  break;
716  case 5:
717  n_gk=26;
718  x_gk=qk51_xgk;
719  w_g=qk51_wg;
720  w_gk=qk51_wgk;
721  break;
722  case 6:
723  n_gk=31;
724  x_gk=qk61_xgk;
725  w_g=qk61_wg;
726  w_gk=qk61_wgk;
727  break;
728  default:
729  O2SCL_ERR("Invalid rule.",exc_einval);
730  break;
731  }
732 
733  f_v1=new double[n_gk];
734  f_v2=new double[n_gk];
735 
736  return;
737  }
738 
739  /** \brief Set the limit for the number of subdivisions of
740  the integration region (default 1000)
741 
742  If the value of \c size is zero, the error handler will
743  be called.
744  */
745  int set_limit(size_t lim) {
746  if (lim==0) {
747  O2SCL_ERR2("Requested zero limit in ",
748  "inte_kronrod_gsl::set_limit().",exc_einval);
749  }
750  w->free();
751  w->allocate(lim);
752  return 0;
753  }
754 
755  /** \brief The base Gauss-Kronrod integration function template
756 
757  Given abcissas and weights, this performs the integration of
758  \c func between \c a and \c b, providing a result with
759  uncertainties.
760 
761  The Gauss-Kronrod rule uses \f$ 2m+1 \f$ abscissae \f$ x_1,
762  \ldots, x_{2m+1} \in (a, b) \f$ to estimate the integral of a
763  function as a linear combination of values,
764  \f[
765  \int_a^b f~dx \; \approx \;
766  Q_m^{GK}f \; \equiv \;
767  \sum_{k=1}^{2m+1} \beta_k f(x_k),
768  \f]
769  where the weights \f$ \beta_1,\ldots, \beta_{2m+1} \f$ are
770  intrinsic to the abscissae. The data are designed so that the
771  even-indexed abscissae yield a coarser estimate,
772  \f[
773  Q_m^{G}f \; \equiv \;
774  \sum_{j=1}^{m} \alpha_j f(x_{2j}),
775  \f]
776  and their difference \f$ |Q_m^Gf - Q_m^{GK}f| \f$ is the "raw"
777  error estimate. The various quantities that the function
778  computes are
779  \f{eqnarray*}{
780  \mathtt{result} &=& Q_m^{GK}f, \\
781  \mathtt{resabs} &=& Q_m^{GK}|f|, \\
782  \mathtt{resasc} &=& Q_m^{GK}(|f| - \mu),
783  \quad \mu \equiv \frac{Q_m^{GK}f}{b - a}.
784  \f}
785  The "absolute" error \c abserr is computed from the raw error
786  value using the function \ref inte_gsl::rescale_error.
787 
788  This function is designed for use with the values given in the
789  o2scl_inte_gk_coeffs namespace.
790 
791  This function never calls the error handler.
792 
793  \future This function, in principle, is an integration
794  object in and of itself, and could be implemented separately
795  as an object of type \ref o2scl::inte.
796  */
797  template<class func2_t> void gauss_kronrod_base
798  (func2_t &func, double a, double b, double *result,
799  double *abserr, double *resabs, double *resasc) {
800 
801  const double center=0.5*(a+b);
802  const double half_length=0.5*(b-a);
803  const double abs_half_length=fabs (half_length);
804 
805  double f_center;
806  f_center=func(center);
807 
808  double result_gauss=0.0;
809  double result_kronrod=f_center*this->w_gk[this->n_gk-1];
810 
811  double result_abs=fabs (result_kronrod);
812  double result_asc=0.0;
813  double mean=0.0, err=0.0;
814 
815  if (this->n_gk % 2 == 0) {
816  result_gauss=f_center*this->w_g[this->n_gk / 2-1];
817  }
818 
819  // Sum up results from left half of interval
820  for (int j=0; j < (this->n_gk-1) / 2; j++) {
821 
822  const int jtw=j*2+1;
823  const double abscissa=half_length*this->x_gk[jtw];
824 
825  double fval1, fval2;
826  fval1=func(center-abscissa);
827  fval2=func(center+abscissa);
828 
829  const double fsum=fval1+fval2;
830  this->f_v1[jtw]=fval1;
831  this->f_v2[jtw]=fval2;
832  result_gauss+=this->w_g[j]*fsum;
833  result_kronrod+=this->w_gk[jtw]*fsum;
834  result_abs+=this->w_gk[jtw]*(fabs(fval1)+fabs(fval2));
835  }
836 
837  // Sum up results from right half of interval
838  for (int j=0; j < this->n_gk / 2; j++) {
839 
840  int jtwm1=j*2;
841  const double abscissa=half_length*this->x_gk[jtwm1];
842 
843  double fval1, fval2;
844  fval1=func(center-abscissa);
845  fval2=func(center+abscissa);
846 
847  this->f_v1[jtwm1]=fval1;
848  this->f_v2[jtwm1]=fval2;
849  result_kronrod+=this->w_gk[jtwm1]*(fval1+fval2);
850  result_abs+=this->w_gk[jtwm1]*(fabs(fval1)+fabs(fval2));
851  }
852 
853  mean=result_kronrod*0.5;
854 
855  result_asc=this->w_gk[this->n_gk-1]*fabs(f_center-mean);
856 
857  for (int j=0;j<this->n_gk-1;j++) {
858  result_asc+=this->w_gk[j]*(fabs(this->f_v1[j]-mean)+
859  fabs(this->f_v2[j]-mean));
860  }
861 
862  /* Scale by the width of the integration region */
863 
864  err=(result_kronrod-result_gauss)*half_length;
865 
866  result_kronrod*=half_length;
867  result_abs*=abs_half_length;
868  result_asc*=abs_half_length;
869 
870  *result=result_kronrod;
871  *resabs=result_abs;
872  *resasc=result_asc;
873  *abserr=rescale_error(err,result_abs,result_asc);
874 
875  return;
876  }
877 
878  /** \brief Integration wrapper for user-specified function type
879  */
880  virtual void gauss_kronrod
881  (func_t &func, double a, double b,
882  double *result, double *abserr, double *resabs, double *resasc) {
883  return gauss_kronrod_base<func_t>
884  (func,a,b,result,abserr,resabs,resasc);
885 
886  }
887 
888  };
889 
890 #ifndef DOXYGEN_NO_O2NS
891 }
892 #endif
893 
894 #endif
o2scl::inte_kronrod_gsl::gauss_kronrod_base
void gauss_kronrod_base(func2_t &func, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
The base Gauss-Kronrod integration function template.
Definition: inte_kronrod_gsl.h:798
o2scl::inte_kronrod_gsl::n_gk
int n_gk
Size of Gauss-Kronrod arrays.
Definition: inte_kronrod_gsl.h:623
o2scl_inte_gk_coeffs::qk61_wgk
static const double qk61_wgk[31]
Weights of the 61-point Kronrod rule.
Definition: inte_kronrod_gsl.h:391
o2scl_inte_gk_coeffs::qk41_wg
static const double qk41_wg[11]
Weights of the 20-point Gauss rule.
Definition: inte_kronrod_gsl.h:214
o2scl::inte_workspace_gsl::maximum_level
size_t maximum_level
Depth of subdivisions reached.
Definition: inte_kronrod_gsl.h:515
o2scl::inte_kronrod_gsl::x_gk
const double * x_gk
Gauss-Kronrod abscissae pointer.
Definition: inte_kronrod_gsl.h:626
o2scl::inte_kronrod_gsl::w
inte_workspace_gsl * w
The integration workspace.
Definition: inte_kronrod_gsl.h:620
o2scl::inte_workspace_gsl::retrieve
int retrieve(double *a, double *b, double *r, double *e) const
Retrieve the ith result from the workspace stack.
o2scl_inte_gk_coeffs::qk31_wg
static const double qk31_wg[8]
Weights of the 15-point Gauss rule.
Definition: inte_kronrod_gsl.h:154
o2scl::inte_workspace_gsl::set_initial_result
int set_initial_result(double result, double error)
Update the workspace with the result and error from the first integration.
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::inte_kronrod_gsl::set_rule
void set_rule(int rule)
Set the Gauss-Kronrod integration rule to be used.
Definition: inte_kronrod_gsl.h:682
o2scl::inte_workspace_gsl::elist
double * elist
Integral error estimates.
Definition: inte_kronrod_gsl.h:523
o2scl::inte_workspace_gsl::nrmax
size_t nrmax
Counter for extrapolation routine.
Definition: inte_kronrod_gsl.h:511
o2scl_inte_gk_coeffs::qk31_wgk
static const double qk31_wgk[16]
Weights of the 31-point Kronrod rule.
Definition: inte_kronrod_gsl.h:167
o2scl_inte_gk_coeffs::qk15_xgk
static const double qk15_xgk[8]
Abscissae of the 15-point Kronrod rule.
Definition: inte_kronrod_gsl.h:57
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl_inte_gk_coeffs::qk61_wg
static const double qk61_wg[15]
Weights of the 30-point Gauss rule.
Definition: inte_kronrod_gsl.h:371
o2scl::inte_workspace_gsl::initialise
int initialise(double a, double b)
Initialize the workspace for an integration with limits a and b.
o2scl::inte_workspace_gsl::order
size_t * order
Linear ordering vector for sort routine.
Definition: inte_kronrod_gsl.h:525
o2scl::inte_workspace_gsl::i
size_t i
Index of current subinterval.
Definition: inte_kronrod_gsl.h:513
o2scl_inte_gk_coeffs::qk21_wg
static const double qk21_wg[5]
Weights of the 10-point Gauss rule.
Definition: inte_kronrod_gsl.h:107
o2scl::inte_workspace_gsl::level
size_t * level
Numbers of subdivisions made.
Definition: inte_kronrod_gsl.h:527
o2scl::inte_workspace_gsl::allocate
int allocate(size_t sz)
Allocate a workspace.
o2scl::inte_workspace_gsl::qpsrt
int qpsrt()
Sort the workspace stack.
o2scl::inte_workspace_gsl::alist
double * alist
Left endpoints of subintervals.
Definition: inte_kronrod_gsl.h:517
o2scl::inte_workspace_gsl::sum_results
double sum_results()
Add up all of the contributions to construct the final result.
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::inte_gsl
GSL integration base.
Definition: inte_gsl.h:62
o2scl_inte_gk_coeffs::qk51_xgk
static const double qk51_xgk[26]
Abscissae of the 51-point Kronrod rule.
Definition: inte_kronrod_gsl.h:255
o2scl::inte_workspace_gsl::update
int update(double a1, double b1, double area1, double error1, double a2, double b2, double area2, double error2)
Determine which new subinterval to add to the workspace stack and perform update.
o2scl::inte_kronrod_gsl::f_v2
double * f_v2
Scratch space.
Definition: inte_kronrod_gsl.h:638
o2scl::inte_kronrod_gsl::gauss_kronrod
virtual void gauss_kronrod(func_t &func, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
Integration wrapper for user-specified function type.
Definition: inte_kronrod_gsl.h:881
o2scl_inte_gk_coeffs::qk41_xgk
static const double qk41_xgk[21]
Abscissae of the 41-point Kronrod rule.
Definition: inte_kronrod_gsl.h:188
o2scl::inte_kronrod_gsl::set_limit
int set_limit(size_t lim)
Set the limit for the number of subdivisions of the integration region (default 1000)
Definition: inte_kronrod_gsl.h:745
o2scl_inte_gk_coeffs::qk21_xgk
static const double qk21_xgk[11]
Abscissae of the 21-point Kronrod rule.
Definition: inte_kronrod_gsl.h:91
o2scl_inte_gk_coeffs::qk41_wgk
static const double qk41_wgk[21]
Weights of the 41-point Kronrod rule.
Definition: inte_kronrod_gsl.h:229
o2scl_inte_gk_coeffs::qk15_wg
static const double qk15_wg[4]
Weights of the 7-point Gauss rule.
Definition: inte_kronrod_gsl.h:69
o2scl::inte_workspace_gsl::append_interval
void append_interval(double a1, double b1, double area1, double error1)
Push a new interval to the workspace stack.
o2scl_inte_gk_coeffs::qk31_xgk
static const double qk31_xgk[16]
Abscissae of the 31-point Kronrod rule.
Definition: inte_kronrod_gsl.h:133
o2scl::inte
Base integration class [abstract base].
Definition: inte.h:51
o2scl_inte_gk_coeffs
A namespace for the GSL adaptive Gauss-Kronrod integration coefficients.
Definition: inte_kronrod_gsl.h:54
o2scl::inte_workspace_gsl::blist
double * blist
Right endpoints of subintervals.
Definition: inte_kronrod_gsl.h:519
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::inte_workspace_gsl::limit
size_t limit
Maximum number of subintervals allocated.
Definition: inte_kronrod_gsl.h:507
o2scl::inte_workspace_gsl::free
int free()
Free allocated workspace memory.
o2scl_inte_gk_coeffs::qk51_wgk
static const double qk51_wgk[26]
Weights of the 51-point Kronrod rule.
Definition: inte_kronrod_gsl.h:304
o2scl::inte_kronrod_gsl
Basic Gauss-Kronrod integration class (GSL)
Definition: inte_kronrod_gsl.h:614
o2scl::inte_kronrod_gsl::w_gk
const double * w_gk
Gauss-Kronrod weight pointer.
Definition: inte_kronrod_gsl.h:632
o2scl::inte_gsl::rescale_error
double rescale_error(double err, const double result_abs, const double result_asc)
QUADPACK's nonlinear rescaling of the absolute-error estimate.
Definition: inte_gsl.h:106
o2scl::inte_workspace_gsl
Integration workspace for the GSL integrators.
Definition: inte_kronrod_gsl.h:498
o2scl::inte_kronrod_gsl::f_v1
double * f_v1
Scratch space.
Definition: inte_kronrod_gsl.h:635
o2scl::inte_workspace_gsl::rlist
double * rlist
Integral approximations on subintervals.
Definition: inte_kronrod_gsl.h:521
o2scl_inte_gk_coeffs::qk21_wgk
static const double qk21_wgk[11]
Weights of the 21-point Kronrod rule.
Definition: inte_kronrod_gsl.h:117
o2scl_inte_gk_coeffs::qk51_wg
static const double qk51_wg[13]
Weights of the 25-point Gauss rule.
Definition: inte_kronrod_gsl.h:286
o2scl_inte_gk_coeffs::qk15_wgk
static const double qk15_wgk[8]
Weights of the 15-point Kronrod rule.
Definition: inte_kronrod_gsl.h:78
o2scl::inte_workspace_gsl::size
size_t size
Current number of subintervals being used.
Definition: inte_kronrod_gsl.h:509
o2scl::inte_kronrod_gsl::w_g
const double * w_g
Gauss weight pointer.
Definition: inte_kronrod_gsl.h:629
o2scl::inte_workspace_gsl::subinterval_too_small
int subinterval_too_small(double a1, double a2, double b2)
Test whether the proposed subdivision falls before floating-point precision.
o2scl::inte_kronrod_gsl::get_rule
int get_rule()
Get the Gauss-Kronrod integration rule.
Definition: inte_kronrod_gsl.h:663
o2scl_inte_gk_coeffs::qk61_xgk
static const double qk61_xgk[31]
Abscissae of the 61-point Kronrod rule.
Definition: inte_kronrod_gsl.h:335

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).