SymPy中含Abs的方程组求解出现Invalid NaN comparison错误

首先修正代码笔误 #

你定义的函数变量名为f1，但后续调用时使用了未定义的f，这是代码中的低级错误，先修正为：

f = c * (1/2) * Abs(2*x - 1)

方法1：用solveset替代solve #

SymPy的solve在处理含绝对值的方程组时，可能因分支逻辑处理不当触发NaN比较错误。改用solveset并指定实数域，能更稳定地求解：

from sympy import *
c = 2
x = Symbol('x', real=True)
f = c * (1/2) * Abs(2*x - 1)

x0 = Symbol('x0', real=True)
x1 = Symbol('x1', real=True)

ec1 = Eq(f.subs(x, x0), x1)
ec2 = Eq(f.subs(x, x1), x0)

# 指定实数域求解
sol = solveset((ec1, ec2), (x0, x1), domain=S.Reals)
print(sol)

输出结果为满足条件的实数解集合，比如{(1/3, 1/3), (1, 1)}。

方法2：手动拆分绝对值分支 #

绝对值函数Abs(2x-1)可拆分为两种分段情况，我们组合所有可能的分支（共4种），分别求解线性方程组并验证解的合法性：

from sympy import *
c = 2
x0 = Symbol('x0', real=True)
x1 = Symbol('x1', real=True)

# 情况1：x0 ≥ 1/2 且 x1 ≥ 1/2
eq1_case1 = Eq(c*(1/2)*(2*x0 - 1), x1)
eq2_case1 = Eq(c*(1/2)*(2*x1 - 1), x0)
sol_case1 = solve((eq1_case1, eq2_case1), (x0, x1))
valid_case1 = {(s[x0], s[x1]) for s in sol_case1 if s[x0] >= 1/2 and s[x1] >= 1/2}

# 情况2：x0 ≥ 1/2 且 x1 < 1/2
eq1_case2 = Eq(c*(1/2)*(2*x0 - 1), x1)
eq2_case2 = Eq(c*(1/2)*(1 - 2*x1), x0)
sol_case2 = solve((eq1_case2, eq2_case2), (x0, x1))
valid_case2 = {(s[x0], s[x1]) for s in sol_case2 if s[x0] >= 1/2 and s[x1] < 1/2}

# 情况3：x0 < 1/2 且 x1 ≥ 1/2
eq1_case3 = Eq(c*(1/2)*(1 - 2*x0), x1)
eq2_case3 = Eq(c*(1/2)*(2*x1 - 1), x0)
sol_case3 = solve((eq1_case3, eq2_case3), (x0, x1))
valid_case3 = {(s[x0], s[x1]) for s in sol_case3 if s[x0] < 1/2 and s[x1] >= 1/2}

# 情况4：x0 < 1/2 且 x1 < 1/2
eq1_case4 = Eq(c*(1/2)*(1 - 2*x0), x1)
eq2_case4 = Eq(c*(1/2)*(1 - 2*x1), x0)
sol_case4 = solve((eq1_case4, eq2_case4), (x0, x1))
valid_case4 = {(s[x0], s[x1]) for s in sol_case4 if s[x0] < 1/2 and s[x1] < 1/2}

# 合并所有有效解
all_solutions = valid_case1.union(valid_case2).union(valid_case3).union(valid_case4)
print(all_solutions)

这种方法逻辑清晰，能完全控制绝对值的分支处理，避免SymPy自动处理时的异常。

方法3：数值求解（SciPy） #

若只需近似数值解，可使用SciPy的fsolve函数，无需符号计算：

from scipy.optimize import fsolve
import numpy as np

def equations(vars):
    x0, x1 = vars
    f = lambda x: 2*(1/2)*np.abs(2*x - 1)
    return [f(x0) - x1, f(x1) - x0]

# 基于不同初始猜测求解
sol1 = fsolve(equations, (0, 0))
sol2 = fsolve(equations, (1, 1))
print("数值解1:", sol1)
print("数值解2:", sol2)

数值求解需提供初始猜测值，适合不需要精确符号解的场景。

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