boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Merge branch 'gauss_konrod' of https://github.com/NAThompson/math into gauss

Browse Source 
# Please enter a commit message to explain why this merge is necessary,
# especially if it merges an updated upstream into a topic branch.
#
# Lines starting with '#' will be ignored, and an empty message aborts
# the commit.
This commit is contained in:
jzmaddock
2017-08-24 11:47:26 +01:00
parent 37c71f73f2
commit e7e915816f
5 changed files with 297 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/math.qbk
View File

@@ -633,6 +633,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22.
[section:quadrature Quadrature]
[include quadrature/trapezoidal.qbk]
[include quadrature/double_exponential.qbk]
[include quadrature/gauss_kronrod.qbk]
[endsect]
[endmathpart] [/mathpart roots Root Finding Algorithms]

43
doc/quadrature/gauss_kronrod.qbk Normal file
View File

@@ -0,0 +1,43 @@
[/
Copyright (c) 2017 Nick Thompson
Use, modification and distribution are subject to the
Boost Software License, Version 1.0. (See accompanying file
LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
]
[section:gauss_kronrod Gauss-Kronrod Quadrature]
Gauss-Kronrod quadrature is an extension of Gaussian quadrature which provides an a posteriori error estimate for the integral.
The idea behind Gaussian quadrature is to choose /n/ nodes and weights in such a way that polynomials of order /2n-1/ are integrated exactly.
However, integration of polynomials is trivial, so it is rarely done via numerical methods.
Instead, transcendental and numerically defined functions are integrated via Gaussian quadrature, and the defining problem becomes how to estimate the remainder.
Gaussian quadrature alone (without some form on interval splitting) cannot answer this question.
It is possible to compute a Gaussian quadrature of order /n/ and another of order (say) /2n+1/, and use the difference as an error estimate.
However, this is not optimal, as the zeros of the Legendre polynomials (nodes of the Gaussian quadrature) never the same for different orders, so /3n+1/ function evaluations must be performed.
Kronrod considered the problem of how to interleave nodes into a Gaussian quadrature in such a way that all previous function evaluations can be reused, while increasing the order of polynomials that can be integrated exactly.
This allows an a posteriori error estimate to be provided while still preserving exponential convergence.
Kronrod discovered that by adding /n+1/ nodes (computed from the zeros of the Legendre-Stieltjes polynomials) to a Gaussian quadrature of order /n/, he could integrate polynomials of order /3n+1/.
[section:caveats Caveats]
Although this code allows a user to use any floating point type, the nodes and weights have only been computed to 100 decimal digits.
We think this should be acceptable for most use cases; if more digits are needed we recommend the user switch to `tanh_sinh` quadrature,
which is truly arbitrary precision.
For essentially all analytic integrands bounded on the domain, the error estimates provided by the routine are woefully pessimistic.
In fact, in this case the error is very nearly equal to the error of the Gaussian quadrature formula of order /n/,
whereas the Kronrod extension converges exponentially beyond the Gaussian estimate.
Very sophisticated method exist to estimate the error, but all require the integrand to lie in a particular function space.
A more sophisticated a posteriori error estimate for an element of a particular function space is left to the user.
[section:gk_refes References]
* Kronrod, Aleksandr Semenovish (1965), ['Nodes and weights of quadrature formulas. Sixteen-place tables], New York: Consultants Bureau
* Dirk P. Laurie, ['Calculation of Gauss-Kronrod Quadrature Rules], Mathematics of Computation, Volume 66, Number 219, 1997
* Gonnet, Pedro, ['A Review of Error Estimation in Adaptive Quadrature], https://arxiv.org/pdf/1003.4629.pdf
[endsect]
[endsect]

90
example/legendre_stieltjes_example.cpp
View File

@@ -8,6 +8,7 @@
#include <string>
#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/math/special_functions/legendre_stieltjes.hpp>
using boost::math::legendre_p;
@@ -16,83 +17,94 @@ using boost::math::legendre_p_prime;
using boost::math::legendre_stieltjes;
using boost::multiprecision::cpp_bin_float_quad;
using boost::multiprecision::cpp_bin_float_100;
using boost::multiprecision::cpp_dec_float_100;
template<class Real>
void gauss_konrod_rule(size_t order)
void gauss_kronrod_rule(size_t order)
{
std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
std::cout << std::fixed;
auto gauss_nodes = boost::math::legendre_p_zeros<Real>(order);
auto E = legendre_stieltjes<Real>(order + 1);
std::vector<Real> gauss_weights(gauss_nodes.size(), std::numeric_limits<Real>::quiet_NaN());
std::vector<Real> gauss_konrod_weights(gauss_nodes.size(), std::numeric_limits<Real>::quiet_NaN());
std::vector<Real> gauss_kronrod_weights(gauss_nodes.size(), std::numeric_limits<Real>::quiet_NaN());
for (size_t i = 0; i < gauss_nodes.size(); ++i)
{
Real node = gauss_nodes[i];
Real lp = legendre_p_prime<Real>(order, node);
gauss_weights[i] = 2/( (1-node*node)*lp*lp);
// P_n(x) = (2n)!/(2^n (n!)^2) pi_n(x), where pi_n is the monic Legendre polynomial.
gauss_konrod_weights[i] = gauss_weights[i] + static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p_prime(order, node)*E(node));
gauss_kronrod_weights[i] = gauss_weights[i] + static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p_prime(order, node)*E(node));
}
std::cout << "Gauss Nodes:\n";
std::cout << "static const std::vector<Real> gauss_nodes {\n";
for (auto const & node : gauss_nodes)
{
std::cout << node << "\n";
std::cout << " boost::lexical_cast<Real>(\"" << node << "\"),\n";
}
std::cout << "};\n\n";
std::cout << "Gauss Weights:\n";
std::cout << "static const std::vector<Real> gauss_weights {\n";
for (auto const & weight : gauss_weights)
{
std::cout << weight << "\n";
std::cout << " boost::lexical_cast<Real>(\"" << weight << "\"),\n";
}
std::cout << "};\n\n";
std::cout << "static const std::vector<Real> gauss_kronrod_weights {\n";
for (auto const & weight : gauss_kronrod_weights)
{
std::cout << " boost::lexical_cast<Real>(\"" << weight << "\"),\n";
}
std::cout << "};\n\n";
auto kronrod_nodes = E.zeros();
std::vector<Real> kronrod_weights(kronrod_nodes.size());
for (size_t i = 0; i < kronrod_weights.size(); ++i)
{
Real node = kronrod_nodes[i];
kronrod_weights[i] = static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p(order, node)*E.prime(node));
}
std::cout << "Gauss-Konrod weights: \n";
for (auto const & w : gauss_konrod_weights)
std::cout << "static const std::vector<Real> kronrod_nodes {\n";
for (auto node : kronrod_nodes)
{
std::cout << w << "\n";
std::cout << " boost::lexical_cast<Real>(\"" << node << "\"),\n";
}
std::cout << "};\n\n";
auto konrod_nodes = E.zeros();
std::vector<Real> konrod_weights(konrod_nodes.size());
for (size_t i = 0; i < konrod_weights.size(); ++i)
std::cout << "static const std::vector<Real> kronrod_weights {\n";
for (auto const & weight : kronrod_weights)
{
Real node = konrod_nodes[i];
konrod_weights[i] = static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p(order, node)*E.prime(node));
}
std::cout << "Konrod nodes:\n";
for (auto node : konrod_nodes)
{
std::cout << node << "\n";
}
std::cout << "Konrod weights: \n";
for (auto const & w : gauss_konrod_weights)
{
std::cout << w << "\n";
std::cout << " boost::lexical_cast<Real>(\"" << weight << "\"),\n";
}
std::cout << "};\n\n";
}
int main()
{
std::cout << "Gauss-Konrod 7-15 Rule:\n";
gauss_konrod_rule<cpp_bin_float_100>(7);
//std::cout << "Gauss-Kronrod 7-15 Rule:\n";
//gauss_kronrod_rule<cpp_dec_float_100>(7);
std::cout << "\n\nGauss-Konrod 10-21 Rule:\n";
gauss_konrod_rule<cpp_bin_float_100>(10);
//std::cout << "\n\nGauss-Kronrod 10-21 Rule:\n";
//gauss_kronrod_rule<cpp_dec_float_100>(10);
std::cout << "\n\nGauss-Konrod 15-31 Rule:\n";
gauss_konrod_rule<cpp_bin_float_100>(15);
std::cout << "\n\nGauss-Kronrod 15-31 Rule:\n";
gauss_kronrod_rule<cpp_dec_float_100>(15);
/*
std::cout << "\n\nGauss-Kronrod 20-41 Rule:\n";
gauss_kronrod_rule<cpp_dec_float_100>(20);
std::cout << "\n\nGauss-Konrod 20-41 Rule:\n";
gauss_konrod_rule<cpp_bin_float_100>(20);
std::cout << "\n\nGauss-Kronrod 25-51 Rule:\n";
gauss_kronrod_rule<cpp_dec_float_100>(25);
std::cout << "\n\nGauss-Konrod 25-51 Rule:\n";
gauss_konrod_rule<cpp_bin_float_100>(25);
std::cout << "\n\nGauss-Kronrod 30-61 Rule:\n";
gauss_kronrod_rule<cpp_dec_float_100>(30);
std::cout << "\n\nGauss-Konrod 30-61 Rule:\n";
gauss_konrod_rule<cpp_bin_float_100>(30);
std::cout << "\n\nGauss-Kronrod 35-71 Rule:\n";
gauss_kronrod_rule<cpp_dec_float_100>(35);
std::cout << "\n\nGauss-Kronrod 40-81 Rule:\n";
gauss_kronrod_rule<cpp_dec_float_100>(40);*/
}

142
include/boost/math/quadrature/detail/gauss_kronrod_detail.hpp Normal file
View File

@@ -0,0 +1,142 @@
// Copyright Nick Thompson 2017.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <vector>
namespace boost { namespace math { namespace quadrature { namespace detail {
// These nodes and weights can be verified at:
// https://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/
template<class Real, class F>
Real g7_k15(const F& f, Real* error = nullptr)
{
static const std::vector<Real> gauss_nodes{
boost::lexical_cast<Real>("0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
boost::lexical_cast<Real>("0.4058451513773971669066064120769614633473820140993701263870432517946638132261256553283126897277465878"),
boost::lexical_cast<Real>("0.7415311855993944398638647732807884070741476471413902601199553519674298746721805137928268323668632471"),
boost::lexical_cast<Real>("0.9491079123427585245261896840478512624007709376706177835487691039130633303548401408057307700279257241")
};
static const std::vector<Real> gauss_weights {
boost::lexical_cast<Real>("0.4179591836734693877551020408163265306122448979591836734693877551020408163265306122448979591836734694"),
boost::lexical_cast<Real>("0.3818300505051189449503697754889751338783650835338627347510834510307055464341297083486846593440448015"),
boost::lexical_cast<Real>("0.2797053914892766679014677714237795824869250652265987645370140326936188104305626768132409429011976188"),
boost::lexical_cast<Real>("0.1294849661688696932706114326790820183285874022599466639772086387246552349720423087156254181629208451")
};
static const std::vector<Real> gauss_kronrod_weights {
boost::lexical_cast<Real>("0.2094821410847278280129991748917142636977620802237043167129980065613751513232564861681690821167594910"),
boost::lexical_cast<Real>("0.1903505780647854099132564024210136828260780754553583558854408803674405807241021267960596460510637759"),
boost::lexical_cast<Real>("0.1406532597155259187451895905102379203998897572479985755617454689331270809309095040809737912241555591"),
boost::lexical_cast<Real>("0.0630920926299785532907006631892042866650711572115507071136055451469839974779648749281991702645044420"),
};
static const std::vector<Real> kronrod_nodes {
boost::lexical_cast<Real>("0.2077849550078984676006894037732449134797844071451706497138457346198669384494352022691034322718369853"),
boost::lexical_cast<Real>("0.5860872354676911302941448382587295984367807506043609513049928931988037360744440746451167449893594210"),
boost::lexical_cast<Real>("0.8648644233597690727897127886409262012109723070740881486014577127670677081325957210358584785960459054"),
boost::lexical_cast<Real>("0.9914553711208126392068546975263285166420443383703347012910874135724417393465340723592450350962684176")
};
static const std::vector<Real> kronrod_weights {
boost::lexical_cast<Real>("0.2044329400752988924141619992346490847165176041807183574244709531204546769854659887934837429200934755"),
boost::lexical_cast<Real>("0.1690047266392679028265834265985502841062449003029442414973400675569568092161902911293670240385535991"),
boost::lexical_cast<Real>("0.1047900103222501838398763225415180174437566542138306118933906513396374632157628952416757162750931133"),
boost::lexical_cast<Real>("0.0229353220105292249637320080589695919935608112757469922675074302547118157879760759461563681681562895"),
};
Real f0 = f(0);
Real I0 = f0*gauss_weights[0];
Real I1 = f0*gauss_kronrod_weights[0];
for (size_t i = 1; i < gauss_nodes.size(); ++i) {
Real yp = f(gauss_nodes[i]);
Real ym = f(-gauss_nodes[i]);
I0 += gauss_weights[i]*(yp + ym);
I1 += gauss_kronrod_weights[i]*(yp + ym);
}
for(size_t i = 0; i < kronrod_weights.size(); ++i) {
I1 += kronrod_weights[i]*(f(kronrod_nodes[i]) + f(-kronrod_nodes[i]));
}
if (error) {
using std::abs;
*error = abs(I1 - I0);
}
return I1;
}
template<class Real, class F>
Real g10_k21(const F& f, Real* error = nullptr)
{
static const std::vector<Real> gauss_nodes {
boost::lexical_cast<Real>("0.1488743389816312108848260011297199846175648594206916957079892535159036173556685213711776297994636912"),
boost::lexical_cast<Real>("0.4333953941292471907992659431657841622000718376562464965027015131437669890777035012251027579501177212"),
boost::lexical_cast<Real>("0.6794095682990244062343273651148735757692947118348094676648171889525585753950749246150785735704803795"),
boost::lexical_cast<Real>("0.8650633666889845107320966884234930485275430149653304525219597318453747551380555613567907289460457707"),
boost::lexical_cast<Real>("0.9739065285171717200779640120844520534282699466923821192312120666965952032346361596257235649562685563"),
};
static const std::vector<Real> gauss_weights {
boost::lexical_cast<Real>("0.2955242247147528701738929946513383294210467170268536013543080297559959382171523292703565957937542167"),
boost::lexical_cast<Real>("0.2692667193099963550912269215694693528597599384608837958005632762421534323191792767642266367092527608"),
boost::lexical_cast<Real>("0.2190863625159820439955349342281631924587718705226770898809565436351999106529512812426839931772021928"),
boost::lexical_cast<Real>("0.1494513491505805931457763396576973324025566396694273678354772687532386547266300109459472646347319519"),
boost::lexical_cast<Real>("0.0666713443086881375935688098933317928578648343201581451286948816134120640840871017767855096850588778"),
};
static const std::vector<Real> gauss_kronrod_weights {
boost::lexical_cast<Real>("0.1477391049013384913748415159720680455237316254852066045181919543988599301673569640573270395918288225"),
boost::lexical_cast<Real>("0.1347092173114733259280540017717068327609919130085597140663666849132029140012128203667695315948827177"),
boost::lexical_cast<Real>("0.1093871588022976418992105903258049602718132998343452200781967582982655037289143216868389943267455384"),
boost::lexical_cast<Real>("0.0750396748109199527670431409161900093952193820009100881736970480484304043428584951788138087306465541"),
boost::lexical_cast<Real>("0.0325581623079647274788189724593897606173889398456626095715375042327141218201654986923816076053846265"),
};
static const std::vector<Real> kronrod_nodes {
boost::lexical_cast<Real>("0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
boost::lexical_cast<Real>("0.2943928627014601981311266031038655661626866251569579186488822917272461116633273788844552317826823736"),
boost::lexical_cast<Real>("0.5627571346686046833390000992726941408430138819419669588603462145877926635321632754971208785416999242"),
boost::lexical_cast<Real>("0.7808177265864168970637175783450423771634075202981571797469485999950560798276142065452697723423899624"),
boost::lexical_cast<Real>("0.9301574913557082260012071800595083462251679099819392423034940686682841598309167305501119457285100788"),
boost::lexical_cast<Real>("0.9956571630258080807355272806890028479212605872194789243633791611175702304677486735715232599691207672"),
};
static const std::vector<Real> kronrod_weights {
boost::lexical_cast<Real>("0.1494455540029169056649364683898212037452363166874728038356085187369896447851184192572103070568954026"),
boost::lexical_cast<Real>("0.1427759385770600807970942731387170608859790565319055556074100474397077044990934002781113170628375643"),
boost::lexical_cast<Real>("0.1234919762620658510779581098310741595123003495286483276446799412097405423897545468968153862236373823"),
boost::lexical_cast<Real>("0.0931254545836976055350654650833663443900188288807600319700850387601777356722007752374141230616158275"),
boost::lexical_cast<Real>("0.0547558965743519960313813002445801763737211140583335575244326158047840989278189753251163015690032981"),
boost::lexical_cast<Real>("0.0116946388673718742780643960621920483962173324819318889275981475256222220580649926518067367049699673"),
};
Real I0 = 0;
Real I1 = 0;
for (size_t i = 0; i < gauss_nodes.size(); ++i) {
Real yp = f(gauss_nodes[i]);
Real ym = f(-gauss_nodes[i]);
I0 += gauss_weights[i]*(yp + ym);
I1 += gauss_kronrod_weights[i]*(yp + ym);
}
I1 += f(0)*kronrod_weights[0];
for(size_t i = 1; i < kronrod_weights.size(); ++i) {
I1 += kronrod_weights[i]*(f(kronrod_nodes[i]) + f(-kronrod_nodes[i]));
}
if (error) {
using std::abs;
*error = abs(I1 - I0);
}
return I1;
}
}}}}

60
test/gauss_kronrod_test.cpp Normal file
View File

@@ -0,0 +1,60 @@
// Copyright Nick Thompson 2017.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MODULE GaussKronrod
#include <boost/test/unit_test.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/math/quadrature/detail/gauss_kronrod_detail.hpp>
using boost::math::constants::half_pi;
using boost::math::quadrature::detail::g7_k15;
using boost::math::quadrature::detail::g10_k21;
using boost::multiprecision::cpp_bin_float_quad;
template<class Real>
void test_g7_k15()
{
std::cout << std::setprecision(std::numeric_limits<Real>::digits10) << std::fixed;
auto f = [](Real x) { return 1/(1+x*x); };
Real error;
Real Q = g7_k15<Real, decltype(f)>(f, &error);
std::cout << Q << std::endl;
std::cout << "Q - pi/2 = " << Q - half_pi<Real>() << std::endl;
std::cout << "Error estimate: " << error << std::endl;
}
template<class Real>
void test_g10_k21()
{
std::cout << std::setprecision(std::numeric_limits<Real>::digits10) << std::fixed;
auto f = [](Real x) { return 1/(1+x*x); };
Real error;
Real Q = g10_k21<Real, decltype(f)>(f, &error);
std::cout << Q << std::endl;
std::cout << "Q - pi/2 = " << Q - half_pi<Real>() << std::endl;
std::cout << "Error estimate: " << error << std::endl;
}
BOOST_AUTO_TEST_CASE(GaussKronrod)
{
test_g7_k15<double>();
test_g7_k15<long double>();
test_g7_k15<cpp_bin_float_quad>();
test_g7_k15<boost::multiprecision::cpp_bin_float_100>();
test_g10_k21<double>();
test_g10_k21<long double>();
test_g10_k21<cpp_bin_float_quad>();
test_g10_k21<boost::multiprecision::cpp_bin_float_100>();
}