如何正确使用SymPy求解方程？关于nsolve初始猜测的疑问

先修正方程表达式 #

你的原始方程存在括号缺失问题，导致表达式和实际要计算的折现现金流（NPV=0）逻辑不符。比如cf[0]/1+x实际应为cf[0]/(1+x)，最后一项的(cf[4]+(cf[5]/x))/1+x**5应为(cf[4] + cf[5]/x)/(1+x)**5（对应第5期的现金流加永续年金的现值折现）。修正后的代码：

from sympy import symbols, Eq, solve, nsolve, re, im

index = 3586

cf = [202.31091967441085,
 215.90288600761673,
 230.4080089272309,
 245.88763753689724,
 262.4072425909973,
 272.4574399822326]

x = symbols('x')
# 修正后的NPV=0方程
eq1 = Eq(
    -index + 
    cf[0]/(1+x) + 
    cf[1]/(1+x)**2 + 
    cf[2]/(1+x)**3 + 
    cf[3]/(1+x)**4 + 
    (cf[4] + cf[5]/x)/(1+x)**5,
    0
)

方法一：用solve筛选有效解 #

solve会返回所有代数解（包括复数和负实数），可根据业务场景（折现率为正实数）筛选：

sol = solve(eq1)
# 筛选虚部为0、实部大于0的解
valid_sols = [s for s in sol if im(s) == 0 and re(s) > 0]
print(valid_sols)
# 输出：[0.104472587998023]

方法二：用nsolve稳定求解 #

nsolve是数值迭代法，会收敛到离初始猜测最近的根。要避免收敛到错误解，推荐指定搜索区间而非单点猜测：

# 指定x在(0,1)区间内搜索（折现率通常在0到1之间）
res = nsolve(eq1, x, (0, 1))
print(res)
# 输出：0.104472587998023

如果需要单点猜测，可先通过绘图确定根的大致位置：

from sympy import plot
# 绘制x在0到0.2区间的函数图像，观察和x轴交点
plot(eq1.lhs, (x, 0, 0.2))

从图像能看到根在0.1附近，设置初始猜测为0.1左右即可稳定收敛到正确解。

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