24 #ifndef O2SCL_INTE_GSL_KRONROD_H 25 #define O2SCL_INTE_GSL_KRONROD_H 33 #include <gsl/gsl_integration.h> 34 #include <o2scl/inte.h> 35 #include <o2scl/inte_gsl.h> 36 #include <o2scl/string_conv.h> 59 0.991455371120812639206854697526329,
60 0.949107912342758524526189684047851,
61 0.864864423359769072789712788640926,
62 0.741531185599394439863864773280788,
63 0.586087235467691130294144838258730,
64 0.405845151377397166906606412076961,
65 0.207784955007898467600689403773245,
66 0.000000000000000000000000000000000
72 0.129484966168869693270611432679082,
73 0.279705391489276667901467771423780,
74 0.381830050505118944950369775488975,
75 0.417959183673469387755102040816327
81 0.022935322010529224963732008058970,
82 0.063092092629978553290700663189204,
83 0.104790010322250183839876322541518,
84 0.140653259715525918745189590510238,
85 0.169004726639267902826583426598550,
86 0.190350578064785409913256402421014,
87 0.204432940075298892414161999234649,
88 0.209482141084727828012999174891714
94 0.995657163025808080735527280689003,
95 0.973906528517171720077964012084452,
96 0.930157491355708226001207180059508,
97 0.865063366688984510732096688423493,
98 0.780817726586416897063717578345042,
99 0.679409568299024406234327365114874,
100 0.562757134668604683339000099272694,
101 0.433395394129247190799265943165784,
102 0.294392862701460198131126603103866,
103 0.148874338981631210884826001129720,
104 0.000000000000000000000000000000000
110 0.066671344308688137593568809893332,
111 0.149451349150580593145776339657697,
112 0.219086362515982043995534934228163,
113 0.269266719309996355091226921569469,
114 0.295524224714752870173892994651338
120 0.011694638867371874278064396062192,
121 0.032558162307964727478818972459390,
122 0.054755896574351996031381300244580,
123 0.075039674810919952767043140916190,
124 0.093125454583697605535065465083366,
125 0.109387158802297641899210590325805,
126 0.123491976262065851077958109831074,
127 0.134709217311473325928054001771707,
128 0.142775938577060080797094273138717,
129 0.147739104901338491374841515972068,
130 0.149445554002916905664936468389821
136 0.998002298693397060285172840152271,
137 0.987992518020485428489565718586613,
138 0.967739075679139134257347978784337,
139 0.937273392400705904307758947710209,
140 0.897264532344081900882509656454496,
141 0.848206583410427216200648320774217,
142 0.790418501442465932967649294817947,
143 0.724417731360170047416186054613938,
144 0.650996741297416970533735895313275,
145 0.570972172608538847537226737253911,
146 0.485081863640239680693655740232351,
147 0.394151347077563369897207370981045,
148 0.299180007153168812166780024266389,
149 0.201194093997434522300628303394596,
150 0.101142066918717499027074231447392,
151 0.000000000000000000000000000000000
157 0.030753241996117268354628393577204,
158 0.070366047488108124709267416450667,
159 0.107159220467171935011869546685869,
160 0.139570677926154314447804794511028,
161 0.166269205816993933553200860481209,
162 0.186161000015562211026800561866423,
163 0.198431485327111576456118326443839,
164 0.202578241925561272880620199967519
170 0.005377479872923348987792051430128,
171 0.015007947329316122538374763075807,
172 0.025460847326715320186874001019653,
173 0.035346360791375846222037948478360,
174 0.044589751324764876608227299373280,
175 0.053481524690928087265343147239430,
176 0.062009567800670640285139230960803,
177 0.069854121318728258709520077099147,
178 0.076849680757720378894432777482659,
179 0.083080502823133021038289247286104,
180 0.088564443056211770647275443693774,
181 0.093126598170825321225486872747346,
182 0.096642726983623678505179907627589,
183 0.099173598721791959332393173484603,
184 0.100769845523875595044946662617570,
185 0.101330007014791549017374792767493
191 0.998859031588277663838315576545863,
192 0.993128599185094924786122388471320,
193 0.981507877450250259193342994720217,
194 0.963971927277913791267666131197277,
195 0.940822633831754753519982722212443,
196 0.912234428251325905867752441203298,
197 0.878276811252281976077442995113078,
198 0.839116971822218823394529061701521,
199 0.795041428837551198350638833272788,
200 0.746331906460150792614305070355642,
201 0.693237656334751384805490711845932,
202 0.636053680726515025452836696226286,
203 0.575140446819710315342946036586425,
204 0.510867001950827098004364050955251,
205 0.443593175238725103199992213492640,
206 0.373706088715419560672548177024927,
207 0.301627868114913004320555356858592,
208 0.227785851141645078080496195368575,
209 0.152605465240922675505220241022678,
210 0.076526521133497333754640409398838,
211 0.000000000000000000000000000000000
217 0.017614007139152118311861962351853,
218 0.040601429800386941331039952274932,
219 0.062672048334109063569506535187042,
220 0.083276741576704748724758143222046,
221 0.101930119817240435036750135480350,
222 0.118194531961518417312377377711382,
223 0.131688638449176626898494499748163,
224 0.142096109318382051329298325067165,
225 0.149172986472603746787828737001969,
226 0.152753387130725850698084331955098
232 0.003073583718520531501218293246031,
233 0.008600269855642942198661787950102,
234 0.014626169256971252983787960308868,
235 0.020388373461266523598010231432755,
236 0.025882133604951158834505067096153,
237 0.031287306777032798958543119323801,
238 0.036600169758200798030557240707211,
239 0.041668873327973686263788305936895,
240 0.046434821867497674720231880926108,
241 0.050944573923728691932707670050345,
242 0.055195105348285994744832372419777,
243 0.059111400880639572374967220648594,
244 0.062653237554781168025870122174255,
245 0.065834597133618422111563556969398,
246 0.068648672928521619345623411885368,
247 0.071054423553444068305790361723210,
248 0.073030690332786667495189417658913,
249 0.074582875400499188986581418362488,
250 0.075704497684556674659542775376617,
251 0.076377867672080736705502835038061,
252 0.076600711917999656445049901530102
258 0.999262104992609834193457486540341,
259 0.995556969790498097908784946893902,
260 0.988035794534077247637331014577406,
261 0.976663921459517511498315386479594,
262 0.961614986425842512418130033660167,
263 0.942974571228974339414011169658471,
264 0.920747115281701561746346084546331,
265 0.894991997878275368851042006782805,
266 0.865847065293275595448996969588340,
267 0.833442628760834001421021108693570,
268 0.797873797998500059410410904994307,
269 0.759259263037357630577282865204361,
270 0.717766406813084388186654079773298,
271 0.673566368473468364485120633247622,
272 0.626810099010317412788122681624518,
273 0.577662930241222967723689841612654,
274 0.526325284334719182599623778158010,
275 0.473002731445714960522182115009192,
276 0.417885382193037748851814394594572,
277 0.361172305809387837735821730127641,
278 0.303089538931107830167478909980339,
279 0.243866883720988432045190362797452,
280 0.183718939421048892015969888759528,
281 0.122864692610710396387359818808037,
282 0.061544483005685078886546392366797,
283 0.000000000000000000000000000000000
289 0.011393798501026287947902964113235,
290 0.026354986615032137261901815295299,
291 0.040939156701306312655623487711646,
292 0.054904695975835191925936891540473,
293 0.068038333812356917207187185656708,
294 0.080140700335001018013234959669111,
295 0.091028261982963649811497220702892,
296 0.100535949067050644202206890392686,
297 0.108519624474263653116093957050117,
298 0.114858259145711648339325545869556,
299 0.119455763535784772228178126512901,
300 0.122242442990310041688959518945852,
301 0.123176053726715451203902873079050
307 0.001987383892330315926507851882843,
308 0.005561932135356713758040236901066,
309 0.009473973386174151607207710523655,
310 0.013236229195571674813656405846976,
311 0.016847817709128298231516667536336,
312 0.020435371145882835456568292235939,
313 0.024009945606953216220092489164881,
314 0.027475317587851737802948455517811,
315 0.030792300167387488891109020215229,
316 0.034002130274329337836748795229551,
317 0.037116271483415543560330625367620,
318 0.040083825504032382074839284467076,
319 0.042872845020170049476895792439495,
320 0.045502913049921788909870584752660,
321 0.047982537138836713906392255756915,
322 0.050277679080715671963325259433440,
323 0.052362885806407475864366712137873,
324 0.054251129888545490144543370459876,
325 0.055950811220412317308240686382747,
326 0.057437116361567832853582693939506,
327 0.058689680022394207961974175856788,
328 0.059720340324174059979099291932562,
329 0.060539455376045862945360267517565,
330 0.061128509717053048305859030416293,
331 0.061471189871425316661544131965264,
332 0.061580818067832935078759824240066
338 0.999484410050490637571325895705811,
339 0.996893484074649540271630050918695,
340 0.991630996870404594858628366109486,
341 0.983668123279747209970032581605663,
342 0.973116322501126268374693868423707,
343 0.960021864968307512216871025581798,
344 0.944374444748559979415831324037439,
345 0.926200047429274325879324277080474,
346 0.905573307699907798546522558925958,
347 0.882560535792052681543116462530226,
348 0.857205233546061098958658510658944,
349 0.829565762382768397442898119732502,
350 0.799727835821839083013668942322683,
351 0.767777432104826194917977340974503,
352 0.733790062453226804726171131369528,
353 0.697850494793315796932292388026640,
354 0.660061064126626961370053668149271,
355 0.620526182989242861140477556431189,
356 0.579345235826361691756024932172540,
357 0.536624148142019899264169793311073,
358 0.492480467861778574993693061207709,
359 0.447033769538089176780609900322854,
360 0.400401254830394392535476211542661,
361 0.352704725530878113471037207089374,
362 0.304073202273625077372677107199257,
363 0.254636926167889846439805129817805,
364 0.204525116682309891438957671002025,
365 0.153869913608583546963794672743256,
366 0.102806937966737030147096751318001,
367 0.051471842555317695833025213166723,
368 0.000000000000000000000000000000000
374 0.007968192496166605615465883474674,
375 0.018466468311090959142302131912047,
376 0.028784707883323369349719179611292,
377 0.038799192569627049596801936446348,
378 0.048402672830594052902938140422808,
379 0.057493156217619066481721689402056,
380 0.065974229882180495128128515115962,
381 0.073755974737705206268243850022191,
382 0.080755895229420215354694938460530,
383 0.086899787201082979802387530715126,
384 0.092122522237786128717632707087619,
385 0.096368737174644259639468626351810,
386 0.099593420586795267062780282103569,
387 0.101762389748405504596428952168554,
388 0.102852652893558840341285636705415
394 0.001389013698677007624551591226760,
395 0.003890461127099884051267201844516,
396 0.006630703915931292173319826369750,
397 0.009273279659517763428441146892024,
398 0.011823015253496341742232898853251,
399 0.014369729507045804812451432443580,
400 0.016920889189053272627572289420322,
401 0.019414141193942381173408951050128,
402 0.021828035821609192297167485738339,
403 0.024191162078080601365686370725232,
404 0.026509954882333101610601709335075,
405 0.028754048765041292843978785354334,
406 0.030907257562387762472884252943092,
407 0.032981447057483726031814191016854,
408 0.034979338028060024137499670731468,
409 0.036882364651821229223911065617136,
410 0.038678945624727592950348651532281,
411 0.040374538951535959111995279752468,
412 0.041969810215164246147147541285970,
413 0.043452539701356069316831728117073,
414 0.044814800133162663192355551616723,
415 0.046059238271006988116271735559374,
416 0.047185546569299153945261478181099,
417 0.048185861757087129140779492298305,
418 0.049055434555029778887528165367238,
419 0.049795683427074206357811569379942,
420 0.050405921402782346840893085653585,
421 0.050881795898749606492297473049805,
422 0.051221547849258772170656282604944,
423 0.051426128537459025933862879215781,
424 0.051494729429451567558340433647099
429 #ifndef DOXYGEN_NO_O2NS 531 int allocate(
size_t sz);
539 int initialise(
double a,
double b);
544 int set_initial_result(
double result,
double error);
552 int retrieve(
double *a,
double *b,
double *r,
double *e)
const;
572 int update(
double a1,
double b1,
double area1,
double error1,
573 double a2,
double b2,
double area2,
double error2);
576 double sum_results();
581 int subinterval_too_small(
double a1,
double a2,
double b2);
584 void append_interval(
double a1,
double b1,
double area1,
double error1);
616 public inte<func_t> {
734 f_v1=
new double[n_gk];
735 f_v2=
new double[n_gk];
798 template<
class func2_t>
void gauss_kronrod_base
799 (func2_t &func,
double a,
double b,
double *result,
800 double *abserr,
double *resabs,
double *resasc) {
802 const double center=0.5*(a+b);
803 const double half_length=0.5*(b-a);
804 const double abs_half_length=fabs (half_length);
807 f_center=func(center);
809 double result_gauss=0.0;
810 double result_kronrod=f_center*this->w_gk[this->n_gk-1];
812 double result_abs=fabs (result_kronrod);
813 double result_asc=0.0;
814 double mean=0.0, err=0.0;
816 if (this->n_gk % 2 == 0) {
817 result_gauss=f_center*this->w_g[this->n_gk / 2-1];
821 for (
int j=0; j < (this->n_gk-1) / 2; j++) {
824 const double abscissa=half_length*this->x_gk[jtw];
827 fval1=func(center-abscissa);
828 fval2=func(center+abscissa);
830 const double fsum=fval1+fval2;
831 this->f_v1[jtw]=fval1;
832 this->f_v2[jtw]=fval2;
833 result_gauss+=this->w_g[j]*fsum;
834 result_kronrod+=this->w_gk[jtw]*fsum;
835 result_abs+=this->w_gk[jtw]*(fabs(fval1)+fabs(fval2));
839 for (
int j=0; j < this->n_gk / 2; j++) {
842 const double abscissa=half_length*this->x_gk[jtwm1];
845 fval1=func(center-abscissa);
846 fval2=func(center+abscissa);
848 this->f_v1[jtwm1]=fval1;
849 this->f_v2[jtwm1]=fval2;
850 result_kronrod+=this->w_gk[jtwm1]*(fval1+fval2);
851 result_abs+=this->w_gk[jtwm1]*(fabs(fval1)+fabs(fval2));
854 mean=result_kronrod*0.5;
856 result_asc=this->w_gk[this->n_gk-1]*fabs(f_center-mean);
858 for (
int j=0;j<this->n_gk-1;j++) {
859 result_asc+=this->w_gk[j]*(fabs(this->f_v1[j]-mean)+
860 fabs(this->f_v2[j]-mean));
865 err=(result_kronrod-result_gauss)*half_length;
867 result_kronrod*=half_length;
868 result_abs*=abs_half_length;
869 result_asc*=abs_half_length;
871 *result=result_kronrod;
874 *abserr=rescale_error(err,result_abs,result_asc);
881 virtual void gauss_kronrod
882 (func_t &func,
double a,
double b,
883 double *result,
double *abserr,
double *resabs,
double *resasc) {
884 return gauss_kronrod_base<func_t>
885 (func,a,b,result,abserr,resabs,resasc);
891 #ifndef DOXYGEN_NO_O2NS static const double qk51_xgk[26]
Abscissae of the 51-point Kronrod rule.
int set_limit(size_t lim)
Set the limit for the number of subdivisions of the integration region (default 1000) ...
double * f_v1
Scratch space.
static const double qk41_wg[11]
Weights of the 20-point Gauss rule.
Integration workspace for the GSL integrators.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
double * f_v2
Scratch space.
double * alist
Left endpoints of subintervals.
double * blist
Right endpoints of subintervals.
static const double qk31_wg[8]
Weights of the 15-point Gauss rule.
static const double qk41_wgk[21]
Weights of the 41-point Kronrod rule.
size_t * level
Numbers of subdivisions made.
invalid argument supplied by user
int n_gk
Size of Gauss-Kronrod arrays.
static const double qk21_wgk[11]
Weights of the 21-point Kronrod rule.
size_t * order
Linear ordering vector for sort routine.
static const double qk41_xgk[21]
Abscissae of the 41-point Kronrod rule.
size_t i
Index of current subinterval.
static const double qk61_xgk[31]
Abscissae of the 61-point Kronrod rule.
int get_rule()
Get the Gauss-Kronrod integration rule.
static const double qk21_xgk[11]
Abscissae of the 21-point Kronrod rule.
static const double qk15_wgk[8]
Weights of the 15-point Kronrod rule.
static const double qk61_wg[15]
Weights of the 30-point Gauss rule.
Base integration class [abstract base].
double * rlist
Integral approximations on subintervals.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
A namespace for the GSL adaptive Gauss-Kronrod integration coefficients.
int allocate(size_t sz)
Allocate a workspace.
static const double qk61_wgk[31]
Weights of the 61-point Kronrod rule.
static const double qk31_xgk[16]
Abscissae of the 31-point Kronrod rule.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
size_t size
Current number of subintervals being used.
static const double qk15_wg[4]
Weights of the 7-point Gauss rule.
const double * w_gk
Gauss-Kronrod weight pointer.
static const double qk15_xgk[8]
Abscissae of the 15-point Kronrod rule.
double * elist
Integral error estimates.
Basic Gauss-Kronrod integration class (GSL)
size_t nrmax
Counter for extrapolation routine.
static const double qk51_wgk[26]
Weights of the 51-point Kronrod rule.
inte_workspace_gsl * w
The integration workspace.
void set_rule(int rule)
Set the Gauss-Kronrod integration rule to be used.
static const double qk21_wg[5]
Weights of the 10-point Gauss rule.
static const double qk31_wgk[16]
Weights of the 31-point Kronrod rule.
size_t limit
Maximum number of subintervals allocated.
static const double qk51_wg[13]
Weights of the 25-point Gauss rule.
size_t maximum_level
Depth of subdivisions reached.
const double * x_gk
Gauss-Kronrod abscissae pointer.
const double * w_g
Gauss weight pointer.
int free()
Free allocated workspace memory.