Python使用SymPy求解高次幂单未知数方程的报错与性能问题

使用SymPy的solveset方法求解相关方程时遇到如下报错：

mpmath.libmp.libhyper.NoConvergence: convergence to root failed; try n < 15 or maxsteps > 50

原始求解代码如下：

from __future__ import division
from sympy import *

n = Symbol('n', real=True, positive=True)
d = Symbol('d', real=True, positive=True)
i = Symbol('i', real=True, positive=True)
o = Symbol('o', real=True, positive=True)
z = Symbol('z', real=True, positive=True)
eq_left = Symbol('eq_left', real=True, positive=True)
eq_right = Symbol('eq_right', real=True, positive=True)

d = 100
i = 0.0001
o = 0.0001

eq_left = (1 + d*(z/n))**n
eq_right = (1 + d*(z/(n+1)))**(n+1)*(1-d*i)*(1-d*o)

for every in range(1, 41):
  n = every
  results = solveset(Eq(eq_left, eq_right), z, domain=S.Reals)
  print(results)

代码逻辑为枚举1到40的整数n值，求解每个n对应的唯一未知量z。实际运行时n取值接近12就会触发前述报错，且幂次越高计算耗时越长，预估n=40时单步求解时长可达1小时。需要解决两个核心问题：

• 如何规避收敛失败报错，正确求解这类含高次幂的方程
• 如何优化求解性能，将批量计算的耗时压缩到可接受范围

问题根源 #

solveset默认优先尝试求符号解析解，遇到高次指数项时会调用mpmath的高精度全局求根逻辑：一方面高次幂计算容易出现数值量级溢出、精度丢失，触发收敛失败；另一方面符号推导的开销随幂次呈指数级上升，求解速度极慢。

具体优化方法 #

• 替换求解逻辑：单未知数固定参数的方程不需要求符号解，直接用数值迭代求根，稳定性和速度远高于符号求解。
• 方程等价变形：由于z>0时方程两边所有项均为正，对等式两边同时取自然对数，将高次幂运算转化为乘法运算，大幅降低计算时的数值量级，从根源减少收敛失败概率，变形过程不会引入增根或失根。
• 限定搜索范围：根据方程的复利类计算属性，z的取值范围极小，给求根算法传入明确的搜索区间或合理初始值，避免全局搜索的冗余开销，同时防止迭代跑飞。

高性能实现方案（推荐） #

直接用SciPy的数值求根接口，总耗时不到0.1秒即可算完所有n对应的z值，无收敛问题：

import math
from scipy.optimize import root_scalar

# 固定常量提前计算，避免循环内重复运算
d = 100
i = 0.0001
o = 0.0001
const_factor = (1 - d*i) * (1 - d*o)

def get_z(n_val):
    # 变形为f(z)=0的求根形式，取对数消除高次幂
    def eq(z):
        left = n_val * math.log(1 + d * z / n_val)
        right = (n_val + 1) * math.log(1 + d * z / (n_val + 1)) + math.log(const_factor)
        return left - right
    # 布伦特法求根，传入明确的搜索区间，稳定性强速度快
    res = root_scalar(eq, bracket=[1e-8, 0.1], method='brentq')
    return res.root

# 批量计算输出
for n in range(1, 41):
    print(f"n={n}, z={get_z(n):.8f}")

纯SymPy实现方案 #

如果不想引入SciPy依赖，可将solveset替换为专用于数值求根的nsolve，调整mpmath参数后全量计算仅需数秒：

from sympy import Symbol, log, nsolve
import mpmath

# 调整数值计算参数，调高最大迭代步数
mpmath.mp.dps = 15
mpmath.mp.maxsteps = 200

z = Symbol('z', real=True, positive=True)
d = 100
i = 0.0001
o = 0.0001
const_factor = (1 - d*i) * (1 - d*o)

for n_val in range(1, 41):
    # 提前做对数变形
    eq = n_val * log(1 + d*z/n_val) - (n_val+1)*log(1 + d*z/(n_val+1)) - log(const_factor)
    # 传入初始猜测值，避免全局搜索
    z_res = nsolve(eq, z, 0.01)
    print(f"n={n_val}, z={float(z_res):.8f}")

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