The rate of convergence of the tanh-sinh integrator depends on the minimum distance of the interval of integration to the nearest singularity in the complex plane. In your case, this distance is zero, so the usual convergence theorems do not apply. This doesn't mean you'll do better with another quadrature method; almost all rapidly convergent quadrature routines have a "Bernstein ellipse" hiding somewhere in the background.

Under normal circumstances (finite distance to nearest singularity in the complex plane), we get exponential convergence.
Let Q_h[f] be the value of the tanh-sinh quadrature computed at stepsize h.
If

|Q_h[f] - Q_{h/2}[f]| ≤ √ε|Q_h[f]|,

i.e., we have *observed* a relative accuracy on the order √ε, then in fact Q_h[f] is accurate to relative accuracy ε.
So the default is for the tanh-sinh integrator to terminate once a relative accuracy of √ε has been observed.

    Nick


‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐
On Friday, August 21, 2020 6:22 PM, Anirban Pal via Boost-users <boost-users@lists.boost.org> wrote:

Thank you Kila. I used your code and yes, there was no problem.

With respect to your statement on tolerance being \sqrt{\epsilon}. Is this true for all functions or only functions involving sqrt.

For example, if I integrate the function below from 1 to 2, the error is of the O(10^-201) rather than O(10^-101) for a precision of cpp_dec_float<200>.

auto f3 = [](BF x) { return log(x-BF("1"))*log(BF("2.0")-x); };
BF Q3 = integrator.integrate(f3, BF("1"), BF("2"));

BF pi = boost::math::constants::pi<BF>();
BF Q3_true = BF(2)-pi*pi/BF(6);

BF err = Q3-Q3_true;

So, I'm wondering if the error is dependent on the function or the type of singularity at the end points?

From the boost math documentation here, it seems like the derivatives of f(x) are relevant for the Newton Raphson scheme.
I'm not sure if the integrator is using similar schemes ..

For the log(x-1)*log(2-x) function, the derivative involves terms of the order of x near the singularity.
For the sqrt function, the derivative involves terms of the order of sqrt(x).

Is this reasoning correct? Or maybe this is too simplistic. I'm not really familiar with the math behind the integrator.


On Fri, 21 Aug 2020 at 07:31, kila suelika via Boost-users <boost-users@lists.boost.org> wrote:
First problem - more precision
Check the stopping condition of boost.math, it says the tolerance is \sqrt{\epsilon}. So if you use 100-digits precision, then tolerance is roughly 10^-50. For more precision, you have to use more digits such as 1000 digits by using BF= number<cpp_dec_float<1000> >; . Full code:

#include<iostream>
#include<boost/math/quadrature/tanh_sinh.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

using namespace std;
using namespace boost::math::quadrature;
using namespace boost::multiprecision;
using BF= number<cpp_dec_float<1000> >;

int main(){
auto f2 = [](BF x) { return x/BF(sqrt(x*x-BF("0.25"))); };
tanh_sinh<BF> integrator;

std::cout << std::setprecision(std::numeric_limits<BF>::max_digits10)
     << "Comp.: " << integrator.integrate(f2, BF("0.5"), sqrt(BF("1.25"))) << std::endl;

    return 0;
};

Be patient, this run about 10 mins in my single core server.

Second problem - error
I didn't reproduce it. Just modify my code and everything is OK.

On Fri, Aug 21, 2020 at 9:18 AM Anirban Pal via Boost-users <boost-users@lists.boost.org> wrote:

Thank you for your response. I’ve run into a few more issues with accuracy. Here is my function, which has a singularity at 0.5 and I'm integrating it from 0.5 to sqrt(1.25):

using BF = boost::multiprecision::cpp_bin_float_100;

auto f1 = [](BF x) { return x/BF(sqrt(x*x-BF("0.25"))); };
BF Q1 = integrator.integrate(f1, BF("0.5"), BF(sqrt(BF("1.25"))));

This gives me a value of Q1 which is close to the true value of 1 with an error of O(10^-51). I am curious if this should be lower, like O(10^-100). (I hope I’m not being greedy!)

Nevertheless, if I try a related function.

auto f2 = [](BF x) { return x/BF(sqrt(x*x-BF("0.25"))) - x; };
BF Q2 = integrator.integrate(f2, BF("0.5"), BF(sqrt(BF("1.25"))));

or an equivalent one:

auto f3 = [](BF x) { return (x - x*BF(sqrt(x*x-BF("0.25"))) )/BF(sqrt(x*x-BF("0.25"))); };
BF Q3 = integrator.integrate(f3, BF("0.5"), BF(sqrt(BF("1.25"))));

I get an error message:

boost_test5: /usr/local/include/boost/math/quadrature/tanh_sinh.hpp:196: boost::math::quadrature::tanh_sinh<Real, Policy>::integrate(F, Real, Real, Real, Real*, Real*, std::size_t*)::<lambda(Real, Real)> [with F = main(int, char**)::<lambda(BF)>; Real = boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<100> >; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>; result_type = boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<100> >]: Assertion `position != a' failed.
Aborted (core dumped)

I'm not sure how I could improve my code. I tried using boost::multiprecision::sqrt but i didn't notice any difference.

Thank you for your attention.



On Mon, 17 Aug 2020 at 04:47, Paul A. Bristow via Boost-users boost-users@lists.boost.org wrote:





> -----Original Message-----
> From: Boost-users <boost-users-bounces@lists.boost.org> On Behalf Of John Maddock via Boost-users
> Sent: 16 August 2020 21:24
> To: Anirban Pal via Boost-users <boost-users@lists.boost.org>
> Cc: John Maddock <jz.maddock@googlemail.com>
> Subject: Re: [Boost-users] limited accuracy for tanh-sinh integrator
>
>
> On 16/08/2020 21:15, Anirban Pal via Boost-users wrote:
> >
> > Hello,
> >
> > I’m using the tanh-sinh integrator to integrate a simple function f(x)
> > = 0.26*x from 3.0 to 4.0.
> > The exact result is 0.91. With the integrator I’m getting a result
> > accurate to only 10^-18 with cpp_bin_float_100 multiprecision.
> >
> You have the double precision constant 0.26 in your code, and since this is an inexact binary value,
> everything that depends on that constant is inherently limited to double precision.  Try constructing it
> as Real(26) / 100.

Everyone falls into this pit, and in my case, repeatedly ☹

You need to be ever vigilant not to use double precision items by mistake.

Using numeric_limits<>::max_digits10 is useful because it shows the digits beyond double at about 17 decimal digits.  Being random should raise suspicions.

If you are only using multiprecision types, then always using decimal digit strings like "1.23456" is simple.

If you can use a fraction like Real(1)/1000 that also works fine.

You may also find the macro BOOST_MATH_TEST_VALUE  in I:\boost\libs\math\test\test_value.hpp helpful if you want to make code portable over different floating-point types from fundamental to multiprecision?

Paul


_______________________________________________
Boost-users mailing list




--
Anirban Pal



_______________________________________________
Boost-users mailing list
_______________________________________________
Boost-users mailing list
Boost-users@lists.boost.org
https://lists.boost.org/mailman/listinfo.cgi/boost-users


--
Anirban Pal