From 74c0de2f01054bcafe8765f01be7e58a9ef09fac Mon Sep 17 00:00:00 2001 From: pabristow Date: Wed, 3 Jul 2019 18:14:20 +0100 Subject: [PATCH] Three fourier integrals examples - work-in-progress --- ...ooura_fourier_integrals_cosine_example.cpp | 78 +++++++++++ example/ooura_fourier_integrals_example.cpp | 77 +++++++++++ ...urier_integrals_multiprecision_example.cpp | 125 ++++++++++++++++++ 3 files changed, 280 insertions(+) create mode 100644 example/ooura_fourier_integrals_cosine_example.cpp create mode 100644 example/ooura_fourier_integrals_example.cpp create mode 100644 example/ooura_fourier_integrals_multiprecision_example.cpp diff --git a/example/ooura_fourier_integrals_cosine_example.cpp b/example/ooura_fourier_integrals_cosine_example.cpp new file mode 100644 index 000000000..2f39a70e0 --- /dev/null +++ b/example/ooura_fourier_integrals_cosine_example.cpp @@ -0,0 +1,78 @@ +// Copyright Paul A. Bristow, 2019 +// Copyright Nick Thompson, 2019 + +// 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) + +// This example requires C++17. + +//#define BOOST_MATH_INSTRUMENT_OOURA // or -DBOOST_MATH_INSTRUMENT_OOURA etc for diagnostic output. + +#include // + +#include // For pi (including for multiprecision types, if used.) + +#include +#include +#include +#include + +int main() +{ + + std::cout.precision(std::numeric_limits::max_digits10); // Show all potentially significant digits. + + using boost::math::quadrature::ooura_fourier_cos; + using boost::math::constants::half_pi; + using boost::math::constants::e; + + // constexpr double double_tol = 10 * std::numeric_limits::epsilon(); // Tolerance. + + //[ooura_fourier_integrals_cosine_example_1 + auto integrator = ooura_fourier_cos(); + // Use the default tolerance root_epsilon and eight levels for type double. + + auto f = [](double x) + { // More complex example function. + return 1 / (x * x + 1); + }; + + double omega = 1; + + auto [result, relative_error] = integrator.integrate(f, omega); + std::cout << "Integral = " << result << ", relative error estimate " << relative_error << std::endl; + + //] [/ooura_fourier_integrals_cosine_example_1] + + //[ooura_fourier_integrals_cosine_example_2 + + constexpr double expected = half_pi() / e(); + std::cout << "pi/(2e) = " << expected << ", difference " << result - expected << std::endl; + //] [/ooura_fourier_integrals_cosine_example_2] + +} // int main() + +/* + +//[ooura_fourier_integrals_example_cosine_output_1 + +Integral = 0.57786367489546109, relative error estimate 6.4177395404415149e-09 +pi/2 = 0.57786367489546087, difference 2.2204460492503131e-16 + +//] [/ooura_fourier_integrals_example_cosine_output_1] + + +//[ooura_fourier_integrals_example_cosine_diagnostic_output_1 + +h = 1.000000000000000, I_h = 0.588268622591776 = 0x1.2d318b7e96dbe00p-1, absolute error estimate = nan +h = 0.500000000000000, I_h = 0.577871642184837 = 0x1.27decab8f07b200p-1, absolute error estimate = 1.039698040693926e-02 +h = 0.250000000000000, I_h = 0.577863671186883 = 0x1.27ddbf42969be00p-1, absolute error estimate = 7.970997954576120e-06 +h = 0.125000000000000, I_h = 0.577863674895461 = 0x1.27ddbf6271dc000p-1, absolute error estimate = 3.708578555361441e-09 +Integral = 5.778636748954611e-01, relative error estimate 6.417739540441515e-09 +pi/2 = 5.778636748954609e-01, difference 2.220446049250313e-16 + +//] [/ooura_fourier_integrals_example_cosine_diagnostic_output_1] + +*/ \ No newline at end of file diff --git a/example/ooura_fourier_integrals_example.cpp b/example/ooura_fourier_integrals_example.cpp new file mode 100644 index 000000000..8574c74e8 --- /dev/null +++ b/example/ooura_fourier_integrals_example.cpp @@ -0,0 +1,77 @@ +// Copyright Paul A. Bristow, 2019 +// Copyright Nick Thompson, 2019 + +// 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) + +// This example requires C++11. + +#define BOOST_MATH_INSTRUMENT_OOURA // or -DBOOST_MATH_INSTRUMENT_OOURA etc for diagnostics. + +#include // + +#include // For pi (including for multiprecision types, if used.) + +#include +#include +#include +#include + +int main() +{ + + std::cout.precision(std::numeric_limits::max_digits10); // Show all potentially significant digits. + + using boost::math::quadrature::ooura_fourier_sin; + using boost::math::constants::half_pi; + + // constexpr double double_tol = 10 * std::numeric_limits::epsilon(); // Tolerance. + +//[ooura_fourier_integrals_example_1 + ooura_fourier_sinintegrator = ooura_fourier_sin(); + // Use the default tolerance root_epsilon and eight levels for type double. + + auto f = [](double x) + { // Simple reciprocal function for sinc. + return 1 / x; + }; + + double omega = 1; + std::pair result = integrator.integrate(f, omega); + std::cout << "Integral = " << result.first << ", relative error estimate " << result.second << std::endl; + +//] [/ooura_fourier_integrals_example_1] + +//[ooura_fourier_integrals_example_2 + + constexpr double expected = half_pi(); + std::cout << "pi/2 = " << expected << ", difference " << result.first - expected << std::endl; +//] [/ooura_fourier_integrals_example_2] + +} // int main() + +/* + +//[ooura_fourier_integrals_example_output_1 + +integral = 1.5707963267948966, relative error estimate 1.2655356398390254e-11 +pi/2 = 1.5707963267948966, difference 0 + +//] [/ooura_fourier_integrals_example_output_1] + + +//[ooura_fourier_integrals_example_diagnostic_output_1 + +ooura_fourier_sin with relative error goal 1.4901161193847656e-08 & 8 levels. +h = 1.000000000000000, I_h = 1.571890732004545 = 0x1.92676e56d853500p+0, absolute error estimate = nan +h = 0.500000000000000, I_h = 1.570793292491940 = 0x1.921f825c076f600p+0, absolute error estimate = 1.097439512605325e-03 +h = 0.250000000000000, I_h = 1.570796326814776 = 0x1.921fb54458acf00p+0, absolute error estimate = 3.034322835882008e-06 +h = 0.125000000000000, I_h = 1.570796326794897 = 0x1.921fb54442d1800p+0, absolute error estimate = 1.987898734512328e-11 +Integral = 1.570796326794897e+00, relative error estimate 1.265535639839025e-11 +pi/2 = 1.570796326794897e+00, difference 0.000000000000000e+00 + +//] [/ooura_fourier_integrals_example_diagnostic_output_1] + +*/ \ No newline at end of file diff --git a/example/ooura_fourier_integrals_multiprecision_example.cpp b/example/ooura_fourier_integrals_multiprecision_example.cpp new file mode 100644 index 000000000..3fdbfeb62 --- /dev/null +++ b/example/ooura_fourier_integrals_multiprecision_example.cpp @@ -0,0 +1,125 @@ +// Copyright Paul A. Bristow, 2019 +// Copyright Nick Thompson, 2019 + +// 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) + +// This example math/example/ooura_fourier_integrals_multiprecision_example/ooura_fourier_integrals_multiprecision_example.cpp requires C++17. + +#define BOOST_MATH_INSTRUMENT_OOURA // or -DBOOST_MATH_INSTRUMENT_OOURA etc for diagnostic output. + +#include // +#include // for cpp_bin_float_quad, cpp_bin_float_50... +#include // For pi (including for multiprecision types, if used.) + +#include +#include +#include +#include +#include + +//typedef boost::multiprecision::cpp_bin_float_50 Real; + +int main() +{ + try + { + typedef boost::multiprecision::cpp_bin_float_quad Real; + + std::cout.precision(std::numeric_limits::max_digits10); // Show all potentially significant digits. + + using boost::math::quadrature::ooura_fourier_cos; + using boost::math::constants::half_pi; + using boost::math::constants::e; + + //[ooura_fourier_integrals_multiprecision_example_1 + + // Use the default parameters for tolerance root_epsilon and eight levels for the type. + //auto integrator = ooura_fourier_cos(); + // Decide a (tight) tolerance. + //const Real tol = 2 * std::numeric_limits::epsilon(); // Tolerance. + const Real tol = 1 * std::numeric_limits::epsilon(); + auto integrator = ooura_fourier_cos(tol, 8); // Loops for 16 + + auto f = [](Real x) + { // More complex example function. + return 1 / (x * x + 1); + }; + + double omega = 1; + + auto [result, relative_error] = integrator.integrate(f, omega); + std::cout << "Integral = " << result << ", relative error estimate " << relative_error << std::endl; + + //] [/ooura_fourier_integrals_multiprecision_example_1] + + //[ooura_fourier_integrals_multiprecision_example_2 + + const Real expected = half_pi() / e(); + std::cout << "pi/(2e) = " << expected << ", difference " << result - expected << std::endl; + //] [/ooura_fourier_integrals_multiprecision_example_2] + } + catch (std::exception ex) + { + // Lacking try& catch blocks, the program will abort, whereas the + //message below from the thrown exception will give some helpful clues as to the cause of the problem. + + std::cout << "\n""Message from thrown exception was:\n " << ex.what() << std::endl; + } + +} // int main() + +/* + +//[ooura_fourier_integrals_example_multiprecision_output_1 + +Integral = 0.5778636748954608589550465916563501587, relative error estimate 4.609814684522163895264277312610830278e-17 +pi/(2e) = 0.5778636748954608659545328919193707407, difference -6.999486300263020581921171645255733758e-18 + + +//] [/ooura_fourier_integrals_example_multiprecision_output_1] + + +//[ooura_fourier_integrals_example_multiprecision_diagnostic_output_1 + +ooura_fourier_cos with relative error goal 3.851859888774471706111955885169854637e-34 & 15 levels. +epsilon for type = 1.925929944387235853055977942584927319e-34 +h = 1.000000000000000000000000000000000, I_h = 0.588268622591776615359568690603776 = 0.5882686225917766153595686906037760, absolute error estimate = nan +h = 0.500000000000000000000000000000000, I_h = 0.577871642184837461311756940493259 = 0.5778716421848374613117569404932595, absolute error estimate = 1.039698040693915404781175011051656e-02 +h = 0.250000000000000000000000000000000, I_h = 0.577863671186882539559996800783122 = 0.5778636711868825395599968007831220, absolute error estimate = 7.970997954921751760139710137450075e-06 +h = 0.125000000000000000000000000000000, I_h = 0.577863674895460885593491133506723 = 0.5778636748954608855934911335067232, absolute error estimate = 3.708578346033494332723601147051768e-09 +h = 0.062500000000000000000000000000000, I_h = 0.577863674895460858955046591656350 = 0.5778636748954608589550465916563502, absolute error estimate = 2.663844454185037302771663314961535e-17 +h = 0.031250000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563484, absolute error estimate = 1.733336949948512267750380148326435e-33 +h = 0.015625000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563479, absolute error estimate = 4.814824860968089632639944856462318e-34 +h = 0.007812500000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563473, absolute error estimate = 6.740754805355325485695922799047246e-34 +h = 0.003906250000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563475, absolute error estimate = 1.925929944387235853055977942584927e-34 +Integral = 5.778636748954608589550465916563475e-01, relative error estimate 3.332844800697411177051445985473052e-34 +pi/(2e) = 5.778636748954608589550465916563481e-01, difference -6.740754805355325485695922799047246e-34 + + +>ooura_fourier_cos with relative error goal 1.925929944387235853055977942584927319e-34 & 15 levels. +1>epsilon for type = 1.925929944387235853055977942584927319e-34 +1>h = 1.000000000000000000000000000000000, I_h = 0.588268622591776615359568690603776 = 0.5882686225917766153595686906037760, absolute error estimate = nan +1>h = 0.500000000000000000000000000000000, I_h = 0.577871642184837461311756940493259 = 0.5778716421848374613117569404932595, absolute error estimate = 1.039698040693915404781175011051656e-02 +1>h = 0.250000000000000000000000000000000, I_h = 0.577863671186882539559996800783122 = 0.5778636711868825395599968007831220, absolute error estimate = 7.970997954921751760139710137450075e-06 +1>h = 0.125000000000000000000000000000000, I_h = 0.577863674895460885593491133506723 = 0.5778636748954608855934911335067232, absolute error estimate = 3.708578346033494332723601147051768e-09 +1>h = 0.062500000000000000000000000000000, I_h = 0.577863674895460858955046591656350 = 0.5778636748954608589550465916563502, absolute error estimate = 2.663844454185037302771663314961535e-17 +1>h = 0.031250000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563484, absolute error estimate = 1.733336949948512267750380148326435e-33 +1>h = 0.015625000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563479, absolute error estimate = 4.814824860968089632639944856462318e-34 +1>h = 0.007812500000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563473, absolute error estimate = 6.740754805355325485695922799047246e-34 +1>h = 0.003906250000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563475, absolute error estimate = 1.925929944387235853055977942584927e-34 +1>h = 0.001953125000000000000000000000000, I_h = 0.577863674895460858955046591656346 = 0.5778636748954608589550465916563463, absolute error estimate = 1.155557966632341511833586765550956e-33 +1>h = 0.000976562500000000000000000000000, I_h = 0.577863674895460858955046591656350 = 0.5778636748954608589550465916563504, absolute error estimate = 4.140749380432557084070352576557594e-33 +1>h = 0.000488281250000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563478, absolute error estimate = 2.600005424922768401625570222489652e-33 +1>h = 0.000244140625000000000000000000000, I_h = 0.577863674895460858955046591656342 = 0.5778636748954608589550465916563418, absolute error estimate = 6.066679324819792937126330519142521e-33 +1>h = 0.000122070312500000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563467, absolute error estimate = 4.911121358187451425292743753591565e-33 +1>h = 0.000061035156250000000000000000000, I_h = 0.577863674895460858955046591656342 = 0.5778636748954608589550465916563424, absolute error estimate = 4.333342374871280669375950370816086e-33 +1>h = 0.000030517578125000000000000000000, I_h = 0.577863674895460858955046591656328 = 0.5778636748954608589550465916563282, absolute error estimate = 1.415558509124618351996143787799922e-32 + + + +//] [/ooura_fourier_integrals_example_multiprecision_diagnostic_output_1] + +*/ \ No newline at end of file