使用Sympy整合复杂的表达式。即使一个小时后程序也不会终止

Question

我正在尝试使用

sympy

计算定积分。下面是我的代码：

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))

我已经预先计算了

d2crho_dpsi2

，就

psi

而言，这是一个复杂的表达式。如下：

(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)))

我运行了代码一小时，但仍然没有得到输出。被积函数是否太复杂而无法对

sympy

进行积分，或者更糟糕的是，所得积分可能没有封闭形式的解？

Answer 1

除非您能以某种方式简化它们，否则不太可能找到像所示复杂的积分的封闭形式。

您的表达式有一个重复的子表达式：

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

其形式为

W(a*exp(a))

，对于

a

的许多值而言，它仅等于

a

（这取决于您所在的 W 函数的哪个分支以及

a

的值）。通过这种替换，表达式简化为零：

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

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

该更换仅适用于

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

但是，超出这一点，如果您正在考虑

W_{-1}

分支，则同样的替换也适用。根据您想要的分支，这里的正确答案可能只是您的表达式等于零（因此易于积分）。

参见： https://en.wikipedia.org/wiki/Lambert_W_function#Inverse