Opened 9 years ago
Last modified 4 months ago
#11164 needs_info defect
Integral of sin(x)/x gives false result.
Reported by: | benreynwar | 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 (23)
comment:1 Changed 9 years ago by
- 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.
comment:2 Changed 8 years ago by
- Status changed from new to needs_review
Might as well doctest and close.
comment:3 Changed 8 years ago by
- 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 8 years ago by
- 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: ↓ 6 Changed 8 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 in reply to: ↑ 5 Changed 8 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 8 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 7 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 in reply to: ↑ 8 Changed 6 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 5 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 5 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 in reply to: ↑ 11 ; follow-up: ↓ 13 Changed 5 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 in reply to: ↑ 12 ; follow-up: ↓ 14 Changed 5 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 in reply to: ↑ 13 ; follow-up: ↓ 15 Changed 5 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 in reply to: ↑ 14 ; follow-up: ↓ 17 Changed 5 years ago by
- 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/
comment:16 Changed 5 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 in reply to: ↑ 15 Changed 5 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 4 years ago by
- Stopgaps set to todo
comment:19 in reply to: ↑ 18 Changed 5 months 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
Edit: upstream said this would be fixed here in 2016.
comment:20 Changed 5 months ago by
- Keywords maxima added
- Report Upstream changed from Reported upstream. No feedback yet. to Reported upstream. Developers acknowledge bug.
comment:21 Changed 4 months ago by
- Owner changed from burcin to (none)
comment:22 Changed 4 months 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.