python—如何调优fsolve以避免“在日志中遇到无效值”

jgwigjjp 于 2021-08-25 发布在 Java

关注(0)|答案(1)|浏览(599)

我只想解这里提到的氨溶液的ph值。
这里详细介绍了完整的方程式，下面是我的数值解：

import numpy as np
import scipy.optimize as optimize

C = np.array([0, 0, 0.1, 0])
N = np.array([
    [1, 1, 0, 0],
    [0, 1, -1, 1]
])
K = np.array([1E-14, 1.77E-5])

def eq1(X):
    return N.transpose() @ np.array(X) + C

def eq2(X):
    return N @ np.log(eq1(X)) - np.log(K)

X = optimize.fsolve(eq2, [1E-14, 1E-14])
print("X=",X)
print("Y=",eq1(X))

但我得到了这些警告：

RuntimeWarning: invalid value encountered in log
  return N @ np.log(eq1(X)) - np.log(K)
/usr/local/lib64/python3.6/site-packages/scipy/optimize/minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)
X= [2.39338797e-06 2.27841240e-05]
Y= [2.39338797e-06 2.51775120e-05 9.99772159e-02 2.27841240e-05]

这个问题对于fsolve来说很棘手，因为在正确的解决方案中 [7.57E-12 1.32E-3] x1太接近零。
我尝试了其他初始值 [1E-7, 1E-7] , [1E-1, 1E-1] ，等等，没有人能避免这一点，直到我尝试了接近正确答案的配对，例如。 [1E-12, 1E-3] .
fsolve是否有任何参数可以避免此类问题？或者有什么方法可以给出一个合适的初始值？

python numpy scipy numerical-methods chemistry

来源：https://stackoverflow.com/questions/68302523/how-can-i-tune-fsolve-to-avoid-invalid-value-encountered-in-log

1条答案

按热度按时间

o75abkj41#

试试这个安排

import numpy as np
import scipy.optimize as optimize

C = np.array([0, 0, 0.1, 0])
N = np.array([
    [1, 1, 0, 0],
    [0, 1, -1, 1]
])

K = np.array([1e-14, 1.77e-5])

x1 = 1e-2
x2 = 2e-2
y1 = 3e-2
y2 = 4e-2
y3 = 5e-2
y4 = 6e-2
xmin = 1e-10
xmax = 1e0
XY = np.array([x1,x2,y1,y2,y3,y4])

def eq1(XY):
    X = np.array([XY[0],XY[1]])
    Y = np.array([XY[2],XY[3],XY[4],XY[5]]) 
    return N.transpose() @ np.array(X) + C - Y

def eq2(XY):
    Y = np.array([XY[2],XY[3],XY[4],XY[5]])
    return N @ np.log(abs(Y)) - np.log(K)

def obj(XY):
    return np.linalg.norm(eq1(XY),2) + np.linalg.norm(eq2(XY),2)

bounds = ((xmin,xmax),(xmin,xmax),(xmin,xmax),(xmin,xmax),(xmin,xmax),(xmin,xmax))

result = optimize.minimize(obj, XY, bounds = bounds, method = 'trust-constr')

print(result.x)

print(obj(result.x))

与剩余验证

print(eq1(result.x))
print(eq2(result.x))

赞(0）回复(0）举报 2021-08-25

我来回答

python—如何调优fsolve以避免“在日志中遇到无效值”

1条答案

相关问题

热门标签

最新问答