代码之家 › 专栏 › 技术社区 › redfol

具有fsolve的方程组

precision scipy python

redfol · 技术社区 · 7 年前

我试图通过使用 scipy.optimize.fsolve
ier=3,“xtol=0.000000太小,近似解不可能进一步改进。”)
ier=4,“从最近五次雅可比矩阵评估的改进来看,迭代没有取得良好的进展。”)

根据我的研究,在我看来,未能找到方程系统的正确解与数据类型有关 float.64 所以我试着用fsolve和 mpfr 由提供的类型 gmpy2 但这导致了以下错误:

现在,这里有一个关于参数的小示例,如果随机起始参数恰好适合,则会得出一个解决方案。然而,如果常数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

所以问题是。这个问题是基于浮点的精度吗。64和如果是,(如何)用python解决?fsolve是一条路吗?我是否需要更改fsolve函数,使其接受不同的数据类型?

1 回复 | 直到 7 年前

Nathan 7 年前

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

scipy.optimize.fsolve 函数基于MINPACK Fortran求解器( http://www.netlib.org/minpack/ ). 该求解器使用牛顿-拉斐逊优化算法来提供解。

x 应该是可逆的。你更关心的是吸引力的盆地。为了收敛,算法的起点需要接近实际解,即在吸引盆中。凸函数总是满足这个条件,但是很容易找到一些该算法性能较差的函数。你的函数就是其中之一,因为你有一小部分的输入参数。

为了解决这个问题,你应该改变起点。这个起点对于具有多个解决方案的函数也变得非常重要: this picture 维基百科的文章向您展示了根据起点找到的解决方案(五种颜色代表五种解决方案);因此,你应该小心你的解决方案,并实际检查解决方案的“物理”方面。

对于数值方面,牛顿-拉斐逊算法需要具有雅可比矩阵(导数矩阵)的值。如果未提供给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 与机器精度相比,这两种情况都是一种解决方案。为了避免这个问题,只需将所有值乘以一个更大的值。

对于牛顿-拉弗森算法,初始化点很重要!
对于该算法,您应该指定如何计算雅可比矩阵!