Three fourier integrals examples - work-in-progress

2026-01-19 04:22:09 +00:00 · 2019-07-03 18:14:20 +01:00
parent d41bc9d068
commit 74c0de2f01
3 changed files with 280 additions and 0 deletions
--- a/example/ooura_fourier_integrals_cosine_example.cpp
+++ 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 <boost/math/quadrature/ooura_fourier_integrals.hpp>  // 
+
+#include <boost/math/constants/constants.hpp>  // For pi (including for multiprecision types, if used.)
+
+#include <cmath>
+#include <iostream>
+#include <limits>  
+#include <iostream>
+
+int main()
+{
+
+	std::cout.precision(std::numeric_limits<double>::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<double>::epsilon(); // Tolerance. 
+
+ //[ooura_fourier_integrals_cosine_example_1
+	auto integrator = ooura_fourier_cos<double>();
+	// 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<double>() / e<double>();
+	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]
+
+*/
--- a/example/ooura_fourier_integrals_example.cpp
+++ 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 <boost/math/quadrature/ooura_fourier_integrals.hpp>  // 
+
+#include <boost/math/constants/constants.hpp>  // For pi (including for multiprecision types, if used.)
+
+#include <cmath>
+#include <iostream>
+#include <limits>  
+#include <iostream>
+
+int main()
+{
+
+	std::cout.precision(std::numeric_limits<double>::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<double>::epsilon(); // Tolerance. 
+
+//[ooura_fourier_integrals_example_1
+	ooura_fourier_sin<double>integrator = ooura_fourier_sin<double>();
+	// 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<double, double> 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<double>();
+	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]
+
+*/
--- a/example/ooura_fourier_integrals_multiprecision_example.cpp
+++ 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 <boost/math/quadrature/ooura_fourier_integrals.hpp>  // 
+#include <boost/multiprecision/cpp_bin_float.hpp> // for cpp_bin_float_quad, cpp_bin_float_50...
+#include <boost/math/constants/constants.hpp>  // For pi (including for multiprecision types, if used.)
+
+#include <cmath>
+#include <iostream>
+#include <limits>  
+#include <iostream>
+#include <exception>
+
+//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<Real>::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<Real>();
+		// Decide a (tight) tolerance.
+		//const Real tol = 2 * std::numeric_limits<Real>::epsilon(); // Tolerance. 
+		const Real tol = 1 * std::numeric_limits<Real>::epsilon();
+		auto integrator = ooura_fourier_cos<Real>(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<Real>() / e<Real>();
+		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]
+
+*/