Sympy solve求解符号方程组仅得2个正确值的问题排查

问题背景 #

使用SymPy的solve函数求解8个gamma变量的符号方程组时，仅gamma6和gamma7的结果符合预期，其余6个变量计算值与目标值偏差极大。

原代码 #

import numpy as np
from sympy import symbols, exp, Eq, solve

gamma0, gamma1, gamma2, gamma3, gamma4, gamma5, gamma6, gamma7 = symbols('gamma0 gamma1 gamma2 gamma3 gamma4 gamma5 gamma6 gamma7')
gamma = [gamma0, gamma1, gamma2, gamma3, gamma4, gamma5, gamma6, gamma7]
x = symbols('x')
a = 20
b = 5

u =[
    [0.388518*gamma0*exp(-2.012490*x) - 0.675612*gamma1*exp(-0.369901*x) + 0.336576*gamma2*exp(0.369901*x) - 0.103503*gamma3*exp(2.012490*x) + 20.168718, 
    0.163057*gamma4*exp(-2.421891*x) - 0.951601*gamma5*exp(-0.951784*x) + 0.086427*gamma6*exp(0.754538*x) - 0.045759*gamma7*exp(2.528248*x) + 2.394974],

    [-0.900008*gamma0*exp(-2.012490*x) - 0.520883*gamma1*exp(-0.369901*x) + 0.398678*gamma2*exp(0.369901*x) - 0.168303*gamma3*exp(2.012490*x) + 20.034468, 
    -0.976653*gamma4*exp(-2.421891*x) - 0.238418*gamma5*exp(-0.951784*x) + 0.122769*gamma6*exp(0.754538*x) - 0.080425*gamma7*exp(2.528248*x) + 2.329339],

    [-0.168303*gamma0*exp(-2.012490*x) - 0.398678*gamma1*exp(-0.369901*x) + 0.520883*gamma2*exp(0.369901*x) - 0.900008*gamma3*exp(2.012490*x) + 19.932869, 
    -0.095726*gamma4*exp(-2.421891*x) - 0.120587*gamma5*exp(-0.951784*x) + 0.326236*gamma6*exp(0.754538*x) - 0.963388*gamma7*exp(2.528248*x) + 2.504797], 

    [-0.103503*gamma0*exp(-2.012490*x) - 0.336576*gamma1*exp(-0.369901*x) + 0.675612*gamma2*exp(0.369901*x) + 0.388518*gamma3*exp(2.012490*x) + 19.892514, 
    -0.101960*gamma4*exp(-2.421891*x) - 0.151889*gamma5*exp(-0.951784*x) + 0.933288*gamma6*exp(0.754538*x) + 0.251633*gamma7*exp(2.528248*x) + 3.335936]
    ]

# Define the symbolic equations
eqns = [
    Eq(u[0][0].subs(x, 0), 0),
    Eq(u[1][0].subs(x, 0), 0),
    Eq(u[0][0].subs(x, b), u[0][1].subs(x, b)),
    Eq(u[1][0].subs(x, b), u[1][1].subs(x, b)),
    Eq(u[2][0].subs(x, b), u[2][1].subs(x, b)),
    Eq(u[3][0].subs(x, b), u[3][1].subs(x, b)),
    Eq(u[2][1].subs(x, a), 0),
    Eq(u[3][1].subs(x, a), 0)
]

# Solve the equations for the gamma variables
gamma_symbols = [gamma0, gamma1, gamma2, gamma3, gamma4, gamma5, gamma6, gamma7]
solutions = solve(eqns, gamma_symbols)

# Print the solutions
for symbol, value in solutions.items():
    print(f"{symbol}: {value}")

预期结果 #

gamma0: 3.3451862507144963721760546900288
gamma1: 30.086443529744653795656405805539
gamma2: -3.3917670835986188504875542362514
gamma3: 0.00013878499955117735347592615338105
gamma4: -775727.3319427781475689057795114
gamma5: -941.09415775171256674990468706789
gamma6: -0.0000010944407599058234322536821517657
gamma7: 0.00000000000000000000013958849298737787410748516872778

实际运行结果 #

gamma0: -2953067.95180723
gamma1: -1383474.12048191      
gamma2: -13623.4036144577      
gamma3: 1.03375979503954       
gamma4: -8212497765.78320      
gamma5: -15039124.0279307      
gamma6: -0.00000109854435279209
gamma7: 1.40162266686094E-22  

问题分析与解决方法 #

核心原因 #

这个方程组本质是线性方程组，但其中包含如exp(2.528248*20)这类超大指数项，导致系数矩阵的条件数极高（病态矩阵）。SymPy的solve默认符号求解，面对低精度浮点数构成的病态系统时，误差会被急剧放大，最终得到偏离预期的结果。

可行解决方案 #

1. 转为数值线性方程组求解（推荐）

将方程组整理为Ax = b的矩阵形式，使用SciPy的最小二乘法求解，该方法对病态矩阵的数值误差有更好的鲁棒性：

import numpy as np
from scipy.linalg import lstsq

# 构建系数矩阵A和右端向量b
A = []
b = []

# 方程1: u[0][0](x=0) = 0
A.append([0.388518, -0.675612, 0.336576, -0.103503, 0, 0, 0, 0])
b.append(-20.168718)

# 方程2: u[1][0](x=0) = 0
A.append([-0.900008, -0.520883, 0.398678, -0.168303, 0, 0, 0, 0])
b.append(-20.034468)

# 计算x=5处的指数值
exp1 = np.exp(-2.012490*5)
exp2 = np.exp(-0.369901*5)
exp3 = np.exp(0.369901*5)
exp4 = np.exp(2.012490*5)
exp5 = np.exp(-2.421891*5)
exp6 = np.exp(-0.951784*5)
exp7 = np.exp(0.754538*5)
exp8 = np.exp(2.528248*5)

# 方程3: u[0][0](x=5) - u[0][1](x=5) = 0
A.append([0.388518*exp1, -0.675612*exp2, 0.336576*exp3, -0.103503*exp4,
          -0.163057*exp5, 0.951601*exp6, -0.086427*exp7, 0.045759*exp8])
b.append(2.394974 - 20.168718)

# 方程4: u[1][0](x=5) - u[1][1](x=5) = 0
A.append([-0.900008*exp1, -0.520883*exp2, 0.398678*exp3, -0.168303*exp4,
          0.976653*exp5, 0.238418*exp6, -0.122769*exp7, 0.080425*exp8])
b.append(2.329339 - 20.034468)

# 方程5: u[2][0](x=5) - u[2][1](x=5) = 0
A.append([-0.168303*exp1, -0.398678*exp2, 0.520883*exp3, -0.900008*exp4,
          0.095726*exp5, 0.120587*exp6, -0.326236*exp7, 0.963388*exp8])
b.append(2.504797 - 19.932869)

# 方程6: u[3][0](x=5) - u[3][1](x=5) = 0
A.append([-0.103503*exp1, -0.336576*exp2, 0.675612*exp3, 0.388518*exp4,
          0.101960*exp5, 0.151889*exp6, -0.933288*exp7, -0.251633*exp8])
b.append(3.335936 - 19.892514)

# 计算x=20处的指数值
exp5_a = np.exp(-2.421891*20)
exp6_a = np.exp(-0.951784*20)
exp7_a = np.exp(0.754538*20)
exp8_a = np.exp(2.528248*20)

# 方程7: u[2][1](x=20) = 0
A.append([0,0,0,0, -0.095726*exp5_a, -0.120587*exp6_a, 0.326236*exp7_a, -0.963388*exp8_a])
b.append(-2.504797)

# 方程8: u[3][1](x=20) = 0
A.append([0,0,0,0, -0.101960*exp5_a, -0.151889*exp6_a, 0.933288*exp7_a, 0.251633*exp8_a])
b.append(-3.335936)

# 转换为数值数组
A = np.array(A, dtype=np.float64)
b = np.array(b, dtype=np.float64)

# 最小二乘法求解
gamma, _, _, _ = lstsq(A, b)

# 输出结果
gamma_names = ['gamma0','gamma1','gamma2','gamma3','gamma4','gamma5','gamma6','gamma7']
for name, val in zip(gamma_names, gamma):
    print(f"{name}: {val}")

2. 提高SymPy求解精度

将所有浮点数替换为SymPy的高精度浮点数（如Float(0.388518, precision=50)），减少低精度带来的误差累积，再用solve求解。

3. 拆分方程组

观察方程组结构，先解x=0处的方程得到gamma0-gamma3的关系，再解x=20处的方程得到gamma4-gamma7的关系，最后代入x=5处的衔接条件求解，降低系统耦合度。

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