使用scipy的fsolve求解含数组变量方程时报错如何解决？

• 核心报错根源是x为长度50的数组，但fsolve要求目标函数返回的残差长度必须和输入未知数的个数一致（本次求解2个未知数，残差需要为2个元素）。你原代码中第一个返回项(xa+xb)-x是长度50的数组，第二个返回项是标量，形状完全不匹配，因此触发类型错误。
• 你的需求是对x数组中的每一个值分别求解对应的xa、xb，不能直接把整个数组丢进方程中一次性计算。

from scipy.optimize import fsolve
import numpy as np
x_arr = np.array([0.01394346, 0.01727623, 0.02139363, 0.02647018, 0.03271463,
       0.04037362, 0.04974189, 0.06113669, 0.07489298, 0.0914198 ,
       0.11119803, 0.13467519, 0.1622032 , 0.19409497, 0.23056969,
       0.27160588, 0.3169306 , 0.36594778, 0.41776317, 0.471372  ,
       0.5256699 , 0.5794641 , 0.63155216, 0.68075466, 0.72615373,
       0.767264  , 0.80389225, 0.8359975 , 0.8637226 , 0.8873593 ,
       0.9072883 , 0.92398256, 0.93792963, 0.9494987 , 0.95898896,
       0.96673596, 0.9730703 , 0.9782413 , 0.982436  , 0.9858274 ,
       0.9885736 , 0.99079704, 0.9925912 , 0.99403596, 0.99519914,
       0.99613607, 0.9968898 , 0.99749595, 0.9979841 , 0.99837816])

# 提前计算固定常数，避免重复运算
R = 8.314
F = 96485
T = 298.15
E0_A = 0.516
E0_B = 0.393
CA = 0.37
CB = 0.63
K = (CA/CB)*np.exp(-(F/(R*T))*(E0_B - E0_A))

# 方程定义，传入单个x值作为参数
def eqs(p, x_single):
    xa, xb = p
    return [
        xa + xb - x_single,
        xa/xb - K
    ]

# 遍历每个x值求解
result = []
for x in x_arr:
    xa_sol, xb_sol = fsolve(eqs, x0=(1, 1), args=(x,))
    result.append((xa_sol, xb_sol))
result = np.array(result)
print(result)

该方程组可以直接推导解析解，完全不需要调用数值求解函数，既避免收敛问题，运行效率也更高：
由第二个方程可得xa = K * xb，代入第一个方程得xb*(K+1) = x，直接向量化计算即可：

K = (0.37/0.63)*np.exp(-(96485/(8.314*298.15))*(0.393-0.516))
xa = K * x_arr / (K + 1)
xb = x_arr / (K + 1)

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