Search code examples
pythonsympysymbolic-math

Numerical integration vs symbolic integration in python

I want to look at a plot of a function intensity(r) versus space, so r (radially symmetric).

However, I derive my intensity from intensity(r) = integrate(integrand(r), (x,0,5)), where integrand = exp(-x**2) * exp(np.pi*1j*(-x)) * besselj(0, r*x) * x.

All the syntax above is using the sympy package, so I first define x,y = symbols('x r') .

I was using symbolic variables because it seemed to be making visually easier, leaving r as a variable until the end, when I plot it and I assign it a numerical value.

However doing that horrible integral with a symbolic variable seems to be taking ages.

  • Is there any way of performing numerical integration with symbolic variables?

  • Is the only other alternative defining the values of r a priori and finding the integral for each one of them?

Solution

  • Aside: when you create a symbolic expression, keep it symbolic. Don't mix in a real float np.pi and a complex float 1j, use SymPy's pi and I instead. 

    from sympy import exp, pi, I, besselj, symbols
x, r = symbols('x r')
integrand = exp(-x**2) * exp(pi*I*(-x)) * besselj(0, r*x) * x

    But yes, it doesn't look like SymPy can integrate the product of a Bessel function with exp(-x**2) * exp(pi*I*(-x)). This happens already with r replaced by 1, so the symbolic nature of r is beside the point.

    To directly answer your question:

    Is there any way of performing numerical integration with symbolic variables?

    There is not, like there is no dry water. It's a contradiction in terms.

    Is the only other alternative defining the values of r a priori and finding the integral for each one of them?

    Yes. It can be done through SymPy (which will call mpmath):

    >>> intensity = lambda r_: Integral(integrand.subs(r, r_), (x, 0, 5)).evalf()
>>> intensity(3)
0.0783849036516177 - 0.125648626220306*I

    It's not quite clear how you plan to plot this function, given that it's complex-valued. Perhaps you meant to plot the absolute value of intensity?

    Anyway, integration with SymPy/mpmath (pure Python) is too slow for plotting. You'd be better off with SciPy's quad for integration. It does not handle complex integrands, so I integrate the real and complex parts separately.

    from scipy.integrate import quad
from scipy.special import jn
integrand = lambda x, r: np.exp(-x**2) * np.exp(np.pi*1j*(-x)) * jn(0, r*x) * x
intensity = lambda r: np.sqrt(quad(lambda x: np.real(integrand(x, r)), 0, 5)[0]**2 + quad(lambda x: np.imag(integrand(x, r)), 0, 5)[0]**2)

    Now intensity(3) is evaluated a lot faster than in the previous version. And we can plot it:

    import matplotlib.pyplot as plt
t = np.linspace(0, 3)
plt.plot(t, np.vectorize(intensity)(t))

    plot