我试图通过在 python 2.7 中使用 scipy.optimize.fsolve
找到方程组的解。目标是计算化学系统的平衡浓度。由于问题的性质,一些常量非常小。现在对于某些组合,我确实得到了适当的解决方案。对于某些参数,我找不到解决方案。要么解是负的,这从物理角度来看是不合理的,要么 fsolve 产生:
ier = 3, 'xtol=0.000000 太小,无法进一步改进\n近似解。')
ier = 4, '迭代没有取得很好的进展,正如最近五次雅可比评估的\n 改进所衡量的那样。')
ier = 5, '迭代没有取得很好的进展,根据最近十次迭代\n 的改进来衡量。')
根据我的研究,在我看来,未能找到方程组的正确解与数据类型 float.64
不够精确有关。正如一位 friend 指出的那样,系统条件不佳,参数相差几个数量级。
所以我尝试将 fsolve 与 gmpy2
模块提供的 mpfr
类型一起使用,但这导致了以下错误:
TypeError:无法根据规则“安全”将数组数据从 dtype('O') 转换为 dtype('float64')
现在这里有一个带有参数的小例子,如果随机化的起始参数恰好合适,它会导致一个解决方案。然而,如果常数 C_HCL 被选择为类似 1e-4 或更大的东西,那么我永远找不到合适的解决方案。
from numpy import *
from scipy.optimize import *
K_1 = 1e-8
K_2 = 1e-8
K_W = 1e-30
C_HCL = 1e-11
C_NAOH = K_W/C_HCL
C_HL = 1e-6
if C_HCL-C_NAOH > 0:
Saeure_Base = C_HCL-C_NAOH+sqrt(K_W)
OH_init = K_W/(Saeure_Base)
elif C_HCL-C_NAOH < 0:
OH_init = C_NAOH-C_HCL+sqrt(K_W)
Saeure_Base = K_W/OH_init
# some randomized start parameters
G1 = random.uniform(0, 2)*Saeure_Base
G2 = random.uniform(0, 2)*OH_init
G3 = random.uniform(1, 2)*C_HL*(sqrt(K_W))/(Saeure_Base+OH_init)
G4 = random.uniform(0.1, 1)*(C_HL - G3)/2
G5 = C_HL - G3 - G4
zGuess = array([G1,G2,G3,G4,G5])
#equation system / 5 variables --> H3O, OH, HL, H2L, L
def myFunction(z):
H3O = z[0]
OH = z[1]
HL = z[2]
H2L = z[3]
L = z[4]
F = empty((5))
F[0] = H3O*L/HL - K_1
F[1] = OH*H2L/HL - K_2
F[2] = K_W - OH*H3O
F[3] = C_HL - HL - H2L - L
F[4] = OH+L+C_HCL-H2L-H3O-C_NAOH
return F
z = fsolve(myFunction,zGuess, maxfev=10000, xtol=1e-15, full_output=1,factor=0.1)
print z
所以问题是。这个问题是基于 float.64 的精度和 如果是,(如何)可以用 python 解决? fsolve 是要走的路吗?我是否需要更改 fsolve 函数以使其接受不同的数据类型?
最佳答案
问题的根源要么是理论上的,要么是数值上的。
scipy.optimize.fsolve
函数基于 MINPACK Fortran 求解器 (http://www.netlib.org/minpack/)。该求解器使用 Newton-Raphson 优化算法来提供解决方案。
当您使用此算法时,存在关于函数平滑度的基本假设。例如,解点 x
处的雅可比矩阵应该是可逆的。您更关心的是吸引力盆地。
为了收敛,算法的起点需要靠近实际解,即在吸引盆地中。凸函数始终满足此条件,但是很容易找到该算法对某些函数表现不佳的函数。您的函数就是其中之一,因为您只有一小部分输入参数。
要解决此问题,您只需更改起点即可。这个起点对于具有多个解决方案的函数也变得非常重要:this picture来自维基百科的文章向您展示了根据起点找到的解决方案(五种解决方案的五种颜色);所以你应该小心你的解决方案并实际检查你的解决方案的“物理”方面。
对于数值方面,Newton-Raphson 算法需要有雅可比矩阵(导数矩阵)的值。如果未将其提供给 MINPACK 求解器,则使用有限差分公式估计雅可比矩阵。需要提供有限差分公式的扰动步长epsfcn=None
,None
仅在fprime
的情况下作为默认值提供(在这种情况下不需要雅可比估计)。所以首先你应该合并它。您也可以通过手动推导函数来直接指定雅可比矩阵。
但是,步长的最小值将是机器精度,也称为 machine epsilon .对于您的问题,您的输入值非常小,这可能是个问题。我建议将它们中的每个乘以相同的值(如 10^6),这相当于更改单位,但会避免舍入误差和机器精度问题。
当您查看您提供的参数 xtol=1e-15
时,这个问题也很重要。在您的错误消息中,它给出了 xtol=0.000000
,因为它低于机器精度并且不能被考虑在内。此外,如果您查看行 F[2] = K_W - OH*H3O
,考虑到机器精度,K_W
是否为 1e- 并不重要15
或 1e-30
。与机器精度相比,0
是针对这两种情况的解决方案。为避免此问题,只需将所有内容乘以一个更大的值即可。
总结一下:
- 对于 Newton-Raphson 算法,初始化点很重要!
- 对于此算法,您应该指定计算雅可比矩阵的方式!
- 在数值计算中,切勿使用小值。您可以轻松地将尺寸更改为不同的东西:它是基本单位转换,例如以克而不是千克为单位。
关于python - 带 fsolve 的方程系统,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/46423069/