Python 的 fsolve 不工作

Question

我目前正在尝试从我的代码（粘贴在下面）中找到 2 个方程的截距。我正在使用 fsolve 并且在一部分中成功地使用了它，但我无法让它在第二部分工作。

令人困惑的是，它没有显示错误，如果您将此代码粘贴到笔记本中，然后运行您会看到 2 个图表，在第一个图表上有一条斜线应该停在均衡线。

不起作用的部分是 def q_eqm(x_q)。感谢您的帮助

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

AC_LK = np.array([4.02232,1206.53,220.291])
AC_HK = np.array([4.0854,1348.77,219.976])

P_Tot = 1 # Bara
N_Size = 11 # 1001 = 0.1% accuracy for xA

xf = 0.7
q = 0.7

xA = np.linspace(0,1,N_Size)
yA = np.linspace(0.00,0.00,N_Size)
T = np.linspace(0.00,0.00,N_Size)

x = np.array([xA[0:N_Size],yA[0:N_Size],T[0:N_Size]]) # x[xA,yA,T]

F = np.empty((1))

def xA_T(N):
    xA_Ant = x[0,N]
    def P_Ant(T):

        PA = pow(10,AC_LK[0]-(AC_LK[1]/(T+AC_LK[2])))*xA_Ant
        PB = pow(10,AC_HK[0]-(AC_HK[1]/(T+AC_HK[2])))*(1-xA_Ant)

        F[0] = P_Tot - (PA + PB)

        return F
        return x
    TGuess = [100]
    T = opt.fsolve(P_Ant,TGuess)

    x[2,N] = T

    return x

for N in range(0,len(xA)):
    xA_T(N)
    x[1,N] = pow(10,AC_LK[0]-(AC_LK[1]/(x[2,N]+AC_LK[2])))*x[0,N]/P_Tot

q_int = ((-q*0)/(1-q)) + (xf/(1-q))

Eqm_Poly = np.polyfit(x[0,0:N_Size], x[1,0:N_Size], 6)
q_Poly = np.polyfit([xf,0], [xf,q_int], 1)

F = np.empty((1))

def q_Eqm(x_q):

    y_q = q_Poly[0]*x_q + q_Poly[1]
    eqm_y = (Eqm_Poly[0]*pow(x_q,6)+Eqm_Poly[1]*pow(x_q,5)+Eqm_Poly[2]*pow(x_q,4)+Eqm_Poly[3]*pow(x_q,3)+Eqm_Poly[4]*pow(x_q,2)+Eqm_Poly[5]*pow(x_q,1)+Eqm_Poly[6]*pow(x_q,0))

    F[0] = y_q - eqm_y 
    return F

x_qGuess = [0]
x_q = opt.fsolve(q_Eqm,x_qGuess)


print(x,Eqm_Poly,x_q,q_int)

plt.plot(x[0,0:N_Size],x[1,0:N_Size],'k-',linewidth=1)
plt.plot([xf,xf],[0,xf],'b-',linewidth=1)
plt.plot([xf,x_q],[xf,(q_Poly[0]*x_q + q_Poly[1])],'r-',linewidth=1)
plt.legend(['Eqm','Feed'])
plt.xlabel('xA')
plt.ylabel('yA')
plt.xlim([0.00, 1])
plt.ylim([0.00, 1])
plt.savefig('x.png')
plt.savefig('x.eps')
plt.show()

plt.plot(x[0,0:N_Size],x[2,0:N_Size],'r--',linewidth=3)
plt.plot(x[1,0:N_Size],x[2,0:N_Size],'b--',linewidth=3)
plt.legend(['xA','yA'])
plt.xlabel('Mol Frac')
plt.ylabel('Temp degC')
plt.xlim([0, 1])
plt.savefig('Txy.png')
plt.savefig('Txy.eps')
plt.show()

Answer 1

原来答案比较简单：

#F = np.empty((1))  # remove this

def q_Eqm(x_q):
    y_q = q_Poly[0]*x_q + q_Poly[1]
    eqm_y = (Eqm_Poly[0]*pow(x_q,6)+Eqm_Poly[1]*pow(x_q,5)+Eqm_Poly[2]*pow(x_q,4)+Eqm_Poly[3]*pow(x_q,3)+Eqm_Poly[4]*pow(x_q,2)+Eqm_Poly[5]*pow(x_q,1)+Eqm_Poly[6]*pow(x_q,0))
    return y_q - eqm_y

原代码定义了一个全局的F，在函数中修改，然后返回。所以在每次迭代中，函数 returns 不同的值，但它们是 相同的 对象。这似乎混淆了 fsolve（我猜它在内部存储对结果而不是值的引用）。删除此 F 并简单地返回减法的结果即可解决问题。

Python 的 fsolve 不工作

Python's fsolve not working

python

numpy

solver

scipy