Opened 12 years ago
Last modified 2 years ago
#11164 needs_info defect
Integral of sin(x)/x gives false result.
Reported by: | Ben Reynwar | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | |
Component: | calculus | Keywords: | maxima, integration |
Cc: | Merged in: | ||
Authors: | Reviewers: | ||
Report Upstream: | Reported upstream. Developers acknowledge bug. | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: | todo |
Description (last modified by )
> eq = sin(x)/x > integrate(eq, x, -1e-6, 1e-6) -3.14159065359
I expected an answer of about 2e-6.
Attachments (1)
Change History (24)
comment:1 Changed 12 years ago by
Component: | PLEASE CHANGE → calculus |
---|---|
Owner: | changed from tbd to Burcin Erocal |
Priority: | minor → major |
Report Upstream: | N/A → Fixed upstream, in a later stable release. |
Changed 11 years ago by
Attachment: | trac_11164_verify_integration_sinx_x.patch added |
---|
add doctest to confirm fix
comment:3 Changed 11 years ago by
Status: | needs_review → 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 11 years ago by
Status: | needs_work → 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: 6 Changed 11 years ago by
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 Changed 11 years ago by
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 11 years ago by
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: 9 Changed 10 years ago by
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 Changed 9 years ago by
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 8 years ago by
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: 12 Changed 8 years ago by
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 follow-up: 13 Changed 8 years ago by
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 computeintegrate(sin(x)/x, x, -1, 1)
, the two terms coming fromexp(%i*x)
andexp(-%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 ofsin(x)/x
, which is an entire function.
How much you wanna bet this is the same underlying problem as in #17328?
comment:13 follow-up: 14 Changed 8 years ago by
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 follow-up: 15 Changed 8 years ago by
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).
comment:15 follow-up: 17 Changed 8 years ago by
Report Upstream: | Fixed upstream, in a later stable release. → 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/
comment:16 Changed 8 years ago by
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 Changed 8 years ago by
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 follow-up: 19 Changed 7 years ago by
Stopgaps: | → todo |
---|
comment:19 Changed 3 years ago by
Replying to jakobkroeker:
We really need Maxima to provide a proper
Si
function
Maxima has provided a proper Si function for a while called expintegral_si(x)
which seems to fit the bill. Maxima doesn't seem to use it though, and the behavior you described still occurs:
sage: maxima_calculus('integrate(sin(x)/x,x)')
-(%i*gamma_incomplete(0,%i*x)-%i*gamma_incomplete(0,-%i*x))/2
comment:20 Changed 3 years ago by
Keywords: | maxima added |
---|---|
Report Upstream: | Reported upstream. No feedback yet. → Reported upstream. Developers acknowledge bug. |
comment:21 Changed 3 years ago by
Owner: | Burcin Erocal deleted |
---|
comment:22 Changed 3 years ago by
Keywords: | integration added |
---|
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.
We'd need doctesting to test both of these once we upgrade Maxima.