C++：实现量化Integration积分测试实例

#include "integrals.hpp"
#include "utilities.hpp"
#include <ql/math/integrals/exponentialintegrals.hpp>
#include <ql/math/integrals/filonintegral.hpp>
#include <ql/math/integrals/segmentintegral.hpp>
#include <ql/math/integrals/simpsonintegral.hpp>
#include <ql/math/integrals/trapezoidintegral.hpp>
#include <ql/math/integrals/kronrodintegral.hpp>
#include <ql/math/integrals/gausslobattointegral.hpp>
#include <ql/math/integrals/discreteintegrals.hpp>
#include <ql/math/integrals/tanhsinhintegral.hpp>
#include <ql/math/integrals/gaussianquadratures.hpp>
#include <ql/math/interpolations/bilinearinterpolation.hpp>
#include <ql/math/distributions/normaldistribution.hpp>
#include <ql/termstructures/volatility/abcd.hpp>
#include <ql/math/integrals/twodimensionalintegral.hpp>
#include <ql/experimental/math/piecewisefunction.hpp>
#include <ql/experimental/math/piecewiseintegral.hpp>using namespace QuantLib;
using boost::unit_test_framework::test_suite;namespace integrals_test {Real tolerance = 1.0e-6;template <class T>void testSingle(const T& I, const std::string& tag,const ext::function<Real (Real)>& f,Real xMin, Real xMax, Real expected) {Real calculated = I(f,xMin,xMax);if (std::fabs(calculated-expected) > integrals_test::tolerance) {BOOST_FAIL(std::setprecision(10)<< "integrating " << tag<< "    calculated: " << calculated<< "    expected:   " << expected);}}template <class T>void testSeveral(const T& I) {testSingle(I, "f(x) = 0", [](Real x) -> Real { return 0.0; }, 0.0, 1.0, 0.0);testSingle(I, "f(x) = 1", [](Real x) -> Real { return 1.0; }, 0.0, 1.0, 1.0);testSingle(I, "f(x) = x", [](Real x) -> Real { return x; }, 0.0, 1.0, 0.5);testSingle(I, "f(x) = x^2",[](Real x) -> Real { return x * x; }, 0.0, 1.0, 1.0/3.0);testSingle(I, "f(x) = sin(x)",[](Real x) -> Real { return std::sin(x); }, 0.0, M_PI, 2.0);testSingle(I, "f(x) = cos(x)",[](Real x) -> Real { return std::cos(x); }, 0.0, M_PI, 0.0);testSingle(I, "f(x) = Gaussian(x)",NormalDistribution(), -10.0, 10.0, 1.0);testSingle(I, "f(x) = Abcd2(x)",AbcdSquared(0.07, 0.07, 0.5, 0.1, 8.0, 10.0), 5.0, 6.0,AbcdFunction(0.07, 0.07, 0.5, 0.1).covariance(5.0, 6.0, 8.0, 10.0));}template <class T>void testDegeneratedDomain(const T& I) {testSingle(I, "f(x) = 0 over [1, 1 + macheps]", [](Real x) -> Real { return 0.0; }, 1.0,1.0 + QL_EPSILON, 0.0);}}void IntegralTest::testSegment() {BOOST_TEST_MESSAGE("Testing segment integration...");using namespace integrals_test;testSeveral(SegmentIntegral(10000));testDegeneratedDomain(SegmentIntegral(10000));
}void IntegralTest::testTrapezoid() {BOOST_TEST_MESSAGE("Testing trapezoid integration...");using namespace integrals_test;testSeveral(TrapezoidIntegral<Default>(integrals_test::tolerance, 10000));testDegeneratedDomain(TrapezoidIntegral<Default>(integrals_test::tolerance, 10000));
}void IntegralTest::testMidPointTrapezoid() {BOOST_TEST_MESSAGE("Testing mid-point trapezoid integration...");using namespace integrals_test;testSeveral(TrapezoidIntegral<MidPoint>(integrals_test::tolerance, 10000));testDegeneratedDomain(TrapezoidIntegral<MidPoint>(integrals_test::tolerance, 10000));
}void IntegralTest::testSimpson() {BOOST_TEST_MESSAGE("Testing Simpson integration...");using namespace integrals_test;testSeveral(SimpsonIntegral(integrals_test::tolerance, 10000));testDegeneratedDomain(SimpsonIntegral(integrals_test::tolerance, 10000));
}void IntegralTest::testGaussKronrodAdaptive() {BOOST_TEST_MESSAGE("Testing adaptive Gauss-Kronrod integration...");using namespace integrals_test;Size maxEvaluations = 1000;testSeveral(GaussKronrodAdaptive(integrals_test::tolerance, maxEvaluations));testDegeneratedDomain(GaussKronrodAdaptive(integrals_test::tolerance, maxEvaluations));
}void IntegralTest::testGaussLobatto() {BOOST_TEST_MESSAGE("Testing adaptive Gauss-Lobatto integration...");using namespace integrals_test;Size maxEvaluations = 1000;testSeveral(GaussLobattoIntegral(maxEvaluations, integrals_test::tolerance));// on degenerated domain [1,1+macheps] an exception is thrown// which is also ok, but not tested here
}void IntegralTest::testTanhSinh() {BOOST_TEST_MESSAGE("Testing tanh-sinh integration...");using namespace integrals_test;testSeveral(TanhSinhIntegral());
}void IntegralTest::testGaussLegendreIntegrator() {BOOST_TEST_MESSAGE("Testing Gauss-Legendre integrator...");using namespace integrals_test;const GaussLegendreIntegrator integrator(64);testSeveral(integrator);testDegeneratedDomain(integrator);
}void IntegralTest::testGaussChebyshevIntegrator() {BOOST_TEST_MESSAGE("Testing Gauss-Chebyshev integrator...");using namespace integrals_test;const GaussChebyshevIntegrator integrator(64);testSingle(integrator, "f(x) = Gaussian(x)",NormalDistribution(), -10.0, 10.0, 1.0);testDegeneratedDomain(integrator);
}void IntegralTest::testGaussChebyshev2ndIntegrator() {BOOST_TEST_MESSAGE("Testing Gauss-Chebyshev 2nd integrator...");using namespace integrals_test;const GaussChebyshev2ndIntegrator integrator(64);testSingle(integrator, "f(x) = Gaussian(x)",NormalDistribution(), -10.0, 10.0, 1.0);testDegeneratedDomain(integrator);
}void IntegralTest::testGaussKronrodNonAdaptive() {BOOST_TEST_MESSAGE("Testing non-adaptive Gauss-Kronrod integration...");using namespace integrals_test;Real precision = integrals_test::tolerance;Size maxEvaluations = 100;Real relativeAccuracy = integrals_test::tolerance;GaussKronrodNonAdaptive gaussKronrodNonAdaptive(precision, maxEvaluations,relativeAccuracy);testSeveral(gaussKronrodNonAdaptive);testDegeneratedDomain(gaussKronrodNonAdaptive);
}void IntegralTest::testTwoDimensionalIntegration() {BOOST_TEST_MESSAGE("Testing two dimensional adaptive ""Gauss-Lobatto integration...");using namespace integrals_test;const Size maxEvaluations = 1000;const Real calculated = TwoDimensionalIntegral(ext::shared_ptr<Integrator>(new TrapezoidIntegral<Default>(integrals_test::tolerance, maxEvaluations)),ext::shared_ptr<Integrator>(new TrapezoidIntegral<Default>(integrals_test::tolerance, maxEvaluations)))(std::multiplies<>(),std::make_pair(0.0, 0.0), std::make_pair(1.0, 2.0));const Real expected = 1.0;if (std::fabs(calculated-expected) > integrals_test::tolerance) {BOOST_FAIL(std::setprecision(10)<< "two dimensional integration: "<< "\n    calculated: " << calculated<< "\n    expected:   " << expected);}
}namespace integrals_test {class sineF {public:Real operator()(Real x) const {return std::exp(-0.5*(x - M_PI_2/100));}};class cosineF {public:Real operator()(Real x) const {return std::exp(-0.5*x);}};}void IntegralTest::testFolinIntegration() {BOOST_TEST_MESSAGE("Testing Folin's integral formulae...");using namespace integrals_test;// Examples taken from// http://www.tat.physik.uni-tuebingen.de/~kokkotas/Teaching/Num_Methods_files/Comp_Phys5.pdfconst Size nr[] = { 4, 8, 16, 128, 256, 1024, 2048 };const Real expected[] = { 4.55229440e-5,4.72338540e-5, 4.72338540e-5,4.78308678e-5,4.78404787e-5, 4.78381120e-5,4.78381084e-5};const Real t = 100;const Real o = M_PI_2/t;const Real tol = 1e-12;for (Size i=0; i < LENGTH(nr); ++i) {const Size n = nr[i];const Real calculatedCosine= FilonIntegral(FilonIntegral::Cosine, t, n)(cosineF(),0,2*M_PI);const Real calculatedSine= FilonIntegral(FilonIntegral::Sine, t, n)(sineF(), o,2*M_PI + o);if (std::fabs(calculatedCosine-expected[i]) > tol) {BOOST_FAIL(std::setprecision(10)<< "Filon Cosine integration failed: "<< "\n    calculated: " << calculatedCosine<< "\n    expected:   " << expected[i]);}if (std::fabs(calculatedSine-expected[i]) > tol) {BOOST_FAIL(std::setprecision(10)<< "Filon Sine integration failed: "<< "\n    calculated: " << calculatedCosine<< "\n    expected:   " << expected[i]);}}
}namespace integrals_test {Real f1(Real x) {return 1.2*x*x+3.2*x+3.1;}Real f2(Real x) {return 4.3*(x-2.34)*(x-2.34)-6.2*(x-2.34) + f1(2.34);}}void IntegralTest::testDiscreteIntegrals() {BOOST_TEST_MESSAGE("Testing discrete integral formulae...");using namespace integrals_test;Array x(6), f(6);x[0] = 1.0; x[1] = 2.02; x[2] = 2.34; x[3] = 3.3; x[4] = 4.2; x[5] = 4.6;std::transform(x.begin(), x.begin()+3, f.begin(),   f1);std::transform(x.begin()+3, x.end(),   f.begin()+3, f2);const Real expectedSimpson =16.0401216 + 30.4137528 + 0.2*f2(4.2) + 0.2*f2(4.6);const Real expectedTrapezoid =0.5*(f1(1.0)  + f1(2.02))*1.02+ 0.5*(f1(2.02) + f1(2.34))*0.32+ 0.5*(f2(2.34) + f2(3.3) )*0.96+ 0.5*(f2(3.3)  + f2(4.2) )*0.9+ 0.5*(f2(4.2)  + f2(4.6) )*0.4;const Real calculatedSimpson =  DiscreteSimpsonIntegral()(x, f);const Real calculatedTrapezoid = DiscreteTrapezoidIntegral()(x, f);const Real tol = 1e-12;if (std::fabs(calculatedSimpson-expectedSimpson) > tol) {BOOST_FAIL(std::setprecision(16)<< "discrete Simpson integration failed: "<< "\n    calculated: " << calculatedSimpson<< "\n    expected:   " << expectedSimpson);}if (std::fabs(calculatedTrapezoid-expectedTrapezoid) > tol) {BOOST_FAIL(std::setprecision(16)<< "discrete Trapezoid integration failed: "<< "\n    calculated: " << calculatedTrapezoid<< "\n    expected:   " << expectedTrapezoid);}
}void IntegralTest::testDiscreteIntegrator() {BOOST_TEST_MESSAGE("Testing discrete integrator formulae...");using namespace integrals_test;testSeveral(DiscreteSimpsonIntegrator(300));testSeveral(DiscreteTrapezoidIntegrator(3000));
}namespace integrals_test{std::vector<Real> x, y;Real pw_fct(const Real t) { return QL_PIECEWISE_FUNCTION(x, y, t); }void pw_check(const Integrator &in, const Real a, const Real b,const Real expected) {Real calculated = in(pw_fct, a, b);if (!close(calculated, expected))BOOST_FAIL(std::setprecision(16)<< "piecewise integration over [" << a << "," << b<< "] failed: "<< "\n   calculated: " << calculated<< "\n   expected:   " << expected<< "\n   difference: " << (calculated - expected));
}
} // empty namespacevoid IntegralTest::testPiecewiseIntegral() {BOOST_TEST_MESSAGE("Testing piecewise integral...");using namespace integrals_test;x = { 1.0, 2.0, 3.0, 4.0, 5.0 };y = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };ext::shared_ptr<Integrator> segment =ext::make_shared<SegmentIntegral>(1);ext::shared_ptr<Integrator> piecewise =ext::make_shared<PiecewiseIntegral>(segment, x);pw_check(*piecewise, -1.0, 0.0, 1.0);pw_check(*piecewise, 0.0, 1.0, 1.0);pw_check(*piecewise, 0.0, 1.5, 2.0);pw_check(*piecewise, 0.0, 2.0, 3.0);pw_check(*piecewise, 0.0, 2.5, 4.5);pw_check(*piecewise, 0.0, 3.0, 6.0);pw_check(*piecewise, 0.0, 4.0, 10.0);pw_check(*piecewise, 0.0, 5.0, 15.0);pw_check(*piecewise, 0.0, 6.0, 21.0);pw_check(*piecewise, 0.0, 7.0, 27.0);pw_check(*piecewise, 3.5, 4.5, 4.5);pw_check(*piecewise, 5.0, 10.0, 30.0);pw_check(*piecewise, 9.0, 10.0, 6.0);
}namespace integrals_test {template <class T>void reportSiCiFail(const std::string& name, T z, T c, T e, Real diff, Real tol) {BOOST_FAIL(std::setprecision(16)<< name << " calculation failed for " << z<< "\n calculated: " << c<< "\n expected:   " << e<< "\n difference: " << diff<< "\n tolerance:  " << tol);}
}void IntegralTest::testExponentialIntegral() {BOOST_TEST_MESSAGE("Testing exponential integrals...");using namespace ExponentialIntegral;// reference values are calculated with Mathematica or Python/mpmathconst Real data[][10] = {{0.0010000000000000000208, 0.0, 0.00099999994444444613193, 0.0, -6.330539864080593754, 0.0, -6.3295393640250381967, 0.0, 6.3315393641361493112, 0.0},{0.00070710678118654756077, 0.00070710678118654751747, 0.00070710682047025644818, 0.00070710674190283627304, -6.3305396140806145873, 0.785397913397448279, -6.329832507338711751, 0.7861055202179185354, 6.3312467208225174236, -0.78469130657697802259},{6.1232339957367660136e-20, 0.0010000000000000000208, 6.1232350162201617204e-20, 0.001000000055555557243, -6.330539364080593754, 1.570796326794896558, -6.3305398640805937539, 1.5717963267393410041, 6.330539864080593754, -1.5697963268504521119},{-0.00070710678118654747417, 0.00070710678118654760407, -0.00070710682047025636158, 0.00070710674190283635964, -6.3305396140806145873, 2.356194740192344837, -6.3312467208225174235, 2.3569013470128150935, 6.3298325073387117511, -2.3554871333718745805},{-0.0010000000000000000208, 1.2246467991473532027e-19, -0.00099999994444444613193, 1.2246465950321797194e-19, -6.330539864080593754, 3.141592653589793116, -6.3315393641361493112, 3.1415926535897931161, 6.3295393640250381967, -3.1415926535897931159},{-0.00070710678118654764736, -0.00070710678118654743088, -0.00070710682047025653477, -0.00070710674190283618645, -6.3305396140806145873, -2.3561947401923450819, -6.3312467208225174237, -2.3569013470128153382, 6.3298325073387117509, 2.3554871333718748256},{-1.8369701987210298041e-19, -0.0010000000000000000208, -1.8369705049015472569e-19, -0.001000000055555557243, -6.330539364080593754, -1.5707963267948968029, -6.3305398640805937541, -1.5717963267393412491, 6.3305398640805937538, 1.5697963268504523568},{0.00070710678118654738758, -0.00070710678118654769066, 0.00070710682047025627499, -0.00070710674190283644623, -6.3305396140806145873, -0.78539791339744852393, -6.3298325073387117512, -0.78610552021791878051, 6.3312467208225174234, 0.78469130657697826735},{0.10000000000000000555, 0.0, 0.099944461108276955702, 0.0, -1.7278683866572965838, 0.0, -1.6228128139692766136, 0.0, 1.8229239584193906159, 0.0},{0.07071067811865475853, 0.0707106781186547542, 0.070749950041603603669, 0.070671382625480302524, -1.7253704697591484327, 0.78289816362892975745, -1.6546990871336681256, 0.8586481132075703999, 1.7960418523846287393, -0.7171482131243632012},{6.123233995736766226e-18, 0.10000000000000000555, 6.1334444895968712791e-18, 0.10005557222505700112, -1.7228683861943336147, 1.5707963267948965577, -1.7278683866572965777, 1.670740787903173514, 1.7278683866572965899, -1.4708518656866196026},{-0.070710678118654749871, 0.07071067811865476286, -0.070749950041603595024, 0.070671382625480311198, -1.7253704697591484321, 2.3586944899608633586, -1.7960418523846287312, 2.4244444404654299234, 1.6546990871336681349, -2.2829445403822227075},{-0.10000000000000000555, 1.2246467991473532452e-17, -0.099944461108276955702, 1.2226067414369268591e-17, -1.7278683866572965838, 3.1415926535897931166, -1.8229239584193906159, 3.1415926535897931277, 1.6228128139692766136, -3.1415926535897931031},{-0.07071067811865476719, -0.070710678118654745541, -0.070749950041603612314, -0.07067138262548029385, -1.7253704697591484333, -2.3586944899608636035, -1.7960418523846287474, -2.4244444404654301511, 1.6546990871336681163, 2.2829445403822229697},{-1.8369701987210298678e-17, -0.10000000000000000555, -1.8400333469093536425e-17, -0.10005557222505700112, -1.7228683861943336147, -1.5707963267948968038, -1.7278683866572966021, -1.6707407879031737577, 1.7278683866572965654, 1.4708518656866198463},{0.070710678118654741211, -0.07071067811865477152, 0.070749950041603586379, -0.070671382625480319872, -1.7253704697591484315, -0.78289816362893000238, -1.6546990871336681441, -0.85864811320757066212, 1.7960418523846287232, 0.71714821312436342884},{0.25, 0.0, 0.2491335703197571641, 0.0, -0.82466306258094565309, 0.0, -0.54254326466191372953, 0.0, 1.0442826344437381945, 0.0},{0.17677669529663688651, 0.17677669529663687569, 0.17738935115398995617, 0.17616173766104897661, -0.80911938627521916259, 0.76977321991145556311, -0.63295764861417017313, 0.97841245803743094036, 0.98528112393626814822, -0.62363375572945104943},{1.5308084989341914715e-17, 0.25, 1.5468043259937507815e-17, 0.2508696848909122325, -0.79341294955282596203, 1.5707963267948965561, -0.82466306258094563794, 1.819929897114653724, 0.82466306258094566823, -1.3216627564751393958},{-0.17677669529663686486, 0.17677669529663689734, -0.17738935115398993475, 0.17616173766104899848, -0.80911938627521915876, 2.3718194336783375529, -0.98528112393626813018, 2.517958897860342088, 0.63295764861417019883, -2.1631801955523621542},{-0.25, 3.0616169978683829431e-17, -0.2491335703197571641, 3.0298246679137807606e-17, -0.82466306258094565309, 3.1415926535897931198, -1.0442826344437381945, 3.1415926535897931431, 0.54254326466191372953, -3.1415926535897930812},{-0.17677669529663690816, -0.17677669529663685404, -0.1773893511539899776, -0.17616173766104895473, -0.80911938627521916642, -2.3718194336783377978, -0.98528112393626816627, -2.5179588978603422901, 0.63295764861417014743, 2.163180195552362442},{-4.5924254968025744146e-17, -0.25, -4.6404129781529084775e-17, -0.2508696848909122325, -0.79341294955282596203, -1.5707963267948968087, -0.82466306258094569853, -1.8199298971146539613, 0.82466306258094560764, 1.3216627564751396331},{0.17677669529663684321, -0.17677669529663691899, 0.17738935115398991333, -0.17616173766104902036, -0.80911938627521915494, -0.769773219911455808, -0.63295764861417022453, -0.97841245803743122809, 0.98528112393626811213, 0.62363375572945125148},{0.5, 0.0, 0.49310741804306668916, 0.0, -0.17778407880661290134, 0.0, 0.45421990486317357992, 0.0, 0.55977359477616081175, 0.0},{0.35355339059327377302, 0.35355339059327375138, 0.3584268697133189464, 0.34860625536259795411, -0.11658254521497154353, 0.72290178026868503268, 0.23202371014762644078, 1.2063214162395304511, 0.46518880057756951253, -0.48946767681289259981},{3.0616169978683829431e-17, 0.5, 3.190788489546069124e-17, 0.50699674981966719583, -0.052776844956493615913, 1.5707963267948965502, -0.17778407880661287198, 2.0639037448379632547, 0.17778407880661293069, -1.0776889087518298763},{-0.35355339059327372973, 0.35355339059327379467, -0.35842686971331890493, 0.34860625536259799919, -0.11658254521497152822, 2.4186908733211080836, -0.46518880057756948275, 2.652124976776900558, -0.2320237101476263804, -1.9352712373502626237},{-0.5, 6.1232339957367658861e-17, -0.49310741804306668916, 5.8712695126944314728e-17, -0.17778407880661290134, 3.141592653589793131, -0.55977359477616081175, 3.1415926535897931642, -0.45421990486317357992, -3.1415926535897930366},{-0.35355339059327381632, -0.35355339059327370808, -0.35842686971331898787, -0.34860625536259790903, -0.11658254521497155883, -2.4186908733211083279, -0.4651888005775695423, -2.6521249767769007193, -0.23202371014762650116, 1.9352712373502629509},{-9.1848509936051488292e-17, -0.5, -9.5723654690421041556e-17, -0.50699674981966719583, -0.052776844956493615913, -1.5707963267948968264, -0.1777840788066129894, -2.0639037448379634696, 0.17778407880661281327, 1.0776889087518300913},{0.35355339059327368643, -0.35355339059327383797, 0.35842686971331886346, -0.34860625536259804427, -0.11658254521497151291, -0.72290178026868527697, 0.23202371014762632001, -1.2063214162395307784, 0.46518880057756945298, 0.48946767681289276116},{1.0, 0.0, 0.94608307036718301494, 0.0, 0.33740392290096813466, 0.0, 1.8951178163559367555, 0.0, 0.21938393439552027368, 0.0},{0.70710678118654754605, 0.70710678118654750275, 0.74519215535365930209, 0.66666481741950631398, 0.56680209825930888934, 0.53562961732242986806, 1.233466915678815284, 1.7803588648261259588, 0.099862719160197444251, -0.28997455411880742612},{6.1232339957367658861e-17, 1.0, 7.1960319005373041642e-17, 1.0572508753757285146, 0.83786694098020824089, 1.5707963267948965247, 0.33740392290096818619, 2.5168793971620796011, -0.33740392290096808314, -0.62471325642771357121},{-0.70710678118654745945, 0.70710678118654758935, -0.74519215535365923063, 0.66666481741950641427, 0.5668020982593089504, 2.605963036267363253, -0.099862719160197405024, 2.8516180994709857664, -1.2334669156788151226, -1.3612337887636670908},{-1.0, 1.2246467991473531772e-16, -0.94608307036718301494, 1.0305047480945999774e-16, 0.33740392290096813466, 3.1415926535897931723, -0.21938393439552027368, 3.1415926535897931934, -1.8951178163559367555, -3.1415926535897929056},{-0.70710678118654763265, -0.70710678118654741616, -0.74519215535365937355, -0.66666481741950621369, 0.56680209825930882828, -2.6059630362673634878, -0.099862719160197483478, -2.8516180994709858582, -1.2334669156788154453, 1.3612337887636674684},{-1.8369701987210297658e-16, -1.0, -2.158809570266204413e-16, -1.0572508753757285146, 0.83786694098020824089, -1.5707963267948969027, 0.33740392290096798009, -2.5168793971620797334, -0.33740392290096828924, 0.62471325642771370354},{0.70710678118654737286, -0.70710678118654767594, 0.74519215535365915917, -0.66666481741950651456, 0.56680209825930901147, -0.53562961732243010279, 1.2334669156788149613, -1.7803588648261263365, 0.099862719160197365796, 0.28997455411880751794},{1.5, 0.0, 1.3246835311721196804, 0.0, 0.47035631719539988668, 0.0, 3.301285449129797838, 0.0, 0.1000195824066326519, 0.0},{1.0606601717798213191, 1.0606601717798212541, 1.1839593855572855609, 0.9194789610047786165, 0.93002583013377829907, 0.22553329329273497242, 1.8495047911385570699, 2.5292224190594471214, -0.010546869128999664075, -0.16130364794487607552},{9.1848509936051488292e-17, 1.5, 1.3038076345658281264e-16, 1.7006525157682152449, 1.600632933361582593, 1.5707963267948964752, 0.47035631719539994775, 2.8954798579670162953, -0.4703563171953998256, -0.24611279562277693453},{-1.0606601717798211892, 1.060660171779821384, -1.1839593855572854849, 0.91947896100477878934, 0.93002583013377843491, 2.9160593602970581693, 0.010546869128999701077, 2.9802890056449171422, -1.8495047911385567612, -0.61237023453034594437},{-1.5, 1.8369701987210297658e-16, -1.3246835311721196804, 1.2215790424776667532e-16, 0.47035631719539988668, 3.1415926535897932298, -0.1000195824066326519, 3.1415926535897932111, -3.301285449129797838, -3.1415926535897926896},{-1.060660171779821449, -1.0606601717798211242, -1.1839593855572856369, -0.91947896100477844366, 0.93002583013377816324, -2.9160593602970583627, 0.010546869128999627072, -2.9802890056449171836, -1.8495047911385573786, 0.6123702345303462898},{-2.7554552980815446488e-16, -1.5, -3.9114229038065365107e-16, -1.7006525157682152449, 1.600632933361582593, -1.5707963267948970514, 0.47035631719539970344, -2.8954798579670163126, -0.47035631719540006991, 0.24611279562277695186},{1.0606601717798210593, -1.0606601717798215139, 1.1839593855572854089, -0.91947896100477896218, 0.93002583013377857075, -0.22553329329273516584, 1.8495047911385564526, -2.5292224190594474668, -0.010546869128999738079, 0.16130364794487611693},{2.0, 0.0, 1.6054129768026948486, 0.0, 0.4229808287748649957, 0.0, 4.9542343560018901634, 0.0, 0.048900510708061119567, 0.0},{1.4142135623730950921, 1.4142135623730950055, 1.6883194933582990243, 1.0649044719839209475, 1.1044891171909028954, -0.19981522706084708457, 2.1693935891748240916, 3.4589310472140426889, -0.039584645206981933184, -0.08229206049744467714},{1.2246467991473531772e-16, 2.0, 2.2208114946641807464e-16, 2.5015674333549756415, 2.4526669226469145219, 1.5707963267948963889, 0.42298082877486505138, 3.1762093035975914933, -0.42298082877486494002, 0.034616650007798203864},{-1.4142135623730949189, 1.4142135623730951787, -1.6883194933582989874, 1.0649044719839212109, 1.1044891171909031294, 3.3414078806506402814, 0.039584645206981962593, 3.0593005930923485567, -2.169393589174823594, 0.31733839362424952895},{-2.0, 2.4492935982947063545e-16, -1.6054129768026948486, 1.1135681831696612616e-16, 0.4229808287748649957, 3.1415926535897932894, -0.048900510708061119567, 3.1415926535897932219, -4.9542343560018901634, -3.1415926535897923336},{-1.4142135623730952653, -1.4142135623730948323, -1.6883194933582990613, -1.064904471983920684, 1.1044891171909026613, -3.3414078806506403646, 0.039584645206981903775, -3.059300593092348566, -2.1693935891748245892, -0.31733839362424937184},{-3.6739403974420595317e-16, -2.0, -6.6624344842268023737e-16, -2.5015674333549756415, 2.4526669226469145219, -1.5707963267948973103, 0.42298082877486482866, -3.1762093035975913914, -0.42298082877486516273, -0.03461665000779830579},{1.4142135623730947457, -1.4142135623730953519, 1.6883194933582989504, -1.0649044719839214744, 1.1044891171909033635, 0.19981522706084700137, 2.1693935891748230965, -3.458931047214042846, -0.039584645206981992002, 0.082292060497444686426},{5.0, 0.0, 1.5499312449446741373, 0.0, -0.19002974965664387862, 0.0, 40.185275355803177455, 0.0, 0.0011482955912753257973, 0.0},{3.5355339059327377302, 3.5355339059327375138, 3.6871508611543429709, -3.157181373909001357, -3.1547681046716106125, -2.1118502909279661892, -6.311949478580612776, 7.3697974788772077195, -0.0024132692373907452624, 0.0045042434314801607547},{3.0616169978683829431e-16, 5.0, 4.5436362172750887551e-15, 20.09321182569722639, 20.092063530105951065, 1.5707963267948920752, -0.19002974965664393734, 3.1207275717395707391, 0.1900297496566438199, -0.020865081850222464588},{-3.5355339059327372973, 3.5355339059327379467, -3.6871508611543449094, -3.1571813739090021642, -3.1547681046716114182, 5.2534429445177613695, 0.0024132692373907438925, 3.1460968970212734025, 6.3119494785806111631, 4.2282048252874106008},{-5.0, 6.1232339957367658861e-16, -1.5499312449446741373, -1.174343542809398959e-16, -0.19002974965664387862, 3.1415926535897932037, -0.0011482955912753257973, 3.1415926535897932376, -40.185275355803177455, -3.1415926535897750631},{-3.5355339059327381632, -3.5355339059327370808, -3.6871508611543410324, 3.1571813739090005499, -3.1547681046716098067, -5.2534429445177574859, 0.0024132692373907466323, -3.1460968970212733959, 6.3119494785806143889, -4.2282048252874183613},{-9.1848509936051488292e-16, -5.0, -1.3630908648516543815e-14, -20.09321182569722639, 20.092063530105951065, -1.5707963267949102514, -0.19002974965664370247, -3.1207275717395708086, 0.19002974965664405477, 0.020865081850222534065},{3.5355339059327368643, -3.5355339059327383797, 3.6871508611543468479, 3.1571813739090029713, -3.154768104671612224, 2.1118502909279700728, -6.3119494785806095501, -7.3697974788771999589, -0.0024132692373907425226, -0.004504243431480167346},{10.0, 0.0, 1.6583475942188740493, 0.0, -0.045456433004455372635, 0.0, 2492.2289762418777591, 0.0, 4.1569689296853242774e-6, 0.0},{7.0710678118654754605, 7.0710678118654750275, -3.7745175303432929615, 62.642575559232345941, 62.642571122904060317, 5.3452347019798827768, 125.28514668213645736, -7.5489559055283299711, 4.4363282856236060522e-6, -0.000079155158306803868409},{6.1232339957367658861e-16, 10.0, 6.7436601941373126554e-13, 1246.1144901994233444, 1246.1144860424544147, 1.5707963267942222532, -0.045456433004455405946, 3.2291439210137707199, 0.045456433004455339323, 0.087551267423977378721},{-7.0710678118654745945, 7.0710678118654758935, 3.7745175303433438135, 62.642575559232397046, 62.642571122904111422, -2.2036420483901403904, -4.4363282856235323221e-6, 3.1415134984314864345, -125.28514668213635515, -10.690548559118021506},{-10.0, 1.2246467991473531772e-15, -1.6583475942188740493, -6.6623371218859982286e-17, -0.045456433004455372635, 3.1415926535897933412, -4.1569689296853242774e-6, 3.1415926535897932385, -2492.2289762418777591, -3.1415926535870957744},{-7.0710678118654763265, -7.0710678118654741616, 3.7745175303432421091, -62.642575559232294835, 62.642571122904009212, 2.2036420483900386859, -4.4363282856236797822e-6, -3.1415134984314864347, -125.28514668213655957, 10.690548559118224914},{-1.8369701987210297658e-15, -10.0, -2.0230980582411937966e-12, -1246.1144901994233444, 1246.1144860424544147, -1.5707963267969197173, -0.045456433004455272699, -3.2291439210137705144, 0.04545643300445547257, -0.087551267423977584235},{7.0710678118654737286, -7.0710678118654767594, -3.7745175303433946658, -62.642575559232448152, 62.642571122904162528, -5.3452347019799844813, 125.28514668213625294, 7.5489559055281265623, 4.4363282856234585915e-6, 0.000079155158306804015139},{15.0, 0.0, 1.6181944437083687391, 0.0, 0.046278677674360439604, 0.0, 234955.85249076830358, 0.0, 1.9186278921478669771e-8, 0.0},{10.606601717798213191, 10.606601717798212541, -471.05515346873571713, -1328.2594907541181068, -1328.2594913017652105, 472.62595127241629118, -2656.5189820558856064, -942.11030841435617353, 5.4764710367805069957e-7, 1.4768856774331466472e-6},{9.1848509936051488292e-16, 15.0, 1.0008479153886857898e-10, 117477.92624539374493, 117477.92624537455865, 1.5707963266948118277, 0.046278677674360479423, 3.1889907705032654049, -0.046278677674360399786, 0.047398116913472073375},{-10.606601717798211892, 10.60660171779821384, 471.05515346873477896, -1328.2594907541203959, -1328.2594913017674995, -469.48435861882555978, -5.4764710367805350438e-7, 3.1415941304754706716, 2656.5189820558810283, -945.2519010679478431},{-15.0, 1.8369701987210297658e-15, -1.6181944437083687391, 7.9637292204322254463e-17, 0.046278677674360439604, 3.1415926535897933315, -1.9186278921478669771e-8, 3.1415926535897932385, -234955.85249076830358, -3.1415926531894540723},{-10.60660171779821449, -10.606601717798211242, 471.05515346873665531, 1328.2594907541158178, -1328.2594913017629215, 469.48435861882743613, -5.4764710367804789476e-7, -3.1415941304754706716, 2656.5189820558901845, 945.25190106794409045},{-2.7554552980815446488e-15, -15.0, -3.0025437461660077385e-10, -117477.92624539374493, 117477.92624537455865, -1.5707963270951509938, 0.046278677674360320148, -3.1889907705032652188, -0.04627867767436055906, -0.047398116913472259445},{10.606601717798210593, -10.606601717798215139, -471.05515346873384079, 1328.2594907541226849, -1328.2594913017697886, -472.62595127241441484, -2656.5189820558764503, 942.11030841435992621, 5.4764710367805630921e-7, -1.4768856774331489463e-6},{20.0, 0.0, 1.5482417010434398402, 0.0, 0.04441982084535331654, 0.0, 25615652.66405658882, 0.0, 9.8355252906498816904e-11, 0.0},{14.142135623730950921, 14.142135623730950055, 24486.68965358318626, 26235.092183360235647, 26235.092183384186671, -24485.118857281672279, 52470.184366744507202, 48973.379307191653857, -2.3951024244493564362e-8, -2.5280915398646225218e-8},{1.2246467991473531772e-15, 20.0, 1.4853900090407493933e-8, 12807826.332028294459, 12807826.332028294361, 1.5707963119409965288, 0.044419820845353372442, 3.1190380278383364344, -0.044419820845353260638, -0.02255462575145675408},{-14.142135623730949189, 14.142135623730951787, -24486.689653583186682, 26235.092183360320531, 26235.092183384271555, 24488.260449935262494, 2.3951024244493652701e-8, 3.1415926283088778398, -52470.184366744337435, 48970.23771453806322},{-20.0, 2.4492935982947063545e-15, -1.5482417010434398402, 1.1180354790034662799e-16, 0.04441982084535331654, 3.1415926535897931885, -9.8355252906498816904e-11, 3.1415926535897932385, -25615652.66405658882, -3.1415925941741928768},{-14.142135623730952653, -14.142135623730948323, -24486.689653583185838, -26235.092183360150763, 26235.092183384101787, -24488.26044993526165, 2.3951024244493476023e-8, -3.1415926283088778398, -52470.184366744676971, -48970.237714538064908},{-3.6739403974420595317e-15, -20.0, -4.4561700271222483452e-8, -12807826.332028294459, 12807826.332028294361, -1.5707963713565968905, 0.044419820845353148834, -3.1190380278383365344, -0.044419820845353484245, 0.022554625751456854031},{14.142135623730947457, -14.142135623730953519, 24486.689653583187103, -26235.092183360405415, 26235.09218338435644, 24485.118857281673122, 52470.184366744167666, -48973.37930719165217, -2.3951024244493741041e-8, 2.528091539864622434e-8},{24.989999999999998437, 0.0, 1.5315374843262580141, 0.0, -0.0072448862399538447391, 0.0, 2977286754.0962403117, 0.0, 5.4047414055481939457e-13, 0.0},{17.67059846185182207, 17.670598461851820988, -885673.16331345207799, -400248.00540652140574, -400248.00540652215797, 885674.7341097792085, -800496.01081304623667, -1771346.3266269055961, 7.5223231588352466724e-10, 3.3560562518531656391e-10},{1.5301961755346176992e-15, 24.989999999999998437, 2.1825789542468514196e-6, 1488643377.0481201559, 1488643377.0481201559, 1.5707941442159423724, -0.0072448862399538534498, 3.1023338111211545727, 0.0072448862399538360284, -0.039258842468638544566},{-17.670598461851819906, 17.670598461851823152, 885673.16331345318248, -400248.0054065240787, -400248.00540652483093, -885671.59251712672318, -7.5223231588352706351e-10, 3.1415926539253988636, 800496.01081304089075, -1771349.4682195569769},{-24.989999999999998437, 3.0603923510692353985e-15, -1.5315374843262580141, -1.7421457417638512695e-17, -0.0072448862399538447391, 3.1415926535897931172, -5.4047414055481939457e-13, 3.1415926535897932385, -2977286754.0962403117, -3.1415839232739762511},{-17.670598461851824234, -17.670598461851818824, 885673.16331345097349, 400248.00540651873277, -400248.005406519485, 885671.5925171245142, -7.5223231588352227096e-10, -3.1415926539253988636, 800496.01081305158263, 1771349.4682195613948},{-4.5905885266038530977e-15, -24.989999999999998437, -6.5477368627405542588e-6, -1488643377.0481201559, 1488643377.0481201559, -1.5708028745317593598, -0.0072448862399538186069, -3.1023338111211548151, 0.0072448862399538708713, 0.039258842468638787005},{17.670598461851817742, -17.670598461851825316, -885673.16331345428695, 400248.00540652675168, -400248.00540652750391, -885674.73410978141745, -800496.01081303554482, 1771346.3266269011781, 7.5223231588352945979e-10, -3.3560562518531458362e-10}};const Real tol = 1e-10;for (const auto& i : data) {const Real x = i[0];const Real y = (std::abs(i[1]) < 1e-12) ? 0.0 : i[1];const std::complex<Real> z(x, y);const std::complex<Real> si = Si(z);std::complex<Real> ref(i[2], i[3]);Real diff = std::abs(si-ref)/std::abs(ref);if (diff > tol|| (std::abs(ref.real()) < tol && std::abs(si.real()) > tol)|| (std::abs(ref.imag()) < tol && std::abs(si.imag()) > tol)) {integrals_test::reportSiCiFail("Si", z, si, ref, diff, tol);}const std::complex<Real> ci = Ci(z);ref = std::complex<Real>(i[4], i[5]);diff = std::min(std::abs(ci-ref), std::abs(ci-ref)/std::abs(ref));if (diff > tol|| (std::abs(ref.real()) < tol && std::abs(ci.real()) > tol)|| (std::abs(ref.imag()) < tol && std::abs(ci.imag()) > tol)) {integrals_test::reportSiCiFail("Ci", z, ci, ref, diff, tol);}const std::complex<Real> ei = Ei(z);ref = std::complex<Real>(i[6], i[7]);diff = std::abs(ei-ref)/std::abs(ref);if (diff > tol|| (std::abs(ref.real()) < tol && std::abs(ei.real()) > tol)|| (std::abs(ref.imag()) < tol && std::abs(ei.imag()) > tol)) {integrals_test::reportSiCiFail("Ei", z, ei, ref, diff, tol);}const std::complex<Real> e1 = E1(z);ref = std::complex<Real>(i[8], i[9]);diff = std::min(std::abs(e1-ref), std::abs(e1-ref)/std::abs(ref));if (std::abs(z) < 10.0)if (diff > tol|| (std::abs(ref.real()) < tol && std::abs(e1.real()) > tol)|| (std::abs(ref.imag()) < tol && std::abs(e1.imag()) > tol)) {integrals_test::reportSiCiFail("E1", z, e1, ref, diff, tol);}}
}void IntegralTest::testRealSiCiIntegrals() {BOOST_TEST_MESSAGE("Testing real Ci and Si...");using namespace ExponentialIntegral;// reference values are calculated with Mathematica or Python/mpmathconst Real data[][3] = {{1e-12, 1e-12, -27.0538054510270153677},{0.1, 0.09994446110827695570, -1.7278683866572965838},{1.0, 0.9460830703671830149, 0.3374039229009681347},{1.9999, 1.6053675097543679041, 0.4230016343635392},{3.9999, 1.758222058430840841, -0.140965355646150101},{4.0001, 1.758184218306157867, -0.140998037827177150},{5.0, 1.5499312449446741373, -0.19002974965664387862},{7.0, 1.4545966142480935906, 0.076695278482184518383,},{10.0, 1.6583475942188740493, -0.045456433004455372635},{15.0, 1.6181944437083687391, 0.046278677674360439604},{20.0, 1.5482417010434398402, 0.04441982084535331654},{24.9, 1.532210740207620024, -0.010788215638781789846},{25.1, 1.5311526281483412938, -0.0028719014454227088097},{30.0, 1.566756540030351111, -0.033032417282071143779},{40.0, 1.5869851193547845068, 0.019020007896208766962},{400.0, 1.5721148692738117518, -0.00212398883084634893},{4000.0, 1.5709788562309441985, -0.00017083030544201591130}};const Real tol = 1e-12;for (const auto& i : data) {Real x = i[0];Real si = Si(x);Real diff = std::fabs(si - i[1]);if (diff > tol) {integrals_test::reportSiCiFail("SineIntegral", x, si, i[1], diff, tol);}const Real ci = Ci(x);diff = std::fabs(ci - i[2]);if (diff > tol) {integrals_test::reportSiCiFail("CosineIntegral", x, ci, i[2], diff, tol);}x = -i[0];si = Si(x);diff = std::fabs(si + i[1]);if (diff > tol) {integrals_test::reportSiCiFail("SineIntegral", x, si, Real(-i[1]), diff, tol);}}
}test_suite* IntegralTest::suite() {auto* suite = BOOST_TEST_SUITE("Integration tests");suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testSegment));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testTrapezoid));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testMidPointTrapezoid));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testSimpson));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussKronrodAdaptive));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussKronrodNonAdaptive));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussLobatto));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussLegendreIntegrator));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussChebyshevIntegrator));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussChebyshev2ndIntegrator));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testTwoDimensionalIntegration));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testFolinIntegration));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testDiscreteIntegrals));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testDiscreteIntegrator));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testPiecewiseIntegral));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testExponentialIntegral));suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testRealSiCiIntegrals));#ifdef QL_BOOST_HAS_TANH_SINHsuite->add(QUANTLIB_TEST_CASE(&IntegralTest::testTanhSinh));
#endifreturn suite;
}

该博文为原创文章，未经博主同意不得转。
本文章博客地址：https://cplusplus.blog.csdn.net/article/details/128364208

C++：实现量化Integration积分测试实例相关推荐

  1. C++：实现量化Libor市场模型测试实例

    C++:实现量化Libor市场模型测试实例 #include "libormarketmodel.hpp" #include "utilities.hpp"#i ...

  2. C++：实现量化ODE模型测试实例

    C++:实现量化ODE模型测试实例 #include "ode.hpp" #include "utilities.hpp" #include <ql/ex ...

  3. C++：实现量化默认概率曲线测试实例

    C++:实现量化默认概率曲线测试实例 #include "defaultprobabilitycurves.hpp" #include "utilities.hpp&qu ...

  4. C++：实现量化GSR模型测试实例

    C++:实现量化GSR模型测试实例 #include "gsr.hpp" #include "utilities.hpp" #include <ql/pr ...

  5. C++：实现量化covariance协方差矩阵测试实例

    C++:实现量化covariance协方差矩阵测试实例 #include "covariance.hpp" #include "utilities.hpp" # ...

  6. C++：实现量化MarketModels市场模型测试实例

    C++:实现量化MarketModels市场模型测试实例 #include <ql/qldefines.hpp> #if !defined(BOOST_ALL_NO_LIB) && ...

  7. C++：实现量化相关的Interpolation插值测试实例

    C++:实现量化相关的Interpolation插值测试实例 #include "interpolations.hpp" #include "utilities.hpp& ...

  8. C++：实现量化期权交易CDS加密货币衍生品测试实例

    C++:实现量化期权交易CDS加密货币衍生品测试实例 #include <ql/qldefines.hpp> #if !defined(BOOST_ALL_NO_LIB) &&am ...

  9. C++：实现量化daycounters 日计数器测试实例

    C++:实现量化daycounters 日计数器测试实例 #include "daycounters.hpp" #include "utilities.hpp" ...

最新文章

  1. 编译Android源码前的一个步骤
  2. 微信小程序实时将less编译为wxss
  3. linux文件目录的管理,Linux文件目录管理
  4. spring 集成mybatis——多数据源切换（附带定时器的配置）
  5. 信安教程第二版-第16章网络安全风险评估技术原理与应用
  6. c语言编程输出年月日,C语言程序设计： 输入年月日 然后输出是星期几
  7. (9)数据结构-双端队列
  8. 学习的四重境界，给上初中侄女，如何学习，如何定义社会人才
  9. ios带嗅探器的浏览器_MAC系统下 有没有像WIN系统的傲游浏览器那样有嗅探功能的呢?...
  10. 华泰证券 python 自动交易软件_为何选用股票自动交易助手
  11. 天地图显示不全的问题
  12. Labelling tools 的环境配置
  13. 深度Linux声卡驱动安装,Deepin Linux 的声卡驱动有点小问题
  14. Windows VScode SSH连接 Bad owner or permissions on C:\\Users\\admin/.ssh/config 错误解决方法
  15. 【Leetcode】1526. Minimum Number of Increments on Subarrays to Form a Target Array（配数学证明）
  16. python pip install fitter 失败解决方案
  17. 工控液晶屏开机白屏怎么回事,开机白屏解决方法？
  18. HTML网页图片使用技巧集锦
  19. kafka是什么？深刻理解kafka
  20. [jzoj 4787] 数格子 {矩阵乘法}

热门文章

  1. windows屏幕捕捉鼠标闪烁问题
  2. 一位大龄程序员所经历的面试的历炼和思考
  3. Es6 set和map
  4. 次元服务器的位置,暮色次元服务器介绍
  5. 用verilog HDL实现数字基带信号的2FSK调制
  6. 统计学习导论（ISLR）（三）：线性回归（超详细介绍）
  7. go语言初学者常见错误
  8. [Compose] 使用 Compose 在 Android 中的脚手架 Scaffold
  9. codevs1373 射命丸文
  10. Flutter 实战开发-网络请求