如何用Python求解非线性方程组？现有SymPy代码无结果求方案

问题描述 #

需要求解以下非线性方程组，目标是找到曲线2y³ + y² - y⁵ - x⁴ + 2x³ - x² = 0上切线斜率为1的点（第二个方程是隐函数求导后得到的dy/dx=1）：

• 曲线方程：2y³ + y² - y⁵ - x⁴ + 2x³ - x² = 0
• 斜率条件：(-4x³ + 6x² - 2x) / (5y⁴ - 6y² - 2y) = 1

使用SymPy的solve函数直接求解未得到结果，但已知存在5个有效解，希望通过Python库完成求解。

解决方案 #

方法1：SymPy符号求解（转化为整式方程） #

原第二个方程为分式形式，直接求解可能受限于符号计算的局限性，建议先转化为整式方程并添加分母不为零的约束，避免引入增根。

修改后的代码：

from sympy import symbols, Eq, solve, simplify

x, y = symbols('x y')
# 定义原曲线方程
eq1 = Eq(2*y**3 + y**2 - y**5 - x**4 + 2*x**3 - x**2, 0)
# 将分式斜率方程转化为整式形式
eq2 = Eq(-4*x**3 + 6*x**2 - 2*x, 5*y**4 - 6*y**2 - 2*y)
# 添加分母不为零的约束条件
constraint = 5*y**4 - 6*y**2 - 2*y != 0

# 求解方程组，返回字典形式的解
solutions = solve([eq1, eq2], [x, y], dict=True)
# 过滤掉违反约束的无效解
valid_solutions = [sol for sol in solutions if constraint.subs(sol)]

# 输出所有有效精确解
for idx, sol in enumerate(valid_solutions, 1):
    print(f"解 {idx}: x = {simplify(sol[x])}, y = {simplify(sol[y])}")

方法2：Scipy数值求解（适合符号解难获取的场景） #

如果符号求解仍无法得到结果，可使用数值求解方法，结合GeoGebra给出的大致解作为初始猜测值，得到近似解：

from scipy.optimize import fsolve
import numpy as np

# 定义方程组，返回两个方程的误差值
def equations(vars):
    x, y = vars
    eq1 = 2*y**3 + y**2 - y**5 - x**4 + 2*x**3 - x**2
    eq2 = (-4*x**3 + 6*x**2 - 2*x)/(5*y**4 - 6*y**2 - 2*y) - 1
    return [eq1, eq2]

# 基于GeoGebra结果设置初始猜测值（可根据实际调整）
initial_guesses = [
    (0, 0),
    (1, 0),
    (0, 1),
    (1, 1),
    (0.5, 0.5)
]

solutions = []
for guess in initial_guesses:
    try:
        # 求解方程组
        sol = fsolve(equations, guess)
        # 验证解的有效性（误差足够小）
        if np.allclose(equations(sol), [0, 0], rtol=1e-6):
            solutions.append(sol)
    except Exception:
        continue

# 去重并输出近似解
unique_solutions = np.unique(np.round(solutions, 6), axis=0)
for idx, (x, y) in enumerate(unique_solutions, 1):
    print(f"解 {idx}: x ≈ {x:.6f}, y ≈ {y:.6f}")

说明 #

• 符号求解适合能获取精确解的场景，转化为整式方程后能降低SymPy的计算难度；
• 数值求解依赖合理的初始猜测值，从GeoGebra获取的大致坐标能有效提升求解成功率；
• 无论哪种方法，都需要过滤掉分母为零的解，避免出现无效的增根。

内容的提问来源于stack exchange，提问作者Sylvain Bérubé