gammalib / gammalib

Commit 6393b152cbaa92af8c436ed78f7631f401147ddb

Authored by Jürgen Knödlseder
1 parent 9c93b03c

Add GIntegral::adaptive_gauss_kronrod() method (#4204)

Inline Side-by-side
Showing 9 changed files with 275 additions and 27 deletions
ChangeLog
  View file @6393b15
1   -2023-01-19
  1 +2023-01-20
2 2  
3 3 * Version 2.1.0 released
4 4 ========================
5 5  
  6 + Add GIntegral::adaptive_gauss_kronrod() method (#4204)
6 7 Add optional unit parameter to GApplication::log_value() methods (#4202)
7 8 Fix segfault on saving empty GModelSpectralTable (#4198)
8 9 Implement response vector cache storage in FITS file (#4159)
... ...
1 1 New Features and Important Changes in GammaLib 2.1.0
2 2  
3   -19 January 2023
  3 +20 January 2023
4 4  
5 5  
6 6 1. Introduction
... ... @@ -22,6 +22,7 @@ The following classes have been renamed:
22 22  
23 23 The following methods have been added:
24 24 - GNdarray::index(int&)
  25 +- GIntegral::adaptive_gauss_kronrod()
25 26 - GModelSpectralTable::nspectra()
26 27 - GModelSpectralTable::scale_energy()
27 28 - GModelSpectralTable::has_energy_scale()
... ... @@ -92,6 +93,9 @@ table model (#4069).
92 93  
93 94 10. Numerics module
94 95 -------------------
  96 +Added GIntegral::adaptive_gauss_kronrod() method to adaptively compute the
  97 +Gauss-Lobatto-Kronrod integral of a kernel function (#4204).
  98 +
95 99 Added GNdarray::index(int&) method to return the index vector from an array element
96 100 index (#4069).
97 101  
... ...
README.md
  View file @6393b15
1 1 GammaLib information
2 2 ====================
3   -* Version: 2.1.0.dev (19 January 2023)
  3 +* Version: 2.1.0.dev (20 January 2023)
4 4  
5 5 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
6 6  
... ...
doc/source/admin/release_history/2.1.rst
  View file @6393b15
... ... @@ -23,6 +23,8 @@ Bug fixes
23 23 Improvements
24 24 ------------
25 25  
  26 +* [`4204 <https://cta-redmine.irap.omp.eu/issues/4204>`_] -
  27 + Add ``GIntegral::adaptive_gauss_kronrod()`` method
26 28 * [`4202 <https://cta-redmine.irap.omp.eu/issues/4202>`_] -
27 29 Add optional unit parameter to ``GApplication::log_value()`` methods
28 30 * [`4159 <https://cta-redmine.irap.omp.eu/issues/4159>`_] -
... ...
include/GIntegral.hpp
  View file @6393b15
1 1 /***************************************************************************
2 2 * GIntegral.hpp - Integration class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2023 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -83,6 +83,7 @@ public:
83 83 const int& n = 1,
84 84 double result = 0.0);
85 85 double adaptive_simpson(const double& a, const double& b) const;
  86 + double adaptive_gauss_kronrod(const double& a, const double& b) const;
86 87 double gauss_kronrod(const double& a, const double& b) const;
87 88 std::string print(const GChatter& chatter = NORMAL) const;
88 89  
... ... @@ -97,6 +98,10 @@ protected:
97 98 const double& fa, const double& fb,
98 99 const double& fc,
99 100 const int& bottom) const;
  101 + double adaptive_gauss_kronrod_aux(const double& a, const double& b,
  102 + const double& fa, const double& fb,
  103 + const double& is,
  104 + const double& toler) const;
100 105 double rescale_error(double err,
101 106 const double& result_abs,
102 107 const double& result_asc) const;
... ... @@ -117,6 +122,7 @@ protected:
117 122 mutable double m_abserr; //!< Absolute integration error
118 123 mutable double m_relerr; //!< Absolute integration error
119 124 mutable std::string m_message; //!< Status message (if result is invalid)
  125 + mutable bool m_terminate; //!< Signals termination of subdivision
120 126 };
121 127  
122 128  
... ...
pyext/GIntegral.i
  View file @6393b15
1 1 /***************************************************************************
2 2 * GIntegral.i - Integration class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2023 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -71,6 +71,7 @@ public:
71 71 double trapzd(const double& a, const double& b,
72 72 const int& n = 1, double result = 0.0);
73 73 double adaptive_simpson(const double& a, const double& b) const;
  74 + double adaptive_gauss_kronrod(const double& a, const double& b) const;
74 75 double gauss_kronrod(const double& a, const double& b) const;
75 76 };
76 77  
... ...
src/numerics/GIntegral.cpp
  View file @6393b15
1 1 /***************************************************************************
2 2 * GIntegral.cpp - Integration class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2023 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -28,6 +28,7 @@
28 28 #include <cmath> // std::abs()
29 29 #include <vector>
30 30 #include <algorithm> // std::sort
  31 +#include <limits> // numeric_limits
31 32 #include "GIntegral.hpp"
32 33 #include "GException.hpp"
33 34 #include "GTools.hpp"
... ... @@ -433,7 +434,7 @@ double GIntegral::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
433 434 m_calls = 0;
434 435 m_has_abserr = false;
435 436 m_has_relerr = false;
436   -
  437 +
437 438 // Continue only if integration range is valid
438 439 if (b > a) {
439 440  
... ... @@ -541,7 +542,7 @@ double GIntegral::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
541 542 gammalib::warning(origin, m_message);
542 543 }
543 544 }
544   -
  545 +
545 546 } // endif: integration range was valid
546 547  
547 548 // Return result
... ... @@ -596,21 +597,21 @@ double GIntegral::trapzd(const double&amp; a, const double&amp; b, const int&amp; n,
596 597 if (a == b) {
597 598 result = 0.0;
598 599 }
599   -
  600 +
600 601 // ... otherwise use trapeziodal rule
601 602 else {
602   -
  603 +
603 604 // Case A: Only a single step is requested
604 605 if (n == 1) {
605   -
  606 +
606 607 // Evaluate integrand at boundaries
607 608 double y_a = m_kernel->eval(a);
608 609 double y_b = m_kernel->eval(b);
609 610 m_calls += 2;
610   -
  611 +
611 612 // Compute result
612 613 result = 0.5*(b-a)*(y_a + y_b);
613   -
  614 +
614 615 } // endif: only a single step was requested
615 616  
616 617 // Case B: More than a single step is requested
... ... @@ -660,20 +661,20 @@ double GIntegral::trapzd(const double&amp; a, const double&amp; b, const int&amp; n,
660 661 double x = a + 0.5*del;
661 662 double sum = 0.0;
662 663 for (int j = 0; j < it; ++j, x+=del) {
663   -
  664 +
664 665 // Evaluate integrand
665 666 double y = m_kernel->eval(x);
666 667 m_calls++;
667 668  
668 669 // Add integrand
669 670 sum += y;
670   -
  671 +
671 672 } // endfor: looped over steps
672 673  
673 674 // Set result
674 675 result = 0.5*(result + (b-a)*sum/tnm);
675 676 }
676   -
  677 +
677 678 } // endelse: trapeziodal rule was applied
678 679  
679 680 // Return result
... ... @@ -723,7 +724,7 @@ double GIntegral::adaptive_simpson(const double&amp; a, const double&amp; b) const
723 724 double fb = m_kernel->eval(b);
724 725 double fc = m_kernel->eval(c);
725 726 m_calls += 3;
726   -
  727 +
727 728 // Compute integral using Simpson's rule
728 729 double S = (h/6.0) * (fa + 4.0*fc + fb);
729 730  
... ... @@ -757,6 +758,99 @@ double GIntegral::adaptive_simpson(const double&amp; a, const double&amp; b) const
757 758  
758 759  
759 760 /***********************************************************************//**
  761 + * @brief Adaptive Gauss-Lobatto-Kronrod integration
  762 + *
  763 + * @param[in] a Left integration boundary.
  764 + * @param[in] b Right integration boundary.
  765 + * @return Integral of kernel over interval [a,b].
  766 + *
  767 + * Integrates the function using an adaptive Gauss-Lobatto method with a
  768 + * Kronrod extension.
  769 + ***************************************************************************/
  770 +double GIntegral::adaptive_gauss_kronrod(const double& a, const double& b) const
  771 +{
  772 + // Set constants
  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();
  780 +
  781 + // Set tolerance, assuring that the tolerance is not smaller than the
  782 + // machine precision
  783 + double tol = (m_eps < 10.0*eps) ? 10.0 * eps : m_eps;
  784 +
  785 + // Initialise integration status information
  786 + m_isvalid = true;
  787 + m_iter = 0;
  788 + m_calls = 0;
  789 + m_has_abserr = false;
  790 + m_has_relerr = false;
  791 + m_terminate = true;
  792 +
  793 + // Allocate result storage array
  794 + double y[13];
  795 +
  796 + // Compute midpoint and step size
  797 + double m = 0.5 * (a+b);
  798 + double h = 0.5 * (b-a);
  799 +
  800 + // Evaluate function at end points
  801 + double fa = m_kernel->eval(a);
  802 + double fb = m_kernel->eval(b);
  803 + m_calls += 2;
  804 +
  805 + // Store end points
  806 + y[0] = fa;
  807 + y[12] = fb;
  808 +
  809 + // Evaluate function at intermediate points
  810 + for (int i = 1; i < 12; ++i) {
  811 + y[i] = m_kernel->eval(m + x[i]*h);
  812 + m_calls++;
  813 + }
  814 +
  815 + // 4-point Gauss-Lobatto formula
  816 + double i2 = (h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));
  817 +
  818 + // 7-point Konrod extension
  819 + double i1 = (h/1470.0) * ( 77.0 * (y[0]+y[12])+
  820 + 432.0 * (y[2]+y[10])+
  821 + 625.0 * (y[4]+y[8])+
  822 + 672.0 * y[6]);
  823 +
  824 + // 13-point Konrod extension
  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]);
  832 +
  833 + // Estimate errors
  834 + double erri1 = std::abs(i1-is);
  835 + double erri2 = std::abs(i2-is);
  836 +
  837 + // Check on convergence
  838 + double r = (erri2 != 0.0) ? erri1/erri2 : 1.0;
  839 + double toler = (r > 0.0 && r < 1.0) ? tol/r : tol;
  840 + if (is == 0.0) {
  841 + is = b-a;
  842 + }
  843 +
  844 + // Call helper
  845 + is = std::abs(is);
  846 + double value = adaptive_gauss_kronrod_aux(a, b, fa, fb, is, toler);
  847 +
  848 + // Return result
  849 + return value;
  850 +}
  851 +
  852 +
  853 +/***********************************************************************//**
760 854 * @brief Gauss-Kronrod integration
761 855 *
762 856 * @param[in] a Left integration boundary.
... ... @@ -929,9 +1023,9 @@ double GIntegral::gauss_kronrod(const double&amp; a, const double&amp; b) const
929 1023 gammalib::str(b)+")";
930 1024 gammalib::warning(origin, m_message);
931 1025 }
932   -
  1026 +
933 1027 } while (false); // end of main loop
934   -
  1028 +
935 1029 // Return result
936 1030 return result;
937 1031 }
... ... @@ -1191,7 +1285,7 @@ double GIntegral::adaptive_simpson_aux(const double&amp; a, const double&amp; b,
1191 1285 double fd = m_kernel->eval(d);
1192 1286 double fe = m_kernel->eval(e);
1193 1287 m_calls += 2;
1194   -
  1288 +
1195 1289 // Compute integral using Simpson's rule for the left and right interval
1196 1290 double h12 = h / 12.0;
1197 1291 double Sleft = h12 * (fa + 4.0*fd + fc);
... ... @@ -1206,14 +1300,14 @@ double GIntegral::adaptive_simpson_aux(const double&amp; a, const double&amp; b,
1206 1300 // if (std::abs(S2 - S) <= 15.0 * m_eps * std::abs(S2)) {
1207 1301 value = S2 + (S2 - S)/15.0;
1208 1302 }
1209   -
  1303 +
1210 1304 // ... else if the maximum recursion depth was reached then compute the
1211 1305 // result and signal result invalidity
1212 1306 else if (bottom <= 0) {
1213 1307 value = S2 + (S2 - S)/15.0;
1214 1308 m_isvalid = false;
1215 1309 }
1216   -
  1310 +
1217 1311 // ... otherwise call this method recursively
1218 1312 else {
1219 1313 value = adaptive_simpson_aux(a, c, 0.5*eps, Sleft, fa, fc, fd, bottom-1) +
... ... @@ -1226,6 +1320,81 @@ double GIntegral::adaptive_simpson_aux(const double&amp; a, const double&amp; b,
1226 1320  
1227 1321  
1228 1322 /***********************************************************************//**
  1323 + * @brief Adaptive Gauss-Lobatto-Kronrod integration helper
  1324 + *
  1325 + * @param[in] a Left integration boundary.
  1326 + * @param[in] b Right integration boundary.
  1327 + * @param[in] fa Function value at the left integration boundary.
  1328 + * @param[in] fb Function value at the right integration boundary.
  1329 + * @param[in] is 13-point Kronrod extension of previous step.
  1330 + * @return Integral of kernel over interval [a,b].
  1331 + *
  1332 + * Integrates the function using an adaptive Gauss-Lobatto method with a
  1333 + * Kronrod extension.
  1334 + ***************************************************************************/
  1335 +double GIntegral::adaptive_gauss_kronrod_aux(const double& a,
  1336 + const double& b,
  1337 + const double& fa,
  1338 + const double& fb,
  1339 + const double& is,
  1340 + const double& toler) const
  1341 +{
  1342 + // Set constants
  1343 + const double alpha = std::sqrt(2.0/3.0);
  1344 + const double beta = 1.0/std::sqrt(5.0);
  1345 +
  1346 + // Initialise value
  1347 + double value = 0.0;
  1348 +
  1349 + // Compute midpoint and step size
  1350 + double m = 0.5*(a+b);
  1351 + double h = 0.5*(b-a);
  1352 +
  1353 + // Compute intermediate points
  1354 + double mll = m - alpha*h;
  1355 + double ml = m - beta*h;
  1356 + double mr = m + beta*h;
  1357 + double mrr = m + alpha*h;
  1358 +
  1359 + // Evaluate function at intermediate points
  1360 + double fmll = m_kernel->eval(mll);
  1361 + double fml = m_kernel->eval(ml);
  1362 + double fm = m_kernel->eval(m);
  1363 + double fmr = m_kernel->eval(mr);
  1364 + double fmrr = m_kernel->eval(mrr);
  1365 + m_calls += 5;
  1366 +
  1367 + // 4-point Gauss-Lobatto formula
  1368 + double i2 = h/6.0*(fa+fb+5.0*(fml+fmr));
  1369 +
  1370 + // 7-point Konrod extension
  1371 + double i1 = h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);
  1372 +
  1373 + // Check for convergence
  1374 + if (std::abs(i1-i2) <= toler*is || mll <= a || b <= mrr) {
  1375 + if ((mll <= a || b <= mrr) && m_terminate) {
  1376 + m_isvalid = false;
  1377 + m_terminate = false;
  1378 + }
  1379 + value = i1;
  1380 + }
  1381 +
  1382 + // ... otherwise subdivide interval
  1383 + else {
  1384 + value = adaptive_gauss_kronrod_aux(a, mll, fa, fmll, is, toler) +
  1385 + adaptive_gauss_kronrod_aux(mll, ml, fmll, fml, is, toler) +
  1386 + adaptive_gauss_kronrod_aux(ml, m, fml, fm, is, toler) +
  1387 + adaptive_gauss_kronrod_aux(m, mr, fm, fmr, is, toler) +
  1388 + adaptive_gauss_kronrod_aux(mr, mrr, fmr, fmrr, is, toler) +
  1389 + adaptive_gauss_kronrod_aux(mrr, b, fmrr, fb, is, toler);
  1390 + }
  1391 +
  1392 + // Return result
  1393 + return value;
  1394 +}
  1395 +
  1396 +
  1397 +/***********************************************************************//**
1229 1398 * @brief Rescale errors for Gauss-Kronrod integration
1230 1399 *
1231 1400 * @param[in] err Error estimate.
... ...
test/test_GNumerics.cpp
  View file @6393b15
1 1 /***************************************************************************
2 2 * test_GNumerics.cpp - test numerics modules *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2020 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2023 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -66,8 +66,12 @@ void TestGNumerics::set(void)
66 66 "Test GIntegrals");
67 67 append(static_cast<pfunction>(&TestGNumerics::test_romberg_integration),
68 68 "Test Romberg integration");
  69 + append(static_cast<pfunction>(&TestGNumerics::test_trapzd_integration),
  70 + "Test Trapezoidal integration");
69 71 append(static_cast<pfunction>(&TestGNumerics::test_adaptive_simpson_integration),
70 72 "Test adaptive Simpson integration");
  73 + append(static_cast<pfunction>(&TestGNumerics::test_adaptive_gauss_kronrod_integration),
  74 + "Test adaptive Gauss-Kronrod integration");
71 75 append(static_cast<pfunction>(&TestGNumerics::test_gauss_kronrod_integration),
72 76 "Test Gauss-Kronrod integration");
73 77  
... ... @@ -577,6 +581,8 @@ void TestGNumerics::test_integral(void)
577 581 // Test integral and integrand allocation
578 582 Gauss integrand(m_sigma);
579 583 GIntegral integral(&integrand);
  584 + test_value(integral.iter(), 0, "Check initial iterations");
  585 + test_value(integral.calls(), 0, "Check initial calls");
580 586  
581 587 // Exit test
582 588 return;
... ... @@ -631,7 +637,7 @@ void TestGNumerics::test_romberg_integration(void)
631 637  
632 638 // Integrate over the entire Gaussian
633 639 double result = integral.romberg(-10.0*m_sigma, 10.0*m_sigma);
634   - test_value(result, 1.0, 1.0e-6, "Check full integration result");
  640 + test_value(result, 1.0, 1.0e-6, "Check full interval integration result");
635 641  
636 642 // Test [-1sigma, 1sigma]
637 643 result = integral.romberg(-m_sigma, m_sigma);
... ... @@ -649,6 +655,34 @@ void TestGNumerics::test_romberg_integration(void)
649 655  
650 656  
651 657 /***********************************************************************//**
  658 + * @brief Test Trapezoidal integration
  659 + ***************************************************************************/
  660 +void TestGNumerics::test_trapzd_integration(void)
  661 +{
  662 + // Set-up integral
  663 + Gauss integrand(m_sigma);
  664 + GIntegral integral(&integrand);
  665 +
  666 + // Integrate over the entire Gaussian
  667 + double result = integral.trapzd(-10.0*m_sigma, 10.0*m_sigma, 100);
  668 + test_value(result, 1.0, 1.0e-6, "Check full interval integration result");
  669 +
  670 + // Test [-1sigma, 1sigma]
  671 + result = integral.trapzd(-m_sigma, m_sigma, 100);
  672 + test_value(result, 0.68268948, 1.0e-6,
  673 + "Check [-1sigma, 1sigma] integration results");
  674 +
  675 + // Test [0.0, 1sigma]
  676 + result = integral.trapzd(0.0, m_sigma, 100);
  677 + test_value(result, 0.34134475, 1.0e-6,
  678 + "Check [0, 1sigma] integration results");
  679 +
  680 + // Return
  681 + return;
  682 +}
  683 +
  684 +
  685 +/***********************************************************************//**
652 686 * @brief Test adaptive Simpson integration
653 687 ***************************************************************************/
654 688 void TestGNumerics::test_adaptive_simpson_integration(void)
... ... @@ -659,7 +693,7 @@ void TestGNumerics::test_adaptive_simpson_integration(void)
659 693  
660 694 // Integrate over the entire Gaussian
661 695 double result = integral.adaptive_simpson(-10.0*m_sigma, 10.0*m_sigma);
662   - test_value(result, 1.0, 1.0e-6, "Check full integration result");
  696 + test_value(result, 1.0, 1.0e-6, "Check full interval integration result");
663 697  
664 698 // Test [-1sigma, 1sigma]
665 699 result = integral.adaptive_simpson(-m_sigma, m_sigma);
... ... @@ -678,6 +712,35 @@ void TestGNumerics::test_adaptive_simpson_integration(void)
678 712  
679 713  
680 714 /***********************************************************************//**
  715 + * @brief Test adaptive Gauss-Kronrod integration
  716 + ***************************************************************************/
  717 +void TestGNumerics::test_adaptive_gauss_kronrod_integration(void)
  718 +{
  719 + // Set-up integral
  720 + Gauss integrand(m_sigma);
  721 + GIntegral integral(&integrand);
  722 +
  723 + // Integrate over the entire Gaussian
  724 + double result = integral.adaptive_gauss_kronrod(-10.0*m_sigma, 10.0*m_sigma);
  725 + test_value(result, 1.0, 1.0e-6, "Check full interval integration result");
  726 +
  727 + // Test [-1sigma, 1sigma]
  728 + result = integral.adaptive_gauss_kronrod(-m_sigma, m_sigma);
  729 + test_value(result, 0.68268948, 1.0e-6,
  730 + "Check [-1sigma, 1sigma] integration results");
  731 +
  732 +
  733 + // Test [0.0, 1sigma]
  734 + result = integral.adaptive_gauss_kronrod(0.0, m_sigma);
  735 + test_value(result, 0.34134475, 1.0e-6,
  736 + "Check [0, 1sigma] integration results");
  737 +
  738 + // Return
  739 + return;
  740 +}
  741 +
  742 +
  743 +/***********************************************************************//**
681 744 * @brief Test Gauss-Kronrod integration
682 745 ***************************************************************************/
683 746 void TestGNumerics::test_gauss_kronrod_integration(void)
... ... @@ -688,7 +751,7 @@ void TestGNumerics::test_gauss_kronrod_integration(void)
688 751  
689 752 // Integrate over the entire Gaussian
690 753 double result = integral.gauss_kronrod(-10.0*m_sigma, 10.0*m_sigma);
691   - test_value(result, 1.0, 1.0e-6, "Check full integration result");
  754 + test_value(result, 1.0, 1.0e-6, "Check full interval integration result");
692 755  
693 756 // Test [-1sigma, 1sigma]
694 757 result = integral.gauss_kronrod(-m_sigma, m_sigma);
... ...
test/test_GNumerics.hpp
  View file @6393b15
1 1 /***************************************************************************
2 2 * test_GNumerics.hpp - test numerics modules *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2023 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -105,7 +105,9 @@ public:
105 105 void test_integral(void);
106 106 void test_integrals(void);
107 107 void test_romberg_integration(void);
  108 + void test_trapzd_integration(void);
108 109 void test_adaptive_simpson_integration(void);
  110 + void test_adaptive_gauss_kronrod_integration(void);
109 111 void test_gauss_kronrod_integration(void);
110 112  
111 113 private:
... ...