Sympy计算含Gamma函数导数未化简polygamma，如何获取浮点结果？

问题 #

我编写了一段Python脚本，用于计算包含Gamma函数的函数的导数。但代入数值后，Sympy并未将polygamma(0, 1.05)这类表达式化简为digamma对应的浮点值，而是保留未求解的polygamma形式。以下是我的代码及输出示例，请问如何修改以得到最终浮点结果？

原代码 #

import mpmath
import time
import sympy
x = sympy.symbols ('x')
s = sympy.symbols ('s')
from sympy import S, I, pi, gamma, lambdify

Original = ((((sympy.pi**(x/2))*(s**x))/sympy.gamma((x/2)+1))-(((2*s)/(x**0.5))**x))
Prime = Original.diff (x)
Prime = lambdify ((x, s), Prime, modules = 'sympy')
for s_times_10 in range (1, 31):
    s = float (int (s_times_10) / 10)
    for x_times_10 in range (1, 151):
        x = float ((int (x_times_10) / 10))
        print ("x: " + str (x) + ", s: " + str (s))
        print (Prime (x, s))
        if (x > 0.3):
            if (Prime (x + 0.1, s) < Prime (x, s)):
                print ("MAXIMUM N LOCATED: " + str (x))
                time.sleep (1)
                break
        print ("=======")
        time.sleep (0.5)

输出示例 #

x: 0.1, s: 0.1
-0.579691734344519 - 0.432005861274674*polygamma(0, 1.05)
=======
x: 0.2, s: 0.1
-0.175935858863424 - 0.371829906705536*polygamma(0, 1.1)
=======
x: 0.3, s: 0.1
0.0107518316667914 - 0.31889065255819*polygamma(0, 1.15)
=======
x: 0.4, s: 0.1
0.098684205215577 - 0.27256963654143*polygamma(0, 1.2)
=======
x: 0.5, s: 0.1
0.133891927091406 - 0.232239660951436*polygamma(0, 1.25)
MAXIMUM N LOCATED: 0.5

解决方案 #

问题根源是lambdify指定了modules='sympy'，这会让Sympy始终保持符号表达式形式，不会自动计算浮点数值。可以通过两种方式修改：

方法1：改用mpmath作为计算模块 #

mpmath对特殊函数的数值计算支持更直接，能自动将polygamma(0, z)（即digamma函数）解析为浮点结果。只需修改lambdify的参数：

Prime = lambdify ((x, s), Prime, modules = 'mpmath')

方法2：用Sympy的N()函数强制求值 #

如果坚持使用Sympy模块，在调用Prime(x, s)后用sympy.N()包裹，强制转换为数值：

result = sympy.N(Prime(x, s))
print(result)

修改后的完整代码（以方法1为例） #

import mpmath
import time
import sympy
x = sympy.symbols ('x')
s = sympy.symbols ('s')
from sympy import S, I, pi, gamma, lambdify

Original = ((((sympy.pi**(x/2))*(s**x))/sympy.gamma((x/2)+1))-(((2*s)/(x**0.5))**x))
Prime = Original.diff(x)
# 改用mpmath作为计算模块
Prime = lambdify((x, s), Prime, modules='mpmath')

for s_times_10 in range(1, 31):
    s_val = float(s_times_10 / 10)
    for x_times_10 in range(1, 151):
        x_val = float(x_times_10 / 10)
        print(f"x: {x_val}, s: {s_val}")
        result = Prime(x_val, s_val)
        print(result)
        if x_val > 0.3:
            if Prime(x_val + 0.1, s_val) < result:
                print(f"MAXIMUM N LOCATED: {x_val}")
                time.sleep(1)
                break
        print("=======")
        time.sleep(0.5)

补充说明 #

• 注意循环中的变量名不要和符号变量x、s重复，避免变量混淆（示例中改为x_val、s_val）。
• mpmath计算特殊函数的效率和精度都能满足需求，是数值计算场景下的更优选择。

内容的提问来源于stack exchange，提问作者TH3TA