Opened 8 years ago

Last modified 3 years ago

#11164 needs_info defect

Integral of sin(x)/x gives false result.

Reported by: benreynwar Owned by: burcin
Priority: major Milestone:
Component: calculus Keywords:
Cc: Merged in:
Authors: Reviewers:
Report Upstream: Reported upstream. No feedback yet. Work issues:
Branch: Commit:
Dependencies: Stopgaps: todo

Description (last modified by tscrim)

> eq = sin(x)/x
> integrate(eq, x, -1e-6, 1e-6)
-3.14159065359

I expected an answer of about 2e-6.

Attachments (1)

trac_11164_verify_integration_sinx_x.patch (1.1 KB) - added by dsm 7 years ago.
add doctest to confirm fix

Download all attachments as: .zip

Change History (19)

comment:1 Changed 8 years ago by kcrisman

  • Component changed from PLEASE CHANGE to calculus
  • Owner changed from tbd to burcin
  • Priority changed from minor to major
  • Report Upstream changed from N/A to Fixed upstream, in a later stable release.

In the future, be sure to pick a component (for instance, calculus or symbolics); that will help people find it more easily.

There are two problems here.

  • The actual problem - the report unfortunately has a typo.
    sage: integrate(eq,x,-1e-6,1e-6)
    -3.14159065359
    
    This is definitely a problem. However, it's fixed in the latest Maxima. Here is the old one, in Sage:
    Maxima 5.22.1 http://maxima.sourceforge.net
    using Lisp ECL 10.4.1
    <snip>
    (%i2) keepfloat:true;
    (%o2)                                true
    (%i3) integrate(sin(x)/x,x,-1.00000000000000e-6,1.00000000000000e-6)
    ;
    (%o3)                         - 3.141590653589793
    
    And here is 5.24.0. I assume this isn't a problem of which Lisp we are using?
    Maxima 5.24.0 http://maxima.sourceforge.net
    using Lisp SBCL 1.0.24
    (%i3) keepfloat:true;
    (%o3)                                true
    (%i8) integrate(sin(x)/x,x,-1.00000000000000e-6,1.00000000000000e-6);
    
    rat: replaced 9.999999999999999e-7 by 1/1000000 = 9.999999999999999e-7
    
    rat: replaced -9.999999999999999e-7 by -1/1000000 = -9.999999999999999e-7
    
    rat: replaced -9.999999999999999e-7 by -1/1000000 = -9.999999999999999e-7
    
    rat: replaced 9.999999999999999e-7 by 1/1000000 = 9.999999999999999e-7
                                     %i                                %i
    (%o8) %i gamma_incomplete(0, - -------) - %i gamma_incomplete(0, -------)
                                   1000000                           1000000
    
  • A problem revealed, but also fixed in the latest. Old:
    (%i4) integrate(sin(x)/x,x,-1e-6,1e6);
    
    Maxima encountered a Lisp error:
    
     #<a FLOATING-POINT-OVERFLOW>
    
    Automatically continuing.
    To enable the Lisp debugger set *debugger-hook* to nil.
    
    
    New:
    (%i10) integrate(sin(x)/x,x,-1e-6,1e6);
    
    rat: replaced 9.999999999999999e-7 by 1/1000000 = 9.999999999999999e-7
    
    rat: replaced -9.999999999999999e-7 by -1/1000000 = -9.999999999999999e-7
    
    rat: replaced -1000000.0 by -1000000/1 = -1000000.0
    
    rat: replaced 1000000.0 by 1000000/1 = 1000000.0
                                                                           %i
                                                  %i gamma_incomplete(0, -------)
             %i gamma_incomplete(0, 1000000 %i)                          1000000
    (%o10) - ---------------------------------- - -------------------------------
                             2                                   2
                                      %i
           %i gamma_incomplete(0, - -------)
                                    1000000    %i gamma_incomplete(0, - 1000000 %i)
         + --------------------------------- + ------------------------------------
                           2                                    2
    

We'd need doctesting to test both of these once we upgrade Maxima.

Changed 7 years ago by dsm

add doctest to confirm fix

comment:2 Changed 7 years ago by dsm

  • Status changed from new to needs_review

Might as well doctest and close.

comment:3 Changed 7 years ago by kcrisman

  • Status changed from needs_review to needs_work

I may have misstated whether this is "fixed".

sage: integrate(sin(x)/x,x,-1e-6,1e-6)
I*Ei(-1/1000000*I) - I*Ei(1/1000000*I)
sage: integrate(sin(x)/x,x,-1e-6,1e-6).n()
3.14159465358979
sage: numerical_integral(sin(x)/x,-1e-6,1e-6)
(1.9999999999998889e-06, nan)

The graph makes it very obvious the second (numerical) answer is correct. Unsurprisingly, the other answer is off by exactly pi...

Presumably the first answer is "correct" in the sense that it is an alternate representation of the Si sine integral function. But I think there are some branch issues - Ei(I*x) only has an alternate formula for x>0.

So we still need to do something with this. Or at the very least warn about branches. But really, sin(x)/x shouldn't have to worry about this!

comment:4 Changed 7 years ago by dsm

  • Status changed from needs_work to needs_info

You're right, of course. It's silly to start with sin(x)/x and then generalize to something which doesn't reduce to the expected value, but far, far sillier to add a doctest enshrining the goofy behaviour as the expected. :)

Does this work on the current maxima trunk now?

comment:5 follow-up: Changed 7 years ago by kcrisman

This is the previous version (5.25.1):

(%i5) integrate(sin(x)/x,x,-1/1000000,1/1000000);

(%o5) %i*gamma_incomplete(0,-%i/1000000)-%i*gamma_incomplete(0,%i/1000000)

W|A agrees that this is about -pi. (To be precise, 2*Si(1/10^6)-pi.)

(%i4) display2d:false;

(%o4) false
(%i5) integrate(sin(x)/x,x);

(%o5) -(%i*gamma_incomplete(0,%i*x)-%i*gamma_incomplete(0,-%i*x))/2

Here is the issue; this function may or may not be the same as the sine integral. The problem seems to be that although Maxima has the sine integral (expintegral_si), it isn't really connected with the rest of Maxima (?).

What do you think a good solution to this is? I suppose we could open a ticket for it on the Maxima site, but I'm not sure whether the Si there is robust enough to handle all this yet.

comment:6 in reply to: ↑ 5 Changed 7 years ago by kcrisman

Here is the issue; this function may or may not be the same as the sine integral.

By which I mean it isn't the same, of course, but also that (as defined) it isn't even always the same constant away. Antiderivatives can differ by a constant, but when the constant is different depending on what x is, that's annoying. I think the Maxima folks would say this is a feature (for symbolic reasons).

comment:7 Changed 7 years ago by kcrisman

By the way, achrzesz has a great one-liner for a related question:

 maxima('integrate(diff((exp(x)-1)/x,x),x),gamma_expand:true,factor').sage()

Though here it doesn't solve the problem, but maybe it will trigger someone's memory as to another Maxima flag to set...

sage: maxima('integrate(sin(x)/x,x,-1/1000000,1/1000000),gamma_expand:true').sage()
I*expintegral_ei(-1/1000000*I) - I*expintegral_ei(1/1000000*I) # currently not translated into Sage, but see #11143
sage: a = I*Ei(-1/1000000*I)-I*Ei(1/1000000*I)
sage: a.n()
3.14159465358979

comment:8 follow-up: Changed 6 years ago by kcrisman

I wonder whether this is potentially fixed in #13364 - I don't know if we ever reported upstream, or whether we should have...

comment:9 in reply to: ↑ 8 Changed 6 years ago by kcrisman

I wonder whether this is potentially fixed in #13364 - I don't know if we ever reported upstream, or whether we should have...

Apparently not yet.

comment:10 Changed 5 years ago by tscrim

  • Description modified (diff)

In 6.3.beta4 (so with #13973), I get:

sage: eq = sin(x) / x
sage: integrate(eq, x, -1e-6, 1e-6)
-I*Ei(1/1000000*I) + I*Ei(-1/1000000*I)
sage: _.n()
3.14159465358979
sage: integrate(eq, x, -1e-12, 1e-12)
-I*Ei(1/1000000000000*I) + I*Ei(-1/1000000000000*I)
sage: _.n()
3.14159265359179

Something bad happens with the integration computation when crossing over 0 as we have:

sage: integrate(eq, x, 1e-12, 1e-6)
-1/2*I*Ei(1/1000000*I) + 1/2*I*Ei(1/1000000000000*I) - 1/2*I*Ei(-1/1000000000000*I) + 1/2*I*Ei(-1/1000000*I)
sage: _.n()
9.99999000050877e-7

comment:11 follow-up: Changed 5 years ago by pbruin

The problem is caused by the integral being computed via symbolic integration; see also comment:5. Maxima integrates exp(x)/x to -gamma_incomplete(0, -x), which is in principle correct, but only gives a certain branch of the integral. The branch cut chosen for the incomplete Gamma function is along the negative real axis. When using this to compute integrate(sin(x)/x, x, -1, 1), the two terms coming from exp(%i*x) and exp(-%i*x) have branch cuts along the negative and positive imaginary axis, respectively, so the sum has a branch cut along the whole imaginary axis and therefore can't be the antiderivative of sin(x)/x, which is an entire function.

comment:12 in reply to: ↑ 11 ; follow-up: Changed 4 years ago by kcrisman

The problem is caused by the integral being computed via symbolic integration; see also comment:5. Maxima integrates exp(x)/x to -gamma_incomplete(0, -x), which is in principle correct, but only gives a certain branch of the integral. The branch cut chosen for the incomplete Gamma function is along the negative real axis. When using this to compute integrate(sin(x)/x, x, -1, 1), the two terms coming from exp(%i*x) and exp(-%i*x) have branch cuts along the negative and positive imaginary axis, respectively, so the sum has a branch cut along the whole imaginary axis and therefore can't be the antiderivative of sin(x)/x, which is an entire function.

How much you wanna bet this is the same underlying problem as in #17328?

comment:13 in reply to: ↑ 12 ; follow-up: Changed 4 years ago by kcrisman

How much you wanna bet this is the same underlying problem as in #17328?

Nah, I was thrown off by the gammas, I think.

What if we don't allow incomplete gamma to become Ei when the first arg is zero - would that solve the branch problem?

comment:14 in reply to: ↑ 13 ; follow-up: Changed 4 years ago by kcrisman

What if we don't allow incomplete gamma to become Ei when the first arg is zero - would that solve the branch problem?

Looks like if that was once possible, it isn't any more (already done during integration in Maxima, I mean).

Last edited 4 years ago by kcrisman (previous) (diff)

comment:15 in reply to: ↑ 14 ; follow-up: Changed 4 years ago by kcrisman

  • Report Upstream changed from Fixed upstream, in a later stable release. to Reported upstream. No feedback yet.

What if we don't allow incomplete gamma to become Ei when the first arg is zero - would that solve the branch problem?

Looks like if that was once possible, it isn't any more (already done during integration in Maxima, I mean).

Or maybe not;

sage: maxima('integrate(sin(x)/x,x)')
-(%i*gamma_incomplete(0,%i*x)-%i*gamma_incomplete(0,-%i*x))/2
sage: maxima('integrate(sin(x)/x,x), gamma_expand:true')
-(%i*expintegral_ei(%i*x)-%i*expintegral_ei(-%i*x)+2*%pi)/2

so maybe we are already setting gamma_expand:true somewhere nowadays. Reported upstream at https://sourceforge.net/p/maxima/bugs/2842/

Last edited 4 years ago by kcrisman (previous) (diff)

comment:16 Changed 4 years ago by kcrisman

See also the documentation for gamma_expand. The changelog says

Implementation of the Incomplete Gamma function

 New Maxima User function: gamma_incomplete(a,z)

 The following features are implemented:

 - Evaluation for real and complex numbers in double float and
   bigfloat precision
 - Special values for gamma_incomplete(a,0) and gamma_incomplete(a,inf)
 - When $gamma_expand T expand the following expressions:
   gamma_incomplete(0,z)
   gamma_incomplete(n+1/2)
   gamma_incomplete(1/2-n)
   gamma_incomplete(n,z)
   gamma_incomplete(-n,z)
   gamma_incomplete(a+n,z)
   gamma_incomplete(a-n,z)
 - Mirror symmetry
 - Derivative wrt the arguments a and z

--------------------------------------------------------------------------------
Implementation of the Generalized Incomplete Gamma function

 New Maxima User function: gamma_incomplete_generalized(a,z1,z2)

 The following features are implemented:

 - Evaluation for real and complex numbers in double float and
   bigfloat precision
 - Special values for:
   gamma_incomplete_generalized(a,z1,0)
   gamma_incomplete_generalized(a,0,z2),
   gamma_incomplete_generalized(a,z1,inf)
   gamma_incomplete_generalized(a,inf,z2)
   gamma_incomplete_generalized(a,0,inf)
   gamma_incomplete_generalized(a,x,x)
 - When $gamma_expand T and n an integer expand
   gamma_incomplete_generalized(a+n,z1,z2)
 - Implementation of Mirror symmetry
 - Derivative wrt the arguments a, z1 and z2

--------------------------------------------------------------------------------
Implementation of the Regularized Incomplete Gamma function

 New Maxima User function: gamma_incomplete_regularized(a,z)

 The following features are implemented:

 - Evaluation for real and complex numbers in double float and
   bigfloat precision
 - Special values for:
   gamma_incomplete_regularized(a,0)
   gamma_incomplete_regularized(0,z)
   gamma_incomplete_regularized(a,inf)
 - When $gamma_expand T and n a positive integer expansions for
   gamma_incomplete_regularized(n+1/2,z)
   gamma_incomplete_regularized(1/2-n,z)
   gamma_incomplete_regularized(n,z)
   gamma_incomplete_regularized(a+n,z)
   gamma_incomplete_regularized(a-n,z)
 - Derivative wrt the arguments a and z
 - Implementation of Mirror symmetry

comment:17 in reply to: ↑ 15 Changed 4 years ago by kcrisman

What if we don't allow incomplete gamma to become Ei when the first arg is zero - would that solve the branch problem?

Looks like if that was once possible, it isn't any more (already done during integration in Maxima, I mean).

Or maybe not;

Drivel. We could try this, though the branch problem is likely still there with incomplete gamma.

sage: maxima_calculus('integrate(sin(x)/x,x)')
-(%i*gamma_incomplete(0,%i*x)-%i*gamma_incomplete(0,-%i*x))/2
sage: maxima_calculus('integrate(sin(x)/x,x)').sage()
-1/2*I*Ei(I*x) + 1/2*I*Ei(-I*x)

because we do

        if x == 0:
            return -Ei(-y)

in sage/functions/other.py

But unfortunately

sage: i*CC(0).gamma_inc(-i/100000)-i*CC(0).gamma_inc(i/100000)
-3.14157265358979

So same problem there. Yuck. We really need Maxima to provide a proper Si function, or somehow magically translate such combos ourselves with pattern-matching, which sounds horrible too.

comment:18 Changed 3 years ago by jakobkroeker

  • Stopgaps set to todo
Note: See TracTickets for help on using tickets.