Python的fsolve无法工作

Question

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

令人困惑的是，它没有显示错误，如果你把这个代码粘贴到你的笔记本上，然后运行它，你会看到2个grph，在第一个图上，有一条线在某个角度，应该停在eqm线上。在

不起作用的部分是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，它在函数中被修改，然后返回。因此，在每次迭代中，函数返回不同的值，但它们是相同的对象。这似乎混淆了fsolve（我想它内部存储对结果的引用，而不是值）。删除这个F并返回减法的结果就解决了这个问题。在