使用SymPy nsolve求解非线性方程组遇数值奇异矩阵错误求助

错误含义解释 #

ZeroDivisionError: matrix is numerically singular 表示nsolve使用的牛顿迭代法中，当前迭代点对应的雅可比矩阵是数值奇异的——即矩阵行列式为0（或无限接近0），无法计算逆矩阵。而牛顿法需要通过逆雅可比矩阵更新迭代值，因此触发除零类错误。这类问题通常由方程组冗余约束、符号与数值计算混用、初始猜测严重偏离真实解、变量间线性依赖等原因导致。

针对你的代码的解决方法 #

1. 替换numpy函数为SymPy原生函数 #

代码中使用np.absolute会导致符号计算与数值计算混用，破坏SymPy的符号推导逻辑，进而导致雅可比矩阵计算异常，需替换为SymPy的sym.Abs。

2. 修正初始猜测值，匹配已知约束 #

根据参数定义Cpz=2、Cz=0，Cpvertdist和Cvertdist的真实值必然为2，但你的初始猜测值为1.6和2.6，严重偏离合理范围，直接导致迭代初期雅可比矩阵奇异。需调整初始值，使其贴合方程约束。

3. 简化冗余变量与方程（可选，提升稳定性） #

方程组中eq3和eq6是完全相同的约束（均等于abs(Cpz-Cz)=2），可直接将Cpvertdist和Cvertdist替换为常量2，减少未知数数量，避免变量依赖引发的奇异问题。

修正后的代码 #

import sympy as sym

# 定义参数值
Cpz = 2
intensity = 0.0499717154066123
length = 10
Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz, Dx, Dy, Dz = 0, length, 0, length, length, 0, 0, 0, 0, length, 0, 0
ni, nm = 1.1, 1

# 定义符号变量，添加正实数约束缩小解空间
thetaC1 = sym.symbols("thetaC1", real=True, positive=True)
thetaC2 = sym.symbols("thetaC2", real=True, positive=True)
Cpx = sym.symbols("Cpx", real=True, positive=True)
Px = sym.symbols("Px", real=True, positive=True)
Pz = sym.symbols("Pz", real=True, positive=True)
Cphorizdist = sym.symbols("Cphorizdist", real=True, positive=True)
Cpdist3D = sym.symbols("Cpdist3D", real=True, positive=True)
Chorizdist = sym.symbols("Chorizdist", real=True, positive=True)
Cdist3D = sym.symbols("Cdist3D", real=True, positive=True)

# 利用已知约束简化方程，删除冗余项
eq1 = sym.Eq(sym.asin((nm/ni)*sym.sin(thetaC2)), thetaC1)
eq2 = sym.Eq(sym.Abs(Px - Cpx), Cphorizdist)
eq4 = sym.Eq(sym.sqrt(Cphorizdist**2 + 2**2), Cpdist3D)
eq5 = sym.Eq(sym.Abs(Cpx - Cx), Chorizdist)
eq7 = sym.Eq(sym.sqrt(Chorizdist**2 + 2**2), Cdist3D)
eq8 = sym.Eq(Cphorizdist / Cpdist3D, sym.sin(thetaC1))
eq9 = sym.Eq(Cpx / Cdist3D, sym.sin(thetaC2))
eq10 = sym.Eq(Chorizdist / 2, sym.tan(thetaC2))
eq11 = sym.Eq(1 / ((Cpdist3D + Cdist3D)**2), intensity)

# 修正初始猜测值，贴合约束逻辑
initial_guesses = [
    0.5,    # thetaC1
    0.6,    # thetaC2
    2.0,    # Cpx
    1.0,    # Cphorizdist
    2.2,    # Cpdist3D
    2.0,    # Chorizdist
    2.2,    # Cdist3D
    3.0,    # Px
    2.0     # Pz
]

# 确保未知数、方程、初始猜测数量一致
result = sym.nsolve(
    (eq1, eq2, eq4, eq5, eq7, eq8, eq9, eq10, eq11),
    (thetaC1, thetaC2, Cpx, Cphorizdist, Cpdist3D, Chorizdist, Cdist3D, Px, Pz),
    initial_guesses,
    verify=False,
    maxsteps=50000,
    method='lm'  # 列文伯格-马夸尔特法，对奇异矩阵鲁棒性更强
)

print("求解结果：")
print(f"thetaC1: {result[0]}")
print(f"thetaC2: {result[1]}")
print(f"Cpx: {result[2]}")
print(f"Cphorizdist: {result[3]}")
print(f"Cpdist3D: {result[4]}")
print(f"Chorizdist: {result[5]}")
print(f"Cdist3D: {result[6]}")
print(f"Px: {result[7]}")
print(f"Pz: {result[8]}")

额外建议 #

• 给符号变量添加positive=True约束，帮助SymPy缩小解空间，提升迭代稳定性。
• 若仍出现奇异问题，可通过method='lm'指定列文伯格-马夸尔特法，该方法对奇异矩阵的鲁棒性远高于默认的牛顿法。

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