Opened 10 years ago

Closed 8 years ago

#12947 closed defect (fixed)

Bug in integrating x*cos(x^3)

Reported by: kcrisman Owned by: burcin
Priority: trivial Milestone: sage-6.3
Component: calculus Keywords:
Cc: mjo Merged in:
Authors: Peter Bruin Reviewers: Ralf Stephan
Report Upstream: N/A Work issues:
Branch: 72706eb (Commits, GitHub, GitLab) Commit: 72706eb6697fba361e138ddaca0765f8415eb8f3
Dependencies: Stopgaps:

Status badges

Description

From this sage-support thread:

sage: numerical_integral(x*cos(x^3),0,.5)
(0.1247560409610376, 1.3850702913602309e-15)
sage: integrate(x*cos(x^3),(x,0,1/2)).n()
-0.0677842754623305

Given

sage: integrate(x*cos(x^3),x)
1/12*((-I*sqrt(3) - 1)*gamma(2/3, -I*x^3) + (I*sqrt(3) - 1)*gamma(2/3, I*x^3))*(x^3)^(1/3)/x
sage: integrate(x*cos(x^3),(x,0,1/2))
-1/3*gamma(2/3) + 1/6*gamma(2/3, -1/8*I) + 1/6*gamma(2/3, 1/8*I)

this is probably something where Maxima is returning too many imaginary things and something isn't cancelling right either there or in Sage proper? We've had trouble with either side of this with incomplete gamma functions before.

Change History (27)

comment:1 Changed 10 years ago by kcrisman

sage: assume(u>0)
sage: integrate(x*cos(x^3),x,0,u)
-1/3*gamma(2/3) + 1/6*gamma(2/3, -I*u^3) + 1/6*gamma(2/3, I*u^3)

To compare the "right" answer, do

sage: f(u) = integrate(x*cos(x^3),x,0,u)
sage: plot(f,(u,0,2*pi))
sage: plot(lambda u: numerical_integral(x*cos(x^3),0,u)[0],(u,0,2*pi))

Note how they don't really look quite similar, not really a branch cut thing? Also, in Maxima

(%i2) display2d:false;
(%o2) false
(%i3) f:integrate(x*cos(x^3),x);
(%o3) (gamma_incomplete(2/3,%i*x^3)+gamma_incomplete(2/3,-%i*x^3))/6
(%i5) diff(f,x);
(%o5) (-3*%i*x*%e^(%i*x^3)/(-1)^(1/6)-3*%i*x*%e^-(%i*x^3)/(-1)^(1/6))/6

comment:2 Changed 10 years ago by kcrisman

Robert Dodier of Maxima suggests the following later in the thread.

(%i5) display2d:false;

(%o5) false
(%i6) integrate(x*cos(x^3),x);

(%o6) (gamma_incomplete(2/3,%i*x^3)+gamma_incomplete(2/3,-%i*x^3))/6
(%i7) domain:complex;

(%o7) complex
(%i8) integrate(x*cos(x^3),x);

(%o8) ((sqrt(3)*%i-1)*gamma_incomplete(2/3,%i*x^3)
       +(-sqrt(3)*%i-1)*gamma_incomplete(2/3,-%i*x^3))
       *(x^3)^(1/3)
       /(12*x)

But this doesn't resolve the original issue.

(%i1) display2d:false;

(%o1) false
(%i2) integrate(x*cos(x^3),x,0,1/2);

(%o2) gamma_incomplete(2/3,%i/8)/6+gamma_incomplete(2/3,-%i/8)/6-gamma(2/3)/3
(%i3) domain:complex;

(%o3) complex
(%i4) integrate(x*cos(x^3),x,0,1/2);

(%o4) gamma_incomplete(2/3,%i/8)/6+gamma_incomplete(2/3,-%i/8)/6-gamma(2/3)/3

comment:3 Changed 10 years ago by kcrisman

Robert points out that in 5.27.0 this is at least a little better, though we still need to change domain.

Maxima 5.27.0 http://maxima.sourceforge.net
using Lisp SBCL 1.0.24
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) display2d:false;

(%o1) false
(%i2) domain:complex;

(%o2) complex
(%i3) integrate (x*cos(x^3), x, 0, 1/2); 

(%o3) %i*gamma_incomplete(2/3,%i/8)/(4*sqrt(3))
 -gamma_incomplete(2/3,%i/8)/12-%i*gamma_incomplete(2/3,-%i/8)/(4*sqrt(3))
 -gamma_incomplete(2/3,-%i/8)/12+gamma(2/3)/6
(%i4) expand(float(%));

(%o4) .1247560409610377

comment:4 Changed 10 years ago by kcrisman

  • Cc mjo added
  • Report Upstream changed from N/A to Fixed upstream, in a later stable release.

But since we already use domain:complex in the initialization code for Maxima for integrals and such, presumably an upgrade would fix this?

comment:5 Changed 10 years ago by mjo

See my last comment on #11238 for why we can't have nice things.

That's the only real issue with maxima-5.27. Some messages are getting mangled e.g. when maxima asks us for more information, but they aren't substantial regressions.

comment:6 Changed 9 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:7 Changed 8 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:8 Changed 8 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:9 Changed 8 years ago by pbruin

  • Priority changed from major to trivial
  • Report Upstream changed from Fixed upstream, in a later stable release. to N/A
  • Work issues set to add doctest

This seems to have been fixed some time ago; it works correctly in both Sage 5.13 and 6.2.

sage: integrate(x*cos(x^3),(x,0,1/2)).n()
0.124756040961038

comment:10 Changed 8 years ago by pbruin

  • Authors set to Peter Bruin
  • Branch set to u/pbruin/12947-maxima_integral
  • Commit set to 35f60813c540f3eaabc45200f3e61f3a227ad1d2
  • Status changed from new to needs_review
  • Work issues add doctest deleted

comment:11 Changed 8 years ago by rws

  • Reviewers set to Ralf Stephan

Funny I get:

File "src/sage/interfaces/maxima_lib.py", line 736, in sage.interfaces.maxima_lib.MaximaLib.sr_integral
Failed example:
    integrate(x*cos(x^3),(x,0,1/2)).n()
Expected:
    0.124756040961038
Got:
    0.124756040961038 + 2.77555756156289e-17*I

Please set positive when you think this is resolved.

comment:12 Changed 8 years ago by git

  • Commit changed from 35f60813c540f3eaabc45200f3e61f3a227ad1d2 to 3d9a07cd45fd002020a1f67cb2394af5a38d76a0

Branch pushed to git repo; I updated commit sha1. New commits:

3d9a07cTrac 12947: add tolerance to doctest

comment:13 Changed 8 years ago by pbruin

On the system I tested this on (GNU/Linux x86_64) this gives the expected answer, both with and without the Maxima upgrade of #13973. I just checked on ARM; the first time I got the same answer as you, but then I merged this branch (did not even rebuild) and only got the expected answer from then on. I added # abs tol 1e-16; can you check if this solves the problem for you?

comment:14 Changed 8 years ago by kcrisman

I would argue that the complex bit means that we are not quite ready to merge this if it happens on some platforms. The question is whether this is a Maxima bug or a Sage/Pynac? simplification bug.

comment:15 follow-up: Changed 8 years ago by rws

Would you equally mind if there is a real digit wrong in the 17-th place? Regardless, it's fine on AMD now.

comment:16 Changed 8 years ago by rws

There are 60(!) lines in the sage code with # abs tol and nearly all have lower tolerance (in the sense of 'larger error')

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

comment:17 Changed 8 years ago by pbruin

The complex bit apparently comes from applying Pynac's evalf() to the following expression:

sage: e = integrate(x*cos(x^3),(x,0,1/2)); e
1/12*I*sqrt(3)*gamma(2/3, 1/8*I) - 1/12*I*sqrt(3)*gamma(2/3, -1/8*I) - 1/12*gamma(2/3, 1/8*I) - 1/12*gamma(2/3, -1/8*I) + 1/6*gamma(2/3)
sage: e.n(53)
0.124756040961038 + 2.77555756156289e-17*I
sage: e._convert({'parent': CC})  # called by .n()
0.124756040961038 + 2.77555756156289e-17*I

(This is on ARM; the imaginary part does not appear on x86_64.)

There are known precision problems in the incomplete Gamma function, so this may be related to #7099, although we only want the default 53 bits of precision here.

comment:18 Changed 8 years ago by pbruin

The presence of the imaginary part does not depend on #7099. At least on ARM it is essentially random (maybe it depends on FPU rounding flags or something similar):

$ ./sage -c 'print integrate(x*cos(x^3),(x,0,1/2)).n()'
0.124756040961038
$ ./sage -c 'print integrate(x*cos(x^3),(x,0,1/2)).n()'
0.124756040961038 - 2.77555756156289e-17*I
$ ./sage -c 'print integrate(x*cos(x^3),(x,0,1/2)).n()'
0.124756040961038 - 1.38777878078145e-17*I
$ ./sage -c 'print integrate(x*cos(x^3),(x,0,1/2)).n()'
0.124756040961038 + 1.38777878078145e-17*I

comment:19 in reply to: ↑ 15 Changed 8 years ago by kcrisman

Would you equally mind if there is a real digit wrong in the 17-th place? Regardless, it's fine on AMD now.

Haha! No, but as I said sometimes the complex bit indicates a different kind of problem, or has in the past. Given Peter's experiment in comment:18, though, I guess we can make that a different ticket.

comment:20 Changed 8 years ago by pbruin

The imaginary part turns out to be "simply" the fact that floating-point addition is not associative, and for some reason the terms in the expression e from comment:17 are added in random order by .n().

sage: e = integrate(x*cos(x^3),(x,0,1/2)); e
1/12*I*sqrt(3)*gamma(2/3, 1/8*I) - 1/12*I*sqrt(3)*gamma(2/3, -1/8*I) - 1/12*gamma(2/3, 1/8*I) - 1/12*gamma(2/3, -1/8*I) + 1/6*gamma(2/3)
sage: v = [term.n() for term in e.operands()]; v
[0.0454319516149362 + 0.166098636946833*I,
 0.0454319516149362 - 0.166098636946833*I,
 -0.0958970927532841 + 0.0262301494946934*I,
 -0.0958970927532841 - 0.0262301494946934*I,
 0.225686323237733]
sage: uniq([sum(p.action(v)) for p in permutations(5)])
[0.124756040961038 - 2.77555756156289e-17*I,
 0.124756040961038,
 0.124756040961038 + 2.77555756156289e-17*I,
 0.124756040961038 - 1.38777878078145e-17*I,
 0.124756040961038,
 0.124756040961038 + 1.38777878078145e-17*I,
 0.124756040961038 - 2.77555756156289e-17*I,
 0.124756040961038,
 0.124756040961038 + 2.77555756156289e-17*I,
 0.124756040961038 - 1.38777878078145e-17*I,
 0.124756040961038 + 1.38777878078145e-17*I,
 0.124756040961038 - 2.77555756156289e-17*I,
 0.124756040961038 + 2.77555756156289e-17*I,
 0.124756040961038]

This gives 14 different outcomes depending on the summation order (some of the differences are hidden due to the number of digits printed). I don't think this can be solved easily, so I propose we stick to the # abs tol fix.

comment:21 Changed 8 years ago by git

  • Commit changed from 3d9a07cd45fd002020a1f67cb2394af5a38d76a0 to 8274e4764945e773078535aabe9f262462416e57

Branch pushed to git repo; I updated commit sha1. New commits:

8274e47Merge branch 'develop' into ticket/12947-numerical_integral

comment:22 Changed 8 years ago by rws

  • Status changed from needs_review to positive_review

As already said ...

comment:23 Changed 8 years ago by vbraun

The # abs tol just changes how the floating point numbers are compared, it doesn't know about complex numbers.

sage -t --long src/sage/interfaces/maxima_lib.py
**********************************************************************
File "src/sage/interfaces/maxima_lib.py", line 763, in sage.interfaces.maxima_lib.MaximaLib.sr_integral
Failed example:
    integrate(x*cos(x^3),(x,0,1/2)).n()  # abs tol 1e-16
Expected:
    0.124756040961038
Got:
    0.124756040961038 - 1.38777878078145e-17*I
**********************************************************************

and on a different machine

sage -t --long src/sage/interfaces/maxima_lib.py
**********************************************************************
File "src/sage/interfaces/maxima_lib.py", line 763, in sage.interfaces.maxima_lib.MaximaLib.sr_integral
Failed example:
    integrate(x*cos(x^3),(x,0,1/2)).n()  # abs tol 1e-16
Expected:
    0.124756040961038
Got:
    0.124756040961038 - 2.77555756156289e-17*I
**********************************************************************

comment:24 Changed 8 years ago by vbraun

  • Status changed from positive_review to needs_work

comment:25 Changed 8 years ago by git

  • Commit changed from 8274e4764945e773078535aabe9f262462416e57 to 72706eb6697fba361e138ddaca0765f8415eb8f3

Branch pushed to git repo; I updated commit sha1. New commits:

72706ebTrac 12947: handle numerical noise correctly

comment:26 Changed 8 years ago by pbruin

  • Status changed from needs_work to positive_review

This should solve it; probably no reason to ask for another review.

comment:27 Changed 8 years ago by vbraun

  • Branch changed from u/pbruin/12947-maxima_integral to 72706eb6697fba361e138ddaca0765f8415eb8f3
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.