Possible Bug in 'integrate'

[wang270020]

I found a bug in intergrate. I did same integral in Sympy, Mathmatics, and matlb. The result of the former is different from the laters.

sympy

d = 1
T = 0.25
K0, Z = sympy.symbols("K0, Z")
Z0 = sympy.cosh(K0 * (Z + d)) / sympy.cosh(K0 * d)
Y0 = sympy.sqrt(2) / 2
F0 = sympy.integrate(Z0 * Y0, (Z, -d, -T)) * 2 / (d - T)

res1:

-1.33333333333333*sqrt(2)*sinh(0.75*K0)/(K0*cosh(K0))

mathmatics:

2/0.75*Integrate[Cosh[k0*(z + 1)]/Cosh[k0*1]*Sqrt[2]/2, {z, -1, -0.25}]

res2:
(1.88562 Sech[k0] Sinh[0.75 k0])/k0

[oscarbenjamin]

Seems to be a bug in meijerg somewhere:

In [55]: integrate(cosh(y*(x + 1)), (x, -1, -0.25), meijerg=True)
Out[55]: 
-sinh(0.75⋅y) 
──────────────
      y       

In [56]: integrate(cosh(y*(x + 1)), (x, -1, -0.25), meijerg=False)
Out[56]: 
⎧sinh(0.75⋅y)                            
⎪────────────  for y > -∞ ∧ y < ∞ ∧ y ≠ 0
⎨     y                                  
⎪                                        
⎩    0.75              otherwise  

[arnabnandikgp]

I went through the governing code for this problem at seems that meijerint_indefinite is being triggered

sympy/sympy/integrals/meijerint.py

Lines 1653 to 1687 in c81c816

	def meijerint_indefinite(f, x):
	"""
	Compute an indefinite integral of ``f`` by rewriting it as a G function.
	Examples
	========
	>>> from sympy.integrals.meijerint import meijerint_indefinite
	>>> from sympy import sin
	>>> from sympy.abc import x
	>>> meijerint_indefinite(sin(x), x)
	-cos(x)
	"""
	f = sympify(f)
	results = []
	for a in sorted(_find_splitting_points(f, x) | {S.Zero}, key=default_sort_key):
	res = _meijerint_indefinite_1(f.subs(x, x + a), x)
	if not res:
	continue
	res = res.subs(x, x - a)
	if _has(res, hyper, meijerg):
	results.append(res)
	else:
	return res
	if f.has(HyperbolicFunction):
	_debug('Try rewriting hyperbolics in terms of exp.')
	rv = meijerint_indefinite(
	_rewrite_hyperbolics_as_exp(f), x)
	if rv:
	if not isinstance(rv, list):
	from sympy.simplify.radsimp import collect
	return collect(factor_terms(rv), rv.atoms(exp))
	results.extend(rv)
	if results:
	return next(ordered(results))

and seems like the for loop is responsible for the error and if the expression is allowed to to run through the HyperbolicFunction block then we do not get the error, maybe their position maybe an issue since we can always first check for the f to contain HyperbolicFunction and then go for the substitution method since we have explicit condition to handle specifically the Hyperbolic functions(??).
Though when I ran similar examples like

 >>> integrate(sinh(y*(x + 1)), (x, -1, -0.25),meijerg=True)
cosh(0.75*y)/y - 1/y 

and for even

>>> integrate(sinh(y*(x + 1)), x,meijerg=True)
sinh(y*(x + 1))/y

sympy gives correct result going through the loop
looking into the issue deeper seems that the source of error comes from _meijerint_indefinite_1 and more specifically from the hyperexpand's handling of meijerg expressions with dummy variables as what is observed that there is no difference between the parameters that are being passed in case of computing
integrate(sinh(y*(x + 1)), x,meijerg=True) except of the expression z being in terms of dummy variables for the definite integration case.

[arnabnandikgp]

For the given expression i.e cosh(y*(x + 1)) the following code is being run in hyperexpand and it seems to me that

sympy/sympy/simplify/hyperexpand.py

Lines 1983 to 1989 in c81c816

	if premult == 1:
	C = C.applyfunc(make_simp(z0))
	r = reduce(lambda s,m: s+m[0]*m[1], zip(C, f.B.subs(f.z, z0)), S.Zero)*premult
	res = r.subs(z0, z)
	if rewrite:
	res = res.rewrite(rewrite)
	return res

during substitution(in line 1986) in case of the definite integral the dummy variable that is being used has to go through some operations like sqrt(_Dummy**2) and what i think is that since the limits for the integral are in the negative domain sympy computes that to -_Dummy and that is the source of error. Though converting sqrt(x**2) to -x is not wrong mathematically I guess we don't want that to happen here.I checked the integral's behaviour for mixed sign limits and both positive limits and they are correct. I think we can further investigate and set some assumptions straight in the initial steps of the calculation of such definite integral but I fear that might bring in more problems of its own(?).

[oscarbenjamin]

The problem goes away of symbols with assumptions are not used:

diff --git a/sympy/integrals/integrals.py b/sympy/integrals/integrals.py
index d07f84dfa0..c3dfa1cfcc 100644
--- a/sympy/integrals/integrals.py
+++ b/sympy/integrals/integrals.py
@@ -500,7 +500,8 @@ def doit(self, **hints):
             else:
                 d = None
             if d:
-                reps[x] = d
+                #reps[x] = d
+                pass
         if reps:
             undo = {v: k for k, v in reps.items()}
             did = self.xreplace(reps).doit(**hints)

This is the code that replaces the integration variable with the negative assumption:

sympy/sympy/integrals/integrals.py

Lines 479 to 511 in c81c816

	# first make sure any definite limits have integration
	# variables with matching assumptions
	reps = {}
	for xab in self.limits:
	if len(xab) != 3:
	# it makes sense to just make
	# all x real but in practice with the
	# current state of integration...this
	# doesn't work out well
	# x = xab[0]
	# if x not in reps and not x.is_real:
	# reps[x] = Dummy(real=True)
	continue
	x, a, b = xab
	l = (a, b)
	if all(i.is_nonnegative for i in l) and not x.is_nonnegative:
	d = Dummy(positive=True)
	elif all(i.is_nonpositive for i in l) and not x.is_nonpositive:
	d = Dummy(negative=True)
	elif all(i.is_real for i in l) and not x.is_real:
	d = Dummy(real=True)
	else:
	d = None
	if d:
	reps[x] = d
	if reps:
	undo = {v: k for k, v in reps.items()}
	did = self.xreplace(reps).doit(**hints)
	if isinstance(did, tuple): # when separate=True
	did = tuple([i.xreplace(undo) for i in did])
	else:
	did = did.xreplace(undo)
	return did

It is possible that this substitution is not valid for meijerg and that it would actually be better to replace all integration variables with symbols that have no assumptions. Algorithms like meijerg and hyperexpand have their own way of handling branches of things like sqrt(x**2) that maybe conflict with the automatic evaluation seen here.

Another possibility is that this is invalid:

In [5]: -sqrt(pi)*I*meijerg([], [], [S.Half], [0], -x**2/4)
Out[5]: 
              ⎛       │   2 ⎞
      ╭─╮1, 0 ⎜       │ -x  ⎟
-ⅈ⋅√π⋅│╶┐     ⎜       │ ────⎟
      ╰─╯0, 2 ⎝1/2  0 │  4  ⎠

In [6]: hyperexpand(_)
Out[6]: sinh(x)

At least according to wikipedia this formula is only valid for -pi < arg(x) <= 0 i.e. x needs to be positive real or it should be in the lower complex plane. For negative x the formula does not hold but it perhaps would with a minus sign. Simple testing suggests:

In [58]: e1 = -sqrt(pi)*I*meijerg([], [], [S.Half], [0], -x**2/4)

In [59]: e2 = Piecewise((sinh(x), Eq(x, 0) | (-pi<arg(x))&(arg(x) <= 0)), (-sinh(x), True))

In [60]: e2
Out[60]: 
⎧sinh(x)   for (arg(x) ≤ 0 ∧ arg(x) > -π) ∨ x = 0
⎨                                                
⎩-sinh(x)                otherwise               

In [61]: for v in [-1, 0, 1, I, -I, 1+I, 1-I, -1-I, -1-I]:
    ...:     print((e1 - e2).subs(x, v).n())
    ...: 
0.e-125
0
0.e-125
0.e-125*I
0.e-125*I
0.e-128 - 0.e-126*I
-0.e-132 + 0.e-129*I
0.e-128 - 0.e-126*I
0.e-128 - 0.e-126*I

Rewriting in the other direction gives a different G function but it also seems incorrect for negative numbers:

In [62]: from sympy.integrals.meijerint import _rewrite_single

In [63]: _rewrite_single(sinh(x), x)
Out[63]: ([(pi**(3/2), 0, meijerg(((), (1,)), ((1/2,), (0, 1)), x**2/4))], True)

In [66]: a3, _, e3 = _rewrite_single(sinh(x), x)[0][0]

In [67]: (a3*e3) - sinh(x)
Out[67]: 
             ⎛          │  2⎞          
 3/2 ╭─╮1, 0 ⎜      1   │ x ⎟          
π   ⋅│╶┐     ⎜          │ ──⎟ - sinh(x)
     ╰─╯1, 3 ⎝1/2  0, 1 │ 4 ⎠          

In [68]: ((a3*e3) - sinh(x)).subs(x, 1).n()
Out[68]: 0.e-125

In [69]: ((a3*e3) - sinh(x)).subs(x, -1).n()
Out[69]: 2.35040238728760

[arnabnandikgp]

I think that in order for the assumptions to be ignored before operations like .subs ,using the following diff

diff --git a/sympy/simplify/hyperexpand.py b/sympy/simplify/hyperexpand.py
index 58d5e9e0c1..7f1f019ca5 100644
--- a/sympy/simplify/hyperexpand.py
+++ b/sympy/simplify/hyperexpand.py
@@ -1983,7 +1983,11 @@ def carryout_plan(f, ops):
         if premult == 1:
             C = C.applyfunc(make_simp(z0))
         r = reduce(lambda s,m: s+m[0]*m[1], zip(C, f.B.subs(f.z, z0)), S.Zero)*premult
-        res = r.subs(z0, z)
+        reps = {s: Dummy(s.name, positive=True)
+                 for s in z.free_symbols if s.is_positive is not True and s.is_Dummy}
+        _z = z.subs(reps)
+        res = r.subs(z0, _z)
+        res = res.subs({r: s for s, r in reps.items()})
         if rewrite:
             res = res.rewrite(rewrite)
         return res

this gives us the solution we are looking for.

[oscarbenjamin]

Setting the symbols to positive is not correct and will cause problems elsewhere.

[arnabnandikgp]

After going through the meijer g function in a bit more detail myself,now I agree with you that replacing all the dummy variables with positive assumption may not be a good idea

It is possible that this substitution is not valid for meijerg and that it would actually be better to replace all integration variables with symbols that have no assumptions.

I think that maybe going forward with this will do, though I am not sure of how the following is implemented

Algorithms like meijerg and hyperexpand have their own way of handling branches of things like sqrt(x**2)

If possible can you please point me out the code for the same.

[asmeurer]

This looks similar to #25910. I would say as in that PR that we should try to fix the underlying problem rather than swapping out symbols with assumptions in meijerint.

[oscarbenjamin]

we should try to fix the underlying problem rather than swapping out symbols with assumptions in meijerint

We might want to do both. The meijerint code already has things to take care of branches and so on and it is possible that assumptions based evaluation messes with that.

[oscarbenjamin]

Fixed by gh-26173

[oscarbenjamin]

In [2]: import sympy

In [3]: d = 1
   ...: T = 0.25
   ...: K0, Z = sympy.symbols("K0, Z")
   ...: Z0 = sympy.cosh(K0 * (Z + d)) / sympy.cosh(K0 * d)
   ...: Y0 = sympy.sqrt(2) / 2
   ...: F0 = sympy.integrate(Z0 * Y0, (Z, -d, -T)) * 2 / (d - T)

In [4]: F0
Out[4]: 
1.33333333333333⋅√2⋅sinh(0.75⋅K₀)
─────────────────────────────────
           K₀⋅cosh(K₀)           

In [5]: F0.n()
Out[5]: 
1.88561808316413⋅sinh(0.75⋅K₀)
──────────────────────────────
         K₀⋅cosh(K₀)