python - 带 fsolve 的方程系统

Question

我试图通过在 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 函数以使其接受不同的数据类型？

Answer 1

问题的根源要么是理论上的，要么是数值上的。

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 是针对这两种情况的解决方案。为避免此问题，只需将所有内容乘以一个更大的值即可。

总结一下:

1. 对于 Newton-Raphson 算法，初始化点很重要!
2. 对于此算法，您应该指定计算雅可比矩阵的方式!
3. 在数值计算中，切勿使用小值。您可以轻松地将尺寸更改为不同的东西:它是基本单位转换，例如以克而不是千克为单位。