Search code examples
pythonsympyintegralnumerical-integration

Sympy strange integral output

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
Solution

  • 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): Desmos comparison