使用least_squares求解方程组遇ValueError：fun需返回至多一维类数组

错误原因 #

你遇到的错误来自scipy.optimize.least_squares调用，核心问题有两点：

• least_squares要求目标函数返回一维数组，但你的equations函数返回了两个长度为100的数组（对应pH_x的100个点），最终形成(2,100)的二维结构，不符合函数要求。
• 逻辑错误：你试图用一次least_squares求解所有pH值对应的nh4和nio2，但该函数仅用于求解单个参数向量，无法批量处理多组独立方程。每个pH值对应一组独立的二元方程组，需要逐个求解。

修正方案 #

1. 重构方程组函数（针对单个pH值） #

将equations改为接收单个pH值和待求解参数的函数，返回两个方程的残差（标量）：

def equations(p, pH):
    nh4, nio2 = p
    h = 10 ** (-pH)
    f = 10 ** (logk - 2 * pH)
    ni2pfree = f/(1 + f / ni_total)
    nh3 = k1 * nh4 / h
    nin4p2 = k2 * (nh4 ** 4) / (h ** 2)
    nin6p2 = k3 * (nh4 ** 6) / (h ** 4)
    return [
        n_total - nh3 - 4*nin4p2 -6*nin6p2 - nh4,
        ni_total - ni2pfree - nin4p2 - nin6p2 - nio2
    ]

2. 遍历所有pH值求解参数 #

用scipy.optimize.fsolve（更适合二元方程组）遍历每个pH值，求解对应参数并保存为数组：

nh4_arr = np.zeros_like(pH_x)
nio2_arr = np.zeros_like(pH_x)

for i, pH in enumerate(pH_x):
    guess = (0.1, 0.1)
    res = scipy.optimize.fsolve(equations, guess, args=(pH,))
    nh4_arr[i] = res[0]
    nio2_arr[i] = res[1]

3. 修正曲线交点求解逻辑 #

原difference函数用绝对值会导致fsolve找极小值而非零点，改成直接返回差值即可：

def difference(pH):
    return interp1(pH) - interp2(pH)

完整修正代码 #

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate, scipy.optimize

k1 = 10 ** (-9.25)
k2 = 10 ** (-18.26)
k3 = 10 ** (-35.91)
logk = 11.96
pH_x = np.linspace(0, 14, 100)
ni_total = 0.1
n_total = 1.2
T_ = 298

# 重构方程组：针对单个pH值求解nh4和nio2
def equations(p, pH):
    nh4, nio2 = p
    h = 10 ** (-pH)
    f = 10 ** (logk - 2 * pH)
    ni2pfree = f/(1 + f / ni_total)
    nh3 = k1 * nh4 / h
    nin4p2 = k2 * (nh4 ** 4) / (h ** 2)
    nin6p2 = k3 * (nh4 ** 6) / (h ** 4)
    return [
        n_total - nh3 - 4*nin4p2 -6*nin6p2 - nh4,
        ni_total - ni2pfree - nin4p2 - nin6p2 - nio2
    ]

# 遍历每个pH值求解参数
nh4_arr = np.zeros_like(pH_x)
nio2_arr = np.zeros_like(pH_x)

for i, pH in enumerate(pH_x):
    guess = (0.1, 0.1)
    res = scipy.optimize.fsolve(equations, guess, args=(pH,))
    nh4_arr[i] = res[0]
    nio2_arr[i] = res[1]

# 计算ni2pfree等变量
f = 10 ** (logk - 2 * pH_x)
ni2pfree = f / (1 + f / ni_total)
nh3 = k1 * nh4_arr / (10 ** (-pH_x))
nin4p2 = k2 * (nh4_arr ** 4) / (10 ** (-pH_x)) ** 2
nin6p2 = k3 * (nh4_arr ** 6) / (10 ** (-pH_x)) ** 4

# 生成曲线
y1 = -0.2405 + 0.0296 * np.log10(ni2pfree)
y2 = 0.1102 - 0.0592 * pH_x

# 插值
interp1 = scipy.interpolate.InterpolatedUnivariateSpline(pH_x, y1)
interp2 = scipy.interpolate.InterpolatedUnivariateSpline(pH_x, y2)

# 找交点：直接返回差值而非绝对值
def difference(pH):
    return interp1(pH) - interp2(pH)

# 多设初始猜测值，避免漏交点
x_candidates = [2, 6, 10]
x_at_crossing = []
for x0 in x_candidates:
    root = scipy.optimize.fsolve(difference, x0=x0)
    if 0 <= root <=14:
        x_at_crossing.append(root[0])

# 去重保证结果唯一
x_at_crossing = np.unique(np.round(x_at_crossing, decimals=4))

# 绘图
plt.plot(pH_x, y1, label='y1')
plt.plot(pH_x, y2, label='y2')
for x in x_at_crossing:
    plt.plot(x, interp1(x), 'cd', ms=7)
plt.legend()
plt.xlabel('pH')
plt.ylabel('Value')
plt.show()

关键说明 #

• 用fsolve替代least_squares处理二元方程组更直接，因为least_squares更适合超定方程组，而这里是恰好两个方程两个未知数。
• 遍历每个pH值求解，确保每个点的参数都是对应pH下的解，符合物理意义。
• 交点求解时多设几个初始猜测值，避免错过多个交点；最后去重保证结果唯一。

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