如何将Sympy返回的含LambertW函数的超越方程根转为显式形式？

我用Sympy求解某特征方程时得到四个特征值解，其中auto2以LambertW函数形式表示，无法呈现为显式初等函数形式，其余三个解均为显式表达式。希望将该含LambertW函数的解转换为类似的闭形式，请问该如何实现？

我的代码 #

# Modules
import sympy

# Parameters: real and positive values
beta = sympy.Symbol('beta', positive = True, real = True)
epsilon = sympy.Symbol('epsilon', positive = True, real = True)
delta = sympy.Symbol('delta', positive = True, real = True)
sigma = sympy.Symbol('sigma', positive = True, real = True)
gamma = sympy.Symbol('gamma', positive = True, real = True)
mu = sympy.Symbol('mu', positive = True, real = True)
psi = sympy.Symbol('psi', positive = True, real = True)
N = sympy.Symbol('N', positive = True, real = True)
B = sympy.Symbol('B', positive = True, real = True)
L = sympy.Symbol('L', positive = True, real = True)
tau = sympy.Symbol('tau', positive = True, real = True)


# Variables: all are real values >= 0
S = sympy.Symbol('S', negative = False, real = True)
V = sympy.Symbol('V', negative = False, real = True)
I = sympy.Symbol('I', negative = False, real = True)
R = sympy.Symbol('R', negative = False, real = True)


# Equilibria
# Dictionary :: output order to variables :: I, R, S, V
sol1, sol2, sol3 = sympy.solve([B + epsilon*V + delta*R - psi*S - mu*S - beta*S*I/N, psi*S - epsilon*V - mu*V - sigma*beta*V*I/N, beta*S*I/N + sigma*beta*V*I/N - mu*I - gamma*I, gamma*I - mu*R - delta*R], S, V, I, R, dict=True)

# Equilibria point trivial (sol1)::
sol_trivial = list(sol1.values())
It, Rt, St, Vt = sol_trivial

# Linear Part
k = beta*St/N + sigma*beta*Vt/N - (mu + gamma)

A = sympy.Matrix([[-mu, epsilon, beta*St/N, delta], [0, -epsilon-mu, -sigma*beta*Vt/N, 0], [0, 0, k, 0], [0, 0, gamma, -mu-delta]])
B = sympy.Matrix([[-sigma, 0, 0, 0], [sigma, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]])
Id = sympy.Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
Delta = (A + L*Id + B*sympy.exp(-L*tau))
det_Delta = Delta.det()

# Eigenvalues (are the solutions (L variable) to characterist equation "det_Delta")
auto1, auto2, auto3, auto4 = sympy.solve(det_Delta, L)

求解得到的特征值 #

auto1 = mu
auto2 = (tau*(epsilon + mu) + LambertW(sigma*tau*exp(-tau*(epsilon + mu))))/tau
auto3 = (-B*beta*epsilon - B*beta*mu - B*beta*psi*sigma + N*epsilon*gamma*mu + N*epsilon*mu**2 + N*gamma*mu**2 + N*gamma*mu*psi + N*mu**3 + N*mu**2*psi)/(N*mu*(epsilon + mu + psi))
auto4 = delta + mu

已知存在正特征值，该平衡点不稳定，仅auto2因含LambertW函数无法以显式初等函数形式呈现。

解决方案 #

1. 明确LambertW函数的定位：LambertW函数本身就是一种标准的闭形式特殊函数，和指数、对数函数一样属于超越函数范畴。你的问题中因为包含exp(-L*tau)延迟项，属于延迟微分方程的特征方程，这类方程的解通常只能用LambertW函数表示，无法进一步化简为初等函数。

2. 尝试Sympy自动化简：可以调用Sympy的simplify函数尝试自动处理，但大概率不会改变表达式形式，因为Sympy已经识别到这是最简闭形式：

   auto2_simplified = sympy.simplify(auto2)
   

3. 生成近似显式表达式：如果需要初等函数形式的近似解，可以针对特定参数范围做级数展开：

   • 当tau*(epsilon + mu)很大时，LambertW函数有渐近展开：LambertW(x) ≈ ln(x) - ln(ln(x))，代入auto2可得近似式；
   • 用Sympy生成泰勒级数展开，比如在tau趋近于0时展开到指定阶数：

     auto2_series = sympy.series(auto2, tau, 0, 3)  # 展开到tau的2次项
     

4. 符号替换简化写法：如果只是希望表达式更简洁，可以定义新符号替换重复出现的子式，比如：

   C = sympy.Symbol('C')
   auto2_simple = auto2.subs(tau*(epsilon + mu), C)
   

   替换后auto2_simple = (C + LambertW(sigma*tau*exp(-C)))/tau，形式更清爽，但本质仍是LambertW函数的表达。

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