含fs方程组

Question

我试图通过在Python2.7中使用scipy.optimize.fsolve来找到方程组的解。目标是计算化学系统的平衡浓度。由于问题的性质，有些常数非常小。对于某些组合，我确实找到了一个合适的解决方案。对于某些参数，我找不到解决方案。要么解决方案是否定的，从物理角度来看这是不合理的，要么fsolve产生了：
ier=3，'xtol=0.000000太小，无法进一步改进近似解决方案。'）
ier=4，“迭代没有取得很好的进展，这是由最近五次雅可比评估的改进来衡量的。”）
ier=5，'迭代没有取得很好的进展，这是通过\n最近十次迭代的改进来衡量的。'）

在我看来，根据我的研究，未能找到方程组的正确解与数据类型float.64不够精确有关。正如一位朋友所指出的那样，系统的状态不是很好，参数有几个量级的不同。 因此，我尝试将fsolve与gmpy2模块提供的mpfr类型一起使用，但结果导致以下错误：
TypeError:无法根据规则“safe”

这里有一个小例子，如果随机启动参数拟合得很好的话，这个例子会导致一个解决方案。但是，如果常数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解决它？F解决方案是正确的吗？我是否需要更改fsolve函数以便它接受不同的数据类型？在

Answer 1

你问题的根源不是理论上的就是数值上的。在

scipy.optimize.fsolve函数基于MINPACK Fortran解算器（http://www.netlib.org/minpack/）。这个解算器使用牛顿-拉夫森优化算法来提供解。在

使用此算法时，有一些关于函数平滑性的基本假设。例如，解点x处的雅可比矩阵应该是可逆的。你更关心的是吸引力的盆地。 为了收敛，算法的起点必须接近实际解，即在吸引域内。凸函数总是满足这一条件，但很容易找到一些算法性能较差的函数。你的函数就是其中之一，因为你有一小部分输入参数。在

为了解决这个问题，你应该改变出发点。对于具有多个解决方案的函数，这个起点也变得非常重要：来自wikipedia文章的this picture向您展示了根据起始点找到的解决方案（五种解决方案对应五种颜色）；因此您应该谨慎对待您的解决方案，并实际检查解决方案的“物理”方面。在

在数值方面，Newton-Raphson算法需要有雅可比矩阵（导数矩阵）的值。如果未提供给MINPACK解算器，则使用有限差分公式估计雅可比。需要提供有限差分公式的扰动步骤epsfcn=None，只有在提供{}的情况下，None作为默认值（在这种情况下不需要雅可比估计）。所以首先你应该把它合并起来。您也可以通过手动派生函数来直接指定jacobian。在

但是，步长的最小值是机器精度，也称为machine epsilon。对于您的问题，您有非常小的输入值，这可能是一个问题。我建议用相同的值（比如10^6）乘以它们，这相当于单位的改变，但可以避免误差和机器精度的问题。在

当您查看您提供的参数xtol=1e-15时，这个问题也很重要。在您的错误消息中，它给出xtol=0.000000，因为它低于机器精度，不能被考虑在内。另外，如果你看你的行F[2] = K_W - OH*H3O，给定机器精度，K_W是1e-15还是{}，这并不重要。0是这两种情况下与机器精度比较的解决方案。为了避免这个问题，只需将所有值乘以一个更大的值。在

所以总结一下：

1. 对于Newton-Raphson算法，初始化点很重要！在
2. 对于这个算法，您应该指定如何计算雅可比！在
3. 在数值计算中，不要使用小数值。你可以很容易地把尺寸改成不同的尺寸：这是基本单位转换，比如用克代替千克。在