Symbolic integration is done through the following steps:

1. The definite integral of a polynomial function on a finite domain [a,b] can be computed very accurately via the fundamental theorem of calculus, using the numpy.polynomial module. For instance, to calculate the integral of the polynomial p(x)=x5 on the interval [-1,1].
2. We could issue the following part integrating both parts.
3. You need to do the following:

In [1]: import numpy as np
In [2]: p = np.poly1d([1,0,0,0,0,0]); 
   ...: print (p)
   ...: print (p.integ())
   5
1 x
       6
0.1667 x
In [3]: p.integ()(1.0) - p.integ()(-1.0)
Out[3]: 0.0

5. In general, obtaining exact values for a definite integral of a generic function is hard and computationally inefficient. This is possible in some cases through symbolic integration with the aid of the Risch algorithm (for elementary functions) and Meijer G-functions (for non-elementary integrals). Both methods can be called with the common routine, integrate, in the SymPy library. The routine is clever enough to decide which algorithm to use, depending on the source function.
6. Let's show you the following example starting with the definite integral of the polynomial from the previous case:

In [4]: from sympy import integrate, symbols
In [5]: x, y = symbols('x y', real=True)
In [6]: integrate(x**5, x)
Out[6]: x**6/6
In [7]: integrate(x**5, (x, -1, 1))
Out[7]: 0