Numerical approximation of a divergent integral

[eviatarbach]

Sage is numerically approximating this integral, even though it's divergent:

sage: integrate(x^3/sqrt(x^7+1), x, 1, oo).n()              
-2.0585298599983344

It seems that if we don't allow Maxima to detect its divergence (numerical_integral passes it directly to GSL), GSL will also fail on simpler divergent integrals:

sage: numerical_integral(1/x, 1, oo)
(65.94931131932763, 8.156214940519742)
sage: numerical_integral(x,1,oo)                    
(-0.4999999993521234, 1.3615531480049015e-09)

See also this ask question:

sage: numerical_integral(e^(-x)/x,0,oo)
(37.191280375549404, 6.239196965189217)

Component: calculus

Issue created by migration from https://trac.sagemath.org/ticket/14274

[eviatarbach]

Description changed:

--- 
+++ 
@@ -3,13 +3,13 @@
 ```
 sage: integrate(x^3/sqrt(x^7+1), x, 1, oo).n()              
 -2.0585298599983344
-sage: numerical_integral(x^3/sqrt(x^7+1), 1, oo)    
-(-2.0585298599983344, 1.1214140727133781e-11)
 ```
 
-`numerical_integral` also fails on simpler divergent integrals:
+It seems that if we don't allow Maxima to detect its divergence (`numerical_integral` passes it directly to GSL), GSL will also fail on simpler divergent integrals:
 
 ```
 sage: numerical_integral(1/x, 1, oo)
 (65.94931131932763, 8.156214940519742)
+sage: numerical_integral(x,1,oo)                    
+(-0.4999999993521234, 1.3615531480049015e-09)
 ```

[kcrisman]

comment:2

Regarding the first one; basically, when we return a noun form

sage: integrate(x^3/sqrt(x^7+1), x, 1, oo)
integrate(x^3/sqrt(x^7 + 1), x, 1, +Infinity)

and call n we do

sage: A._convert(RR)
-2.0585298599985333

because we do

sage: sage.symbolic.integration.integral.DefiniteIntegral._evalf_??

        from sage.gsl.integration import numerical_integral
        # The gsl routine returns a tuple, which also contains the error.
        # We only return the result.
        return numerical_integral(f, a, b)[0]

So these are both manifestations of the same thing.

---

So... is it user error to numerically integrate a divergent integral? I certainly don't know that we should be checking every numerical integral for divergence, particularly since Maxima apparently can't (yet) do the first one in any case!

[eviatarbach]

comment:3

Maybe the options being passed to GSL could be changed? It seems absurd that it should give a numerical answer for numerical_integral(x^3,1,oo), for example. It apparently has an error for divergence (http://www.gnu.org/software/gsl/manual/html_node/Numerical-integration-error-codes.html).

[eviatarbach]

comment:4

In fact, with the code at http://www.physics.ohio-state.edu/~ntg/780/gsl_examples/qagiu_test.cpp, it returns errors for all these integrals.

[kcrisman]

comment:5

If you present a patch, I think we'd be very interested in reviewing it. Silly to return nonsense in these cases.