Opened 10 years ago

Closed 8 years ago

# Bug in integrating x*cos(x^3)

Reported by: Owned by: kcrisman burcin trivial sage-6.3 calculus mjo Peter Bruin Ralf Stephan N/A 72706eb 72706eb6697fba361e138ddaca0765f8415eb8f3

### Description

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

### 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),(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
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

• 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
• 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:

 ​3d9a07c `Trac 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: ↓ 19 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 interval')

Version 1, edited 8 years ago by rws (previous) (next) (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:

 ​8274e47 `Merge branch 'develop' into ticket/12947-numerical_integral`

### comment:22 Changed 8 years ago by rws

• Status changed from needs_review to positive_review

### 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:

 ​72706eb `Trac 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.