Python: Integrate product of exponential and gaussian function

Question

I have a function y which is a product of y1 = exp(-x) and exp(-x^2):

y1 = exp(-x)

y2 = exp(-x**2)

y = y1*y2 = exp(-x)*exp(-x**2) = exp(-x **2-x)

Integrating the functions y1 or y2 works fine using sympy:

>>> import sympy as sy
>>> import numpy as np
>>> x = sy.Symbol('x')
>>> y1 = sy.exp(-x)
>>> sy.integrate(y1, (x, 0, np.inf))
1

and

>>> y2 = sy.exp(-x**2)
>>> sy.integrate(y2, (x, 0, np.inf))
0.5*sqrt(pi)

However, whenever I'm trying to integrate the product y1*y2 = y, the integral is not accepted:

>>> y = y1*y2
>>> sy.integrate(y, (x, 0, np.inf))
Integral(exp(-x)*exp(-x**2), (x, 0, np.inf))

Perhaps I'm blind and missing something obvious.

Answer 1

If sympy cannot evaluate it symbolically, then it just returns an Integral object. You need to either use a more powerful symbolic calculator like maxima or Wolfram|Alpha (link to solution)), or else resort to numerical integration. Here are a few options for numerically integrating it:

import sympy as sy
import numpy as np
import scipy as sp
import scipy.integrate

x = sy.Symbol('x')
y1 = sy.exp(-x)
y2 = sy.exp(-x**2)
y = y1*y2

# Pure scipy without sympy
print(sp.integrate.quad(lambda x: np.exp(-x) * np.exp(-(x**2)), 0, np.inf))

# Mix of scipy and sympy
print(sp.integrate.quad(lambda x_arg: y.subs(x, x_arg).evalf(), 0, np.inf))

# Without using scipy
print(sy.mpmath.quad(sy.lambdify([x], y), [0, np.inf]))