欠定方程组求解：Sympy报错问题及随机近似解需求

针对变量数远多于约束方程的情况，核心思路是固定大部分自由变量为接近0的小正数，将问题转化为2方程2变量的定解问题，再用数值方法求解，具体步骤和代码如下：

核心思路 #

由于系统欠定，我们可以选择保留2个变量作为待求解变量，将其余所有自由变量设为接近0的随机小正数（满足正数约束），这样原方程组就简化为定解系统，可通过数值求解器得到可行解。

具体实现步骤 #

1. 生成自由变量的小正数初始值 #

为满足变量为正数的约束，给自由变量生成1e-7~1e-6量级的随机正数，既接近0，又避免后续计算中出现分母为0的问题。

2. 代入自由变量，简化方程组 #

将自由变量的取值代入原方程，得到仅含2个待求解变量的新方程组，整理为残差等于0的形式（即目标值 - 表达式 = 0）。

3. 用数值求解器求解定解方程组 #

使用scipy.optimize.root代替sympy的nsolve，它更适合处理纯数值的非线性方程组，且对初始值的兼容性更好。

代码示例 #

import sympy as sp
import numpy as np
from scipy.optimize import root

# ---------------------- 替换为你的实际参数 ----------------------
# 定义常量（c0、c20为负数，其余非负）
c0 = -2
c1 = 5
c2 = 1
c3 = 2
c4 = 3  # x2的系数
c5 = 1  # x3的系数
# ... 补充其余常量

c20 = -1
c21 = 4
c22 = 1
c23 = 2  # x2的系数
c24 = 0.5  # x3的系数
# ... 补充其余c2x常量

# 定义目标值（即eq1、eq2的左边值）
target_eq1 = 1.2
target_eq2 = 0.6

# 定义变量：x0到x10共11个变量（可根据实际数量调整）
x = sp.symbols('x0:11', positive=True)
x0, x1 = x[0], x[1]
free_vars = x[2:]  # 自由变量x2~x10
# -------------------------------------------------------------

# 生成自由变量的小正数随机值
free_vals = np.random.uniform(1e-7, 1e-6, size=len(free_vars))
subs_dict = {var: val for var, val in zip(free_vars, free_vals)}

# 定义符号形式的残差方程（整理为目标值 - 表达式 = 0）
# 构建eq1的分母：c2*x0 + c3*x1 + 其他自由变量项
denominator_eq1 = c2*x0 + c3*x1 + sum([c4*x[2], c5*x[3]])  # 补充其余自由变量的系数项
residual_eq1 = target_eq1 - (c0 + c1 / denominator_eq1)

# 构建eq2的分母：c21*x0 + c22*x1 + 其他自由变量项
denominator_eq2 = c21*x0 + c22*x1 + sum([c23*x[2], c24*x[3]])  # 补充其余自由变量的系数项
residual_eq2 = target_eq2 - (c20 + c21*x0 / denominator_eq2)

# 将符号方程转为数值函数
residual_func1 = sp.lambdify((x0, x1), residual_eq1.subs(subs_dict), 'numpy')
residual_func2 = sp.lambdify((x0, x1), residual_eq2.subs(subs_dict), 'numpy')

# 组合为scipy所需的向量残差函数
def residual(vars):
    x0_val, x1_val = vars
    return [
        residual_func1(x0_val, x1_val),
        residual_func2(x0_val, x1_val)
    ]

# 设置x0、x1的初始猜测值（正数即可）
initial_guess = [1.0, 1.0]
# 调用求解器，使用LM方法适合非线性方程组
result = root(residual, initial_guess, method='lm')

# 输出结果
if result.success:
    x0_sol, x1_sol = result.x
    if x0_sol > 1e-8 and x1_sol > 1e-8:  # 确保解为正数
        full_sol = subs_dict.copy()
        full_sol[x0] = x0_sol
        full_sol[x1] = x1_sol
        print("可行解：")
        for var in sorted(full_sol.keys(), key=lambda v: str(v)):
            print(f"{var}: {full_sol[var]:.8f}")
    else:
        print("求解得到非正解，尝试调整初始值或自由变量取值")
else:
    print(f"求解失败：{result.message}，可重新生成自由变量值或调整初始猜测")

注意事项 #

• 若求解失败，可尝试：
  • 重新生成自由变量的随机值（调整uniform的范围）
  • 修改x0、x1的初始猜测值（比如设为0.5或2.0）
  • 检查常量和方程形式是否正确，确保分母不会出现0
• 所有变量的正数约束通过初始值设置和结果校验保证，避免出现非正解

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