Opened 6 years ago

Last modified 4 years ago

#14274 new defect

Numerical approximation of a divergent integral

Reported by: eviatarbach Owned by: burcin
Priority: major Milestone: sage-6.4
Component: calculus Keywords:
Cc: Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by tmonteil)

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)

Change History (11)

comment:1 Changed 6 years ago by eviatarbach

  • Description modified (diff)

comment:2 Changed 6 years ago by kcrisman

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!

comment:3 Changed 6 years ago by eviatarbach

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).

comment:4 Changed 6 years ago by eviatarbach

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.

comment:5 Changed 6 years ago by kcrisman

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

comment:6 Changed 6 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:7 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:8 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:9 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:10 Changed 5 years ago by rws

Edit: sorry wrong ticket

Last edited 5 years ago by rws (previous) (diff)

comment:11 Changed 4 years ago by tmonteil

  • Description modified (diff)
Note: See TracTickets for help on using tickets.