为何fsolve对两个等价的数值求解函数返回不同结果？

我需要对单变量方程进行数值求解，编写了两个逻辑上等价的函数func1和func2——指定变量值时两者返回的方程残差完全相同，但用fsolve自动求解时却得到不同结果。

相关数据与函数定义 #

基础数据 #

import numpy as np
from scipy.optimize import fsolve

dHi = np.array([-125790,49080,4.2])     #  [n butane :  1,3 butadiene  :   H2]
dGi = np.array([-16570,124520,17.6])
V = np.array([-1,1,2])
No = np.array([1,0,0])

dH = sum(V*dHi)
dG = sum(V*dGi)

函数1定义 #

def func1(e):
    R=8.314
    T1=298
    T2=925
    Nt=1+2*e
    Ni=np.array([1-e,e,2*e])
    
    lnk1=(-dG/(R*T1))
    lnk2=-(dH/R)*(1/T2 - 1/T1)+lnk1
    k2=np.exp(lnk2)
    
    A1=np.prod((Ni/Nt)**V)-k2
    
    return A1

函数2及辅助函数定义 #

def N(e):
    return No+e*V

def func2(e):
    R=8.314
    T1=298
    T2=925
    Nt=1+2*e
     
    lnk1=(-dG/(R*T1))
    lnk2=-(dH/R)*(1/T2 - 1/T1)+lnk1
    k2=np.exp(lnk2)
    
    A1=np.prod(((N(e))/Nt)**V)-k2
    
    return A1

验证与求解结果 #

验证Ni与N(e)等价 #

当传入标量e时，两者输出完全一致：

e=0.1
Ni=np.array([1-e,e,2*e])
print(Ni,N(e))

输出：

[0.9 0.1 0.2] [0.9 0.1 0.2]

fsolve求解结果 #

print(fsolve(func1,0.03), fsolve(func2,0.03))

输出：

[0.10045184] [0.03108138]

验证解的正确性 #

e=0.03108138是方程的正确解，两个函数代入后残差一致：

print(func1(0.03108138),func2(0.03108138))

输出：

1.2794325793047054e-11 1.2794325793047054e-11

问题根源 #

fsolve在迭代求解过程中，会传递数组类型的参数（用于数值微分计算），此时两个函数的内部运算出现本质差异：

• func1中Ni = np.array([1-e,e,2*e])当e是形状为(1,)的数组时，会生成二维数组（例如array([[0.97], [0.03], [0.06]])）；
• func2中N(e) = No + e*V依赖numpy的广播机制，当e为数组时，仍会生成一维数组（例如array([0.97, 0.03, 0.06])）。

在计算np.prod((Ni/Nt)**V)时，二维数组的prod()默认对所有元素求积，而一维数组的prod()是按元素依次求积，两者结果完全不同。这导致fsolve获取错误的函数值和导数信息，最终收敛到错误的解。

修复方案 #

修改func1中Ni的生成逻辑，使其与N(e)保持一致，确保无论e是标量还是数组，都生成一维数组：

def func1(e):
    R=8.314
    T1=298
    T2=925
    Nt=1+2*e
    # 改用广播运算生成一维数组，和N(e)逻辑一致
    Ni = No + e*V
    
    lnk1=(-dG/(R*T1))
    lnk2=-(dH/R)*(1/T2 - 1/T1)+lnk1
    k2=np.exp(lnk2)
    
    A1=np.prod((Ni/Nt)**V)-k2
    
    return A1

修改后，两个函数通过fsolve求解会得到一致的正确结果。

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