Using Sympy to integrate a complicated expression. Program does not terminate even after an hour

Question

I am trying to compute a definite integral using sympy. Below is my code:

import sympy as sp
import math

eps, psi = sp.symbols("eps psi", positive=True)
y = 1 / (math.sqrt(8)*math.pi**2)*sp.integrate(1/sp.sqrt(eps - psi) * d2crho_dpsi2, (psi, 0, eps))

I have pre-computed d2crho_dpsi2 and it is a complicated expression in terms of psi. It is as follows:

(110826921669719*(-110826921669719*psi**2/(500000000000000*(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)**2*(-74440218185373*psi/500000000000000 - LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000))**2) - 20*psi/((1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)**3*(-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)))*(-110826921669719*psi**3*log(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)/(25000000000000*(-74440218185373*psi/500000000000000 - LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000))**3) + 74440218185373*psi**2/(2500000000000*(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)*(-74440218185373*psi/500000000000000 - LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000))**2) + 134335984549343*psi/(200000000000*(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)**2*(-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)))/(500000000000000000*(-110826921669719*psi**2*log(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)/(5000000000000000*(-74440218185373*psi/500000000000000 - LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000))**2) + psi/((1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)*(-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)))**2) + (65999841839257*psi**3/(100000000000000*(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)**2*(-74440218185373*psi/500000000000000 - LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000))**3) + 886615373357751*psi**2/(100000000000000*(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)**3*(-74440218185373*psi/500000000000000 - LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000))**2) + 600*psi/((1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)**4*(-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)))/(-74440218185373*psi**2*log(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)/(50000000000000*(-74440218185373*psi/500000000000000 - LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000))**2) + 134335984549343*psi/(2000000000000*(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)*(-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000))))/(-74440218185373*psi**2*log(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)/(50000000000000*(-74440218185373*psi/500000000000000 - LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000))**2) + 134335984549343*psi/(2000000000000*(1 + (-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)/psi)*(-psi - 134335984549343*LambertW(-74440218185373*psi*exp(-74440218185373*psi/500000000000000)/500000000000000)/20000000000000)))

I ran the code for one hour but still did not get an output. Is the integrand too complicated for sympy to integrate, or worse still, the resulting integral may not have a closed form solution?

Answer 1

It's very unlikely to find a closed form for integrals as complicated as the one shown unless you can simplify them somehow.

Your expression has a repeating subexpression:

In [87]: d2crho_dpsi2.atoms(LambertW)
Out[87]: 
⎧ ⎛                   -74440218185373⋅ψ  ⎞⎫
⎪ ⎜                   ────────────────── ⎟⎪
⎪ ⎜                    500000000000000   ⎟⎪
⎨ ⎜-74440218185373⋅ψ⋅ℯ                   ⎟⎬
⎪W⎜──────────────────────────────────────⎟⎪
⎪ ⎝           500000000000000            ⎠⎪
⎩                                         ⎭

This is of the form W(a*exp(a)) which is just equal to a for many values of a (it depends which branch of the W function you are on and the value of a). With that replacement the expression simplifies to zero:

In [88]: a = Wild('a')

In [89]: factor(d2crho_dpsi2).replace(LambertW(a*exp(a)), a)
Out[89]: 0

That replacement is only valid for

In [101]: solve(-74440218185373*psi/500000000000000 > -1).n()
Out[101]: ψ < 6.71679922746716

However beyond that point the same replacement would apply if you were thinking about the W_{-1} branch. Depending on which branch you intended it is possible that the correct answer here is just that your expression is equal to zero (hence easy to integrate).

See: https://en.wikipedia.org/wiki/Lambert_W_function#Inverse