使用scipy fsolve求解含对数方程组失败问题排查

问题原因 #

• 初始值严重偏离可行域：你给的x[2]初始值是43，但要求x[2]是10的负几次方量级，fsolve迭代时容易出现常数/x[2]趋近于0或者x[2]变为非正数的情况，直接触发对数函数的定义域错误（log仅接受正数输入）。
• fsolve无内置约束支持：该函数会尝试任意实数作为迭代值，不会遵守你要求的“变量非零正数、x[1]∈[1,2]、x[2]为10负几次方”的约束，导致计算过程中出现非法输入。

解决方案 #

1. 修正初始值 #

直接把初始值调整到符合约束的范围内，让fsolve从可行域附近开始迭代：

import numpy as np
from scipy.optimize import fsolve

def func(x):
    return [x[0] * 0.00096 + x[1]*0.0259*np.log(0.00096/x[2]) - 1.2505,
            x[0] * 0.00301 + x[1]*0.0259*np.log(0.00301/x[2]) - 1.2829,
            x[0] * 0.000077 + x[1]*0.0259*np.log(0.000077/x[2]) - 1.1789,]

# 调整初始值：x[2]改为1e-5（10的负5次方），x[1]取1.5（在1-2之间）
root = fsolve(func, [66, 1.5, 1e-5])
print(root)

2. 变量替换规避定义域问题 #

将x[2]替换为对数形式，确保x[2]始终为正数：

import numpy as np
from scipy.optimize import fsolve

def func_transformed(y):
    # y[0]=x0, y[1]=x1, y[2]=log10(x2)，则x2=10^y[2]
    x2 = 10 ** y[2]
    return [y[0] * 0.00096 + y[1]*0.0259*np.log(0.00096/x2) - 1.2505,
            y[0] * 0.00301 + y[1]*0.0259*np.log(0.00301/x2) - 1.2829,
            y[0] * 0.000077 + y[1]*0.0259*np.log(0.000077/x2) - 1.1789,]

# y[2]初始值对应x2=1e-5，即log10(1e-5)=-5
root_transformed = fsolve(func_transformed, [66, 1.5, -5])
# 还原x2
x_solution = [root_transformed[0], root_transformed[1], 10**root_transformed[2]]
print(x_solution)

3. 使用带约束的优化器 #

如果需要严格遵守变量范围约束，可改用scipy.optimize.minimize将方程组转化为最小二乘问题，并添加约束：

import numpy as np
from scipy.optimize import minimize

def residuals(x):
    # 计算残差的平方和
    res = [x[0] * 0.00096 + x[1]*0.0259*np.log(0.00096/x[2]) - 1.2505,
           x[0] * 0.00301 + x[1]*0.0259*np.log(0.00301/x[2]) - 1.2829,
           x[0] * 0.000077 + x[1]*0.0259*np.log(0.000077/x[2]) - 1.1789,]
    return np.sum(np.square(res))

# 约束条件：x>0，x[1]∈[1,2]
constraints = [{'type': 'ineq', 'fun': lambda x: x[0]},
               {'type': 'ineq', 'fun': lambda x: x[1] - 1},
               {'type': 'ineq', 'fun': lambda x: 2 - x[1]},
               {'type': 'ineq', 'fun': lambda x: x[2]}]

initial_guess = [66, 1.5, 1e-5]
result = minimize(residuals, initial_guess, constraints=constraints)
print(result.x)

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