Sympy strange integral output

Question

I'm trying to solve an integral with sympy. But it gives me a wrong solution. Why?

import sympy
from sympy import Integral, exp, oo

x, y = sympy.symbols("x y", real=True)
b, u, l, t = sympy.symbols("b u l t ", real=True, positive=True)
Fortet = Integral(exp(-l * t) * (sympy.sqrt(2 * sympy.pi * t)) ** (-1) * exp(-((b - u * t - y) ** 2) / (2 * t)),
                  (t, 0, oo))
Fortet.doit()

Result (wrong):

Piecewise((-(-b/2 + y)*sqrt(2*l +
u**2)*(-sqrt(pi)*sinh(sqrt(2)*sqrt(b)*sqrt(l +
u**2/2)*sqrt(polar_lift(1 + y**2/(b*polar_lift(b -
2*y))))*sqrt(polar_lift(b - 2*y))) +
sqrt(pi)*cosh(sqrt(2)*sqrt(b)*sqrt(l + u**2/2)*sqrt(polar_lift(1 +
y**2/(b*polar_lift(b - 2*y))))*sqrt(polar_lift(b - 2*y))))*exp(b*u -
u*y)/(sqrt(pi)*(b - 2*y)*(l + u**2/2)), Abs(arg(1 +
y**2/(b*polar_lift(b - 2*y))) + arg(b - 2*y)) <= pi/2),
(Integral(sqrt(2)*exp(-l*t)*exp(-(b - t*u -
y)**2/(2*t))/(2*sqrt(pi)*sqrt(t)), (t, 0, oo)), True))

Expected (correct) solution:

Solution = (exp((-u)*(b - y)) * exp(sympy.sqrt(u**2 + 2*l)*(b-y)))/(sympy.sqrt(2*l + u**2)) #RIGHT solution

I made an edit. Try to run to code on Python. It really gives back an awkward function. — Prettymath77, May 20 '20 at 11:39
@Prettymath77 Always make sure the code you post is copy-pastable, with imports included. Also, this might be a sympy bug, so best post it there. https://github.com/sympy/sympy/issues/ — Nico Schlömer, May 20 '20 at 14:03

Chris du Plessis · Accepted Answer · 2020-06-03T15:46:26.243

Both results are in fact the same. The first one is probably slightly more correct. You tend see these polar_lift functions whenever SymPy tries to do something like square rooting something when it does not know the signs of the things inside (after integrating)

A polar_lift does not appear below, but this basic Gaussian example shows that SymPy tries to be as general as possible:

from sympy import *
x = Symbol("x", real=True)
y = Symbol("y", real=True)
s = Symbol("s", real=True)  # , positive=True

gaussian = exp(-((x-y)**2)/(2*(s**2)))
nfactor = simplify(integrate(gaussian, (x,-oo,oo)))
print(nfactor)

You need s to be declared as positive: s = Symbol("s", real=True, positive=True). A similar thing happens with these kinds of polar_lift(b - 2*y) functions in your example. It also happens with the question I reference below.

I have no idea why, but N(polar_lift(x)) for any float or int x gives x again; yet, SymPy does not simplify nicely with symbolic x. Turns out if you keep on simplifying, you get nicer and nicer looking answers. I couldn't find anything about polar_lift related to pure math so I don't know what it actually does.

Remember for the simple example above how it gave a piece-wise? Same thing here. So we just take the first piece since the second piece is an un-evaluated integral.

In the code below, I use this question to remove the piece-wise function and then I simplify twice. And finally, I manually remove the polar_lift.

import sympy as sp

x, y = sp.symbols("x y", real=True)
b, u, l, t = sp.symbols("b u l t ", real=True, positive=True)

Fortet = sp.integrate(sp.exp(-l * t) * (sp.sqrt(2 * sp.pi * t)) ** (-1) *
                      sp.exp(-((b - u * t - y) ** 2) / (2 * t)),
                      (t, 0, sp.oo), conds='none')
incorrect = Fortet.simplify().simplify()

correct = eval(str(incorrect).replace("polar_lift", ""))
correct = correct.factor()
print(correct)

The result is:

exp(b*u)*exp(-u*y)*exp(-sqrt(2*l + u**2)*sqrt(b**2 - 2*b*y + y**2))/sqrt(2*l + u**2)

That is close enough to your expression. I couldn't make SymPy simplify the sqrt(b**2 - 2*b*y + y**2) to Abs(b-y) no matter how hard I tried.

Note that either SymPy is still wrong or you are wrong since the powers in the numerator are opposite. So I checked on the Desmos for a numeric answer (top one is yours):

Sympy strange integral output

1 Answers1