Wrong integration result involving square root of absolute value

[rodluger]

The function

sin(x / sqrt(abs(x))

is odd, so its definite integral over the range (-1, 1) should be identically zero (see WolframAlpha result here).

Sympy yields a different, nonzero value:

>>> Integral(sin(x / sqrt(abs(x))), (x, -1, 1)).doit()  
-4*cos(1) + 4*sin(1)

>>> Integral(sin(x / sqrt(abs(x))), (x, -1, 1)).doit().evalf()
1.20467471575903

The issue is related to the divergence of the numerator as x -> 0, and possibly the fact that the absolute value is not differentiable at that point. Note that sympy has no issues when the sqrt(abs(x)) is in the numerator.

I'm using sympy 1.7.1.

[oscarbenjamin]

The integral can be checked numerically with .n():

In [5]: Iwhole = Integral(sin(x/sqrt(abs(x))), (x, -1, 1))

In [6]: Iwhole
Out[6]: 
1                 
⌠                 
⎮     ⎛   x   ⎞   
⎮  sin⎜───────⎟ dx
⎮     ⎜  _____⎟   
⎮     ⎝╲╱ │x│ ⎠   
⌡                 
-1                

In [7]: Iwhole.n()
Out[7]: 0.e-100

In [8]: Iwhole.doit()
Out[8]: -4⋅cos(1) + 4⋅sin(1)

In [9]: Iwhole.doit().n()
Out[9]: 1.20467471575903

It seems that the positive part is fine but the negative part is wrong:

In [11]: Ipos = Integral(sin(x/sqrt(abs(x))), (x, 0, 1))

In [12]: Ineg = Integral(sin(x/sqrt(abs(x))), (x, -1, 0))

In [13]: Ipos.n()
Out[13]: 0.602337357879514

In [14]: Ipos.doit().n()
Out[14]: 0.602337357879514

In [15]: Ineg.n()
Out[15]: -0.602337357879514

In [16]: Ineg.doit().n()
Out[16]: 0.602337357879514

Internally integrate splits this into two integrals over the positive and negative parts to eliminate the abs so the wrong result comes from this antiderivative from heurisch:

In [17]: integrate(sin(x/sqrt(-x)), x, heurisch=True)
Out[17]: -2⋅ⅈ⋅√x⋅cosh(√x) + 2⋅ⅈ⋅sinh(√x)

That goes wrong in the second line here:

sympy/sympy/integrals/heurisch.py

Lines 745 to 746 in aa22709

	antideriv = solution.subs(rev_mapping)
	antideriv = cancel(antideriv).expand(force=True)

What we have there is:

(Pdb) l
743  	
744  	    if solution is not None:
745  	        antideriv = solution.subs(rev_mapping)
746  	        antideriv = cancel(antideriv).expand(force=True)
747  	
748  ->	        if antideriv.is_Add:
749  	            antideriv = antideriv.as_independent(x)[1]
750  	
751  	        return indep*antideriv
752  	    else:
753  	        if retries >= 0:
(Pdb) p solution.subs(rev_mapping)
(2*x**2*cos(x/sqrt(-x)) + 2*(-x)**(3/2)*sin(x/sqrt(-x)))/(x*sqrt(-x))
(Pdb) p antideriv
-2*I*sqrt(x)*cosh(sqrt(x)) + 2*I*sinh(sqrt(x))

The first result for solution is correct but after expand it is wrong:

In [5]: e1 = (2*x**2*cos(x/sqrt(-x)) + 2*(-x)**(S(3)/2)*sin(x/sqrt(-x)))/(x*sqrt(-x))

In [6]: e1.diff(x).cancel()
Out[6]: 
   ⎛  x   ⎞
sin⎜──────⎟
   ⎜  ____⎟
   ⎝╲╱ -x ⎠

In [7]: e2 = -2*I*sqrt(x)*cosh(sqrt(x)) + 2*I*sinh(sqrt(x))

In [8]: e2.diff(x).cancel()
Out[8]: -ⅈ⋅sinh(√x)

In [10]: e1.subs(x, -1).n()
Out[10]: 0.602337357879514

In [11]: e2.subs(x, -1).n()
Out[11]: -0.602337357879514

What we have there is:

In [19]: e1
Out[19]: 
   2    ⎛  x   ⎞         3/2    ⎛  x   ⎞
2⋅x ⋅cos⎜──────⎟ + 2⋅(-x)   ⋅sin⎜──────⎟
        ⎜  ____⎟                ⎜  ____⎟
        ⎝╲╱ -x ⎠                ⎝╲╱ -x ⎠
────────────────────────────────────────
                    ____                
                x⋅╲╱ -x                 

In [20]: e1.subs(x, -1).n()
Out[20]: 0.602337357879514

In [21]: cancel(e1) == e1
Out[21]: True

In [22]: e1.expand()
Out[22]: 
       ⎛  x   ⎞                
2⋅x⋅cos⎜──────⎟                
       ⎜  ____⎟                
       ⎝╲╱ -x ⎠        ⎛  x   ⎞
─────────────── - 2⋅sin⎜──────⎟
       ____            ⎜  ____⎟
     ╲╱ -x             ⎝╲╱ -x ⎠

In [23]: e1.expand().subs(x, -1).n()
Out[23]: 0.602337357879514

In [24]: e1.expand(force=True)
Out[24]: -2⋅ⅈ⋅√x⋅cosh(√x) + 2⋅ⅈ⋅sinh(√x)

In [25]: e1.expand(force=True).subs(x, -1).n()
Out[25]: -0.602337357879514

So it seems that using expand with force=True is leading to an incorrect "simplification" of the output result.

The basic difference is:

In [30]: sin(x/sqrt(-x))
Out[30]: 
   ⎛  x   ⎞
sin⎜──────⎟
   ⎜  ____⎟
   ⎝╲╱ -x ⎠

In [31]: sin(x/sqrt(-x)).expand(force=True)
Out[31]: -ⅈ⋅sinh(√x)

In [32]: sin(x/sqrt(-x)).subs(x, -1)
Out[32]: -sin(1)

In [33]: sin(x/sqrt(-x)).expand(force=True).subs(x, -1)
Out[33]: sin(1)

The docstring for expand says

sympy/sympy/core/function.py

Lines 2537 to 2538 in aa22709

	If the ``force`` hint is used, assumptions about variables will be ignored
	in making the expansion.

The reason that this matters is:

In [1]: e3 = x/sqrt(-x)

In [2]: e3
Out[2]: 
  x   
──────
  ____
╲╱ -x 

In [3]: e3.expand()
Out[3]: 
  x   
──────
  ____
╲╱ -x 

In [4]: e3.expand(force=True)
Out[4]: -ⅈ⋅√x

This happens in the power_base hint.

I think that is the documented and expected behaviour for expand_power_base. The force=True parameter is a convenience for when you don't care about always having a fully correct result e.g.:

In [6]: expand_power_base((x*y)**z)
Out[6]: 
     z
(x⋅y) 

In [7]: expand_power_base((x*y)**z, force=True)
Out[7]: 
 z  z
x ⋅y 

The problem then is the fact that heurisch should not be using force=True.

With this diff:

diff --git a/sympy/integrals/heurisch.py b/sympy/integrals/heurisch.py
index c08decb763..6fdc8bb0d4 100644
--- a/sympy/integrals/heurisch.py
+++ b/sympy/integrals/heurisch.py
@@ -743,7 +743,7 @@ def find_non_syms(expr):
 
     if solution is not None:
         antideriv = solution.subs(rev_mapping)
-        antideriv = cancel(antideriv).expand(force=True)
+        antideriv = cancel(antideriv).expand()
 
         if antideriv.is_Add:
             antideriv = antideriv.as_independent(x)[1]

we get the expected result:

In [1]: integrate(sin(x / sqrt(abs(x))), (x, -1, 1))
Out[1]: 0

[oscarbenjamin]

That leads to some test failures e.g.:

In [18]: integrate(sin(log(x**2))) == x*sin(2*log(x))/5 - 2*x*cos(2*log(x))/5
Out[18]: False

However that has improved the result I think to make it correct for negative x:

In [19]: I = Integral(sin(log(x**2)), (x, -1, 0))

In [20]: I
Out[20]: 
0                 
⌠                 
⎮     ⎛   ⎛ 2⎞⎞   
⎮  sin⎝log⎝x ⎠⎠ dx
⌡                 
-1                

In [21]: I.n()
Out[21]: -0.400000000000000

In [22]: I.doit()
Out[22]: -2/5

On master we instead have:

In [1]: I = Integral(sin(log(x**2)), (x, -1, 0))

In [2]: I.n()
Out[2]: -0.400000000000000

In [3]: I.doit()
Out[3]: 
  2⋅cosh(2⋅π)   ⅈ⋅sinh(2⋅π)
- ─────────── + ───────────
       5             5     

In [4]: _.n()
Out[4]: -107.098704593499 + 53.5489788082033⋅ⅈ

The other tests failures are for:

integrate(sqrt(x**2/((y - x)*(y + x))), x)
Integral(li(y*x**2), x).doit()
integrate(1/(x**2 + y), x)

Those should be investigated similarly for negative x or y to see if the answers on master are correct. If they are not fully valid then the tests can just be changed and the diff shown above can be applied along with the OP example as a regression test.

I'm marking this as easy to fix if anyone wants to do those things.

[Jatin2020-24]

That leads to some test failures e.g.:

In [18]: integrate(sin(log(x**2))) == x*sin(2*log(x))/5 - 2*x*cos(2*log(x))/5
Out[18]: False

However that has improved the result I think to make it correct for negative x:

In [19]: I = Integral(sin(log(x**2)), (x, -1, 0))

In [20]: I
Out[20]: 
0                 
⌠                 
⎮     ⎛   ⎛ 2⎞⎞   
⎮  sin⎝log⎝x ⎠⎠ dx
⌡                 
-1                

In [21]: I.n()
Out[21]: -0.400000000000000

In [22]: I.doit()
Out[22]: -2/5

On master we instead have:

In [1]: I = Integral(sin(log(x**2)), (x, -1, 0))

In [2]: I.n()
Out[2]: -0.400000000000000

In [3]: I.doit()
Out[3]: 
  2⋅cosh(2⋅π)   ⅈ⋅sinh(2⋅π)
- ─────────── + ───────────
       5             5     

In [4]: _.n()
Out[4]: -107.098704593499 + 53.5489788082033⋅ⅈ

The other tests failures are for:

integrate(sqrt(x**2/((y - x)*(y + x))), x)
Integral(li(y*x**2), x).doit()
integrate(1/(x**2 + y), x)

Those should be investigated similarly for negative x or y to see if the answers on master are correct. If they are not fully valid then the tests can just be changed and the diff shown above can be applied along with the OP example as a regression test.

I'm marking this as easy to fix if anyone wants to do those things.

For function integrate(sqrt(x**2/((y - x)*(y + x))),(x,-1,0))
On Master, we get: (see WolframAlpha result here)

In[1]: e = integrate(sqrt(x**2/((y - x)*(y + x))),(x,-1,0))

Out[1]: 
       ____           ________        ________
 2     ╱ 1      2     ╱   1           ╱   1    
y ⋅   ╱  ──  - y ⋅   ╱  ──────  +    ╱  ────── 
     ╱    2         ╱    2          ╱    2     
   ╲╱    y        ╲╱    y  - 1    ╲╱    y  - 1 

After improving results we get: (eliminating force=True from sympy/sympy/integrals/heurisch.py)

In[1]: e = integrate(sqrt(x**2/((y - x)*(y + x))),(x,-1,0))

Out[1]: 
       ____           ________        ________
 2     ╱ 1      2     ╱   1           ╱   1    
y ⋅   ╱  ──  - y ⋅   ╱  ──────  +    ╱  ────── 
     ╱    2         ╱    2          ╱    2     
   ╲╱    y        ╲╱    y  - 1    ╲╱    y  - 1 

We get same results for function integrate(sqrt(x**2/((y - x)*(y + x))),(x,-1,0)) on master and after improvements.

For function integrate(1/(x**2 + y), (x,-1,0))
On master, we have:

In[1]: e = integrate(1/(x**2 + y), (x,-1,0))

Out[1]: 
    _____    ⎛       _____⎞       _____    ⎛      _____⎞       _____    ⎛   
     ╱ -1      ⎜      ╱ -1  ⎟      ╱ -1      ⎜     ╱ -1  ⎟      ╱ -1      ⎜   
    ╱  ─── ⋅log⎜-y⋅  ╱  ─── ⎟     ╱  ─── ⋅log⎜y⋅  ╱  ─── ⎟     ╱  ─── ⋅log⎜- y
  ╲╱    y      ⎝   ╲╱    y  ⎠   ╲╱    y      ⎝  ╲╱    y  ⎠   ╲╱    y      ⎝   
- ─────────────────────────── + ────────────────────────── + ─────────────────
               2                            2                               2 

     _____    ⎞       _____    ⎛      _____    ⎞
    ╱ -1      ⎟      ╱ -1      ⎜     ╱ -1      ⎟
⋅  ╱  ───  - 1⎟     ╱  ─── ⋅log⎜y⋅  ╱  ───  - 1⎟
 ╲╱    y      ⎠   ╲╱    y      ⎝  ╲╱    y      ⎠
─────────────── - ──────────────────────────────
                                2               

After improving results we get:

In[1]: e = integrate(1/(x**2 + y), (x,-1,0))

Out[1]: 
      _____    ⎛       _____⎞       _____    ⎛      _____⎞       _____    ⎛   
     ╱ -1      ⎜      ╱ -1  ⎟      ╱ -1      ⎜     ╱ -1  ⎟      ╱ -1      ⎜   
    ╱  ─── ⋅log⎜-y⋅  ╱  ─── ⎟     ╱  ─── ⋅log⎜y⋅  ╱  ─── ⎟     ╱  ─── ⋅log⎜- y
  ╲╱    y      ⎝   ╲╱    y  ⎠   ╲╱    y      ⎝  ╲╱    y  ⎠   ╲╱    y      ⎝   
- ─────────────────────────── + ────────────────────────── + ─────────────────
               2                            2                               2 

     _____    ⎞       _____    ⎛      _____    ⎞
    ╱ -1      ⎟      ╱ -1      ⎜     ╱ -1      ⎟
⋅  ╱  ───  - 1⎟     ╱  ─── ⋅log⎜y⋅  ╱  ───  - 1⎟
 ╲╱    y      ⎠   ╲╱    y      ⎝  ╲╱    y      ⎠
─────────────── - ──────────────────────────────
                                2               

We get same results for function integrate(1/(x**2 + y), (x,-1,0)) on master and after improvements.

For Function integrate(1/(x**2 + y), (x,-1,0))
On master, we have:

In[1]: e = integrate(1/(x**2 + y), (x,-1,0))
In[2]: e
Out[2]: 0             
⌠             
⎮    ⎛ 2  ⎞   
⎮  li⎝x ⋅y⎠ dx
⌡             
-1            

In[3]: e.doit()
Out[3]:
 
⎧          ⎛3⋅log(y)        ⎞                            
⎪        Ei⎜──────── + 3⋅ⅈ⋅π⎟                            
⎪          ⎝   2            ⎠                            
⎨li(y) - ────────────────────  for y > -∞ ∧ y < ∞ ∧ y ≠ 0
⎪                 √y                                     
⎪                                                        
⎩             0                        otherwise         

After improvement, we got:

In[1]: e = integrate(1/(x**2 + y), (x,-1,0))
In[2]: e
Out[2]: 
0             
⌠             
⎮    ⎛ 2  ⎞   
⎮  li⎝x ⋅y⎠ dx
⌡             
-1  
In[3]: e.doit()
Out[3]:

⎧          ⎛3⋅log(y)⎞                            
⎪        Ei⎜────────⎟                            
⎪          ⎝   2    ⎠                            
⎨li(y) - ────────────  for y > -∞ ∧ y < ∞ ∧ y ≠ 0
⎪             √y                                 
⎪                                                
⎩         0                    otherwise         

We get different results on master and after improvements for function integrate(1/(x**2 + y), (x,-1,0))

[oscarbenjamin]

The question is whether those changed results are more correct or not. If you substitute different values for y (especially negative or non-real values) and numerically evaluate the integrals can it be shown that the current master output is incorrect for some cases and that the new output is an improvement?

For example the first result on master seems incorrect for some negative values of y e.g.:

In [15]: e = Integral(1/(x**2 + y), (x,-1,0))

In [16]: e.doit().subs(y, -S.Half).n()
Out[16]: -1.24645048028046 + 2.22144146907918⋅ⅈ

It seems that numerical evaluation fails for this integral though:

In [17]: e.subs(y, -S.Half).n()
Out[17]: 0.e+0

Does the suggested change lead to a result that is correct for y=-1/2?

Since heurisch is just being used to compute the antiderivative here maybe it's better to change the limits so that they don't include the singularity e.g.:

In [59]: e = Integral(1/(x**2 + y), (x,-S(1)/10,0))

In [60]: e.doit().subs(y, -S.Half).n()
Out[60]: -0.201349565519501 + 0.e-20⋅ⅈ

In [61]: e.subs(y, -S.Half).n()
Out[61]: -0.201349565519501