apparent numerical integration bug in sage

[Anthony Kable]

I have run the code

var('x')
f1(x)=1/sqrt(x^3+2)
f2(x)=1/sqrt(x^4+2)
r1=RR(integrate(f1(x),(x,1,10^(10))))
r2=RR(integrate(f2(x),(x,1,10^(10))))
s1=RR(integrate(f1(x),(x,1,10^(11))))
s2=RR(integrate(f2(x),(x,1,10^(11))))
integrals_for_comparison=[[r1,r2],[s1,s2]]
integrals_for_comparison

on Sage Version 6.2, Release Date: 2014-05-06 on both a mac running OSX 10.9 and an old hp compaq running ubuntu 14.04, using
the notebook interface in both cases with Safari as the browser on the mac and Firefox as the browser on the compaq.

In both cases, I obtained the result

[[1.82808026966879, 0.881569060342641], [1.82809394511342,
2.32149165988606e-7]].

The integrals to 10^(10) are being evaluated more or less correctly, as is the integral of f1 to 10^(11), but the
integral of f2 to 10^(11) is wrong by about seven orders of magnitude.

Despite my question title, I do not actually know enough about how sage handles this computation to know whether it is a
numerical integration bug or a bug in the numerical evaluation of elliptic functions, but something is obviously wrong.

[kcrisman]

f1(x)=1/sqrt(x^3+2)
f2(x)=1/sqrt(x^4+2)
r1=RR(integrate(f1(x),(x,1,10^(10))))
r2=RR(integrate(f2(x),(x,1,10^(10))))
s1=RR(integrate(f1(x),(x,1,10^(11))))
s2=RR(integrate(f2(x),(x,1,10^(11))))

Note that probably using something like

sage: numerical_integral(f2,1,10^8)

(0.8815690504421161, 3.309409685784312e-09)

would be better here, since Sage isn't doing the symbolic integration in any case (Maxima is fairly weak on this kind of integral).

The integrals to 10^(10) are being evaluated more or less correctly, as is the integral of f1 to 10^(11), but the
integral of f2 to 10^(11) is wrong by about seven orders of magnitude.

sage: numerical_integral(f2,1,10^8)

(0.8815690504421161, 3.309409685784312e-09)

sage: numerical_integral(f2,1,10^9)

(0.8815690594421439, 2.7280605832086615e-08)

sage: numerical_integral(f2,1,10^10)

(0.8815690603426408, 6.194229607849825e-07)

sage: numerical_integral(f2,1,4.5*10^10)

(0.8815690604198958, 2.5079492928729825e-11)

sage: numerical_integral(f2,1,10^11)

(2.3214916598860602e-07, 4.5569931705290324e-07)

Yeah, that is annoying.  Something similar happens with our other "standard" numerical integration procedure:

sage: f2.nintegrate(x,1,1000000)

(0.8815680604421181, 1.6586910832421555e-12, 819, 0)

sage: f2.nintegrate(x,1,10000000)

(-1.0000000002652395e-07, 9.9265675869388e-15, 861, 5)

Interestingly, going to infinity avoids this:

sage: numerical_integral(f2,1,oo)

(0.8815690604419927, 5.603223703062511e-08)

Or use 

sage: numerical_integral(f2,10^10,10^11)

(9.000000000000001e-11, 9.99200722162641e-25)

This uses the Gnu Scientific Library.  Perhaps this is a limitation of how it is constructed.  But I have a feeling someone else who knows more about the internals of quadrature methods will have more info, so for now I've opened http://trac.sagemath.org/ticket/16905 for this bug.  Thank you!

[Robert Dodier]

On 2014-08-29, kcrisman <kcri...@gmail.com> wrote:

> sage: numerical_integral(f2,1,10^8)
> (0.8815690504421161, 3.309409685784312e-09)
> sage: numerical_integral(f2,1,10^9)
> (0.8815690594421439, 2.7280605832086615e-08)
> sage: numerical_integral(f2,1,10^10)
> (0.8815690603426408, 6.194229607849825e-07)
> sage: numerical_integral(f2,1,4.5*10^10)
> (0.8815690604198958, 2.5079492928729825e-11)
> sage: numerical_integral(f2,1,10^11)
> (2.3214916598860602e-07, 4.5569931705290324e-07)

> This uses the Gnu Scientific Library. Perhaps this is a limitation of how
> it is constructed. But I have a feeling someone else who knows more about
> the internals of quadrature methods will have more info, so for now I've
> opened http://trac.sagemath.org/ticket/16905 for this bug. Thank you!

Sage punts numerical integrals to QUADPACK or a translation of it,
right? QUADPACK is based on Gauss-Kronrod rules which are essentially
Gaussian integration + an efficient refinement scheme. To determine
whether to refine the estimate, the integrand is evaluated at some
points. The difficulty with 1/sqrt(x^4 + 2) is that it is nearly zero
for much of the interval (1, 10^11). If it is close enough to zero
at the evaluation points, the refinement heuristic can be fooled.
Probably other quadrature methods can be fooled in a similar way.

I guess this is a bug in the sense that the result is far from correct
even though the code is a correct implementation of the method.
I wonder what other heuristic could be invented -- maybe since we
can see that integrals on shorter intervals are much larger and
the integrand is nonnegative, the integral on the whole interval
must be at least that large -- perhaps that can be formalized.

I am trying the examples with QUADPACK functions in Maxima --
I find that quad_qags(f2, x, 1, 10^11) fails (with error=5, "integral
is probably divergent or slowly convergent") but
quad_qag(f2, x, 1, 10^11, 4) succeeds, likewise quad_qagi(f2, x, 1, inf)
succeeds. If Sage is indeed calling QUADPACK, perhaps at least the
error number can be reported?

For what it's worth,

Robert Dodier

[Harald Schilly]

On Friday, August 29, 2014 7:15:27 PM UTC+2, Robert Dodier wrote:

 QUADPACK ...

I've tried this in mpmath's quad. I think it works there, but maybe I've overlooked the actual problem.

sage: import mpmath as mp
sage: f1 = lambda _ : 1. / mp.sqrt(_^3 + 2)
sage: f2 = lambda _ : 1. / mp.sqrt(_^4 + 2)
sage: mp.quad(f1, [1, 10^9])
mpf('1.8280370240753521')
sage: mp.quad(f1, [1, 10^10])
mpf('1.8280802692664053')
sage: mp.quad(f1, [1, 4.5 * 10^10])
mpf('1.8280908397677311')
sage: mp.quad(f1, [1, 10^11])
mpf('1.8280939410898087')
sage: mp.quad(f1, [1, 10^12])
mpf('1.8280982294340853')
sage: mp.quad(f2, [1, 10^9])
mpf('0.88156905940188868')
sage: mp.quad(f2, [1, 10^10])
mpf('0.88156905993967871')
sage: mp.quad(f2, [1, 4.5 * 10^10])
mpf('0.88156905860896162')
sage: mp.quad(f2, [1, 10^11])
mpf('0.88156905640783767')

-- H

[kcrisman]

Sage punts numerical integrals to QUADPACK or a translation of it,

What does GSL use?  I forgot that Scipy also has quadrature, in addition to Maxima... wealth of riches.

 

right? QUADPACK is based on Gauss-Kronrod rules which are essentially
Gaussian integration + an efficient refinement scheme. To determine
whether to refine the estimate, the integrand is evaluated at some
points. The difficulty with 1/sqrt(x^4 + 2) is that it is nearly zero
for much of the interval (1, 10^11). If it is close enough to zero
at the evaluation points, the refinement heuristic can be fooled.
Probably other quadrature methods can be fooled in a similar way.

I guess this is a bug in the sense that the result is far from correct
even though the code is a correct implementation of the method.
I wonder what other heuristic could be invented -- maybe since we
can see that integrals on shorter intervals are much larger and
the integrand is nonnegative, the integral on the whole interval
must be at least that large -- perhaps that can be formalized.

I am trying the examples with QUADPACK functions in Maxima --
I find that quad_qags(f2, x, 1, 10^11) fails (with error=5, "integral
is probably divergent or slowly convergent") but
quad_qag(f2, x, 1, 10^11, 4) succeeds, likewise quad_qagi(f2, x, 1, inf)
succeeds. If Sage is indeed calling QUADPACK, perhaps at least the
error number can be reported?

That would be the "nintegrate()" method on symbolic expressions, and indeed

sage: f2.nintegrate(x,1,1000000)

(0.8815680604421181, 1.6586910832421555e-12, 819, 0)

sage: f2.nintegrate(x,1,10000000)

(-1.0000000002652395e-07, 9.9265675869388e-15, 861, 5)

so the "5" in the latter is the error code, which we even document as

      - ``5`` - integral is probably divergent or slowly

        convergent

Also, gp seems to handle it.

sage: gp.eval('intnum(x=1,100000000000,1/sqrt(x^4+2))')

'0.88156906043147435374520375967552406680'

sage: gp.eval('intnum(x=1,1000000000000,1/sqrt(x^4+2))')

'0.88156906044791138558085421922579474969'

Harald, thanks.  I think now we have six or seven ways to do quadrature in Sage!

According to https://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html#Numerical-Integration GSL does indeed reimplement QUADPACK, so I'm not sure what we should do here.

[Robert Dodier]

On 2014-08-29, kcrisman <kcri...@gmail.com> wrote:

> sage: gp.eval('intnum(x=1,100000000000,1/sqrt(x^4+2))')
> '0.88156906043147435374520375967552406680'
> sage: gp.eval('intnum(x=1,1000000000000,1/sqrt(x^4+2))')
> '0.88156906044791138558085421922579474969'

Hmm, what method does PARI/GP use? The documentation for intnum
doesn't seem to mention any algorithms. ... I just looked at the
source code (intnum.c) and I can't tell what's going on. There is
some code for Romberg's method (intnumromb) but it's not called from
intnum.

> Harald, thanks. I think now we have six or seven ways to do quadrature in
> Sage!

Yes, but most of them are QUADPACK, right?

best

Robert Dodier

[Peter Bruin]

Hello Robert,

Hmm, what method does PARI/GP use? The documentation for intnum
doesn't seem to mention any algorithms. ... I just looked at the
source code (intnum.c) and I can't tell what's going on. There is
some code for Romberg's method (intnumromb) but it's not called from
intnum

From the PARI user's manual, chapter 3: "Starting with version 2.2.9 the ``double exponential'' univariate integration method is implemented in intnum and its variants. Romberg integration is still available under the name intnumromb, but superseded."

Peter

[Harald Schilly]

On Sat, Aug 30, 2014 at 12:35 AM, Robert Dodier <robert...@gmail.com> wrote:
>> Harald, thanks. I think now we have six or seven ways to do quadrature in
>> Sage!
>
> Yes, but most of them are QUADPACK, right?

The whole idea of mpmath is to be standalone, right?
https://github.com/fredrik-johansson/mpmath/blob/master/mpmath/calculus/quadrature.py#L385
Fredrik could explain it in more detail …

--H

[John Cremona]

Pari's numerical integration does implement some very sophisticated algorithms, as implemented by Henri Cohen, who has given talks about how good they are but which stress how important it is to use the right variant, depending in the analytic properties of the function being integrated.  If you use the wrong one your get garbage with no warnings -- as HC happily demonstrated in his talks.  [This information may be out of date: it is possible that pari does try to intelligently determine which variant to use.  I don't know.]