带解的方程组 [英] equation system with fsolve

问题描述

我尝试通过使用python 2.7中的scipy.optimize.fsolve来找到方程组的解决方案.目的是计算化学系统的平衡浓度.由于问题的性质，某些常数非常小.现在，对于某些组合，我确实获得了适当的解决方案.对于某些参数，我找不到解决方案.解决方案是否定的，从物理角度来看这是不合理的，或者解决方法是产生以下结果:
ier = 3，'xtol = 0.000000太小，无法进一步改善近似\ n解.')
ier = 4，根据最近五次Jacobian评估对\ n的改进，迭代没有取得很好的进展.")
ier = 5，根据最近十次迭代的\ n改进，该迭代进展不佳."

I try to find a solution for a system of equations by using scipy.optimize.fsolve in python 2.7. The goal is to calculate equilibrium concentrations for a chemical system. Due to the nature of the problem, some of the constants are very small. Now for some combinations i do get a proper solution. For some parameters i don't find a solution. Either the solutions are negative, which is not reasonable from a physical point of view or fsolve produces:
ier = 3, 'xtol=0.000000 is too small, no further improvement in the approximate\n solution is possible.')
ier = 4, 'The iteration is not making good progress, as measured by the \n improvement from the last five Jacobian evaluations.')
ier = 5, 'The iteration is not making good progress, as measured by the \n improvement from the last ten iterations.')

在我看来，根据我的研究，未能找到方程式系统的正确解的原因与数据类型float.64不够精确.正如一位朋友指出的那样，该系统在参数大小不同的情况下条件不佳. 因此，我尝试将fsolve与gmpy2模块提供的mpfr类型一起使用，但是导致以下错误:
TypeError:根据规则安全"，无法将数组数据从dtype('O')转换为dtype('float64')

It seems to me, based on my research, that the failure to find proper solutions of the equation system is connected to the datatype float.64 not being precise enough. As a friend pointed out, the system is not well conditioned with parameters differing in several magnitudes. So i tried to use fsolve with the mpfr type provided by the gmpy2 module but that resulted in the following error:
TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'

现在这是一个带有参数的小示例，如果随机化的起始参数恰好很好，则可以得出解决方案.但是，如果将常量C_HCL选择为1e-4或更大，则我永远找不到合适的解决方案.

Now here is a small example with parameter which lead to a solution if the randomized starting parameters fit happen to be good. However if the constant C_HCL is chosen to be something like 1e-4 or bigger then i never find a proper solution.

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函数，使其接受其他数据类型?

So the questions are. Is this problem based on the precision of float.64 and if yes , (how) can it be solved with python? Is fsolve the way to go? Would i need to change the fsolve function so it accepts a different data type?

推荐答案

问题的根源是理论上的还是数字上的.

The root of your problem is either theoretical or numerical.

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

The scipy.optimize.fsolvefunction is based on the MINPACK Fortran solver (http://www.netlib.org/minpack/). This solver use a Newton-Raphson optimisation algorithm to provide the solution.

使用此算法时，对函数的平滑度有一些基本假设.例如，在求解点x上的雅可比矩阵应该是可逆的.您最关心的是吸引盆地. 为了收敛，算法的起点必须在实际解附近，即在吸引盆地中.凸函数始终满足此条件，但是很容易找到该算法表现不佳的某些函数.您的函数就是其中之一，因为您只有一部分输入参数.

There are underlying assumptions about the smoothness of the function when you use this algorithm. For example, the jacobian matrix at the solution point x is supposed to be invertible. The one you are more concerned about is the basins of attraction. In order to converge, the starting point of the algorithm needs to be near the actual solution, i.e. in the basins of attraction. This condition is always met for convex functions, however it is easy to find some functions for which this algorithm behaves badly. Your function is one of this as you have a fraction of your inputs parameters.

要解决此问题，您应该更改起点.对于具有多种解决方案的功能，此起点也变得非常重要:维基百科文章中的这张图片向您展示了根据起点找到的解决方案(五种解决方案中的五种颜色)；因此您应该谨慎使用解决方案，并实际检查解决方案的物理"方面.

To address this issue you should just change the starting point. This starting point becomes also very important for functions with multiple solutions: this picture from the wikipedia article shows you the solution found depending of the starting point (five colours for five solutions); so you should be careful with your solution and actually check the "physical" aspects of your solution.

对于数值方面，Newton-Raphson算法需要具有雅可比矩阵(导数矩阵)的值.如果未将其提供给MINPACK求解器，则将使用有限差分公式估算雅可比.需要提供epsfcn=None的有限差分公式的摄动步骤，仅在提供fprime的情况下(此处无需进行雅可比估计)，None在此处为默认值.所以首先您应该将其合并.您也可以通过手动派生函数直接指定jacobian.

For the numerical aspects, the Newton-Raphson algorithm needs to have the value of the jacobian matrix (the derivatives matrix). If it is not provided to the MINPACK solver, the jacobian is estimated with a finite-difference formula. The perturbation step for the finite difference formula need to be provided epsfcn=None, the None being here as default value only in the case where fprimeis provided (there is no need for the jacobian estimation in this case). So first you should incorporate that. You could also specify directly the jacobian by derivating your function by hand.

但是，步长的最小值将是机器精度，也称为机器epsilon .对于您的问题，您的输入值非常小，这可能是一个问题.我建议将每个人都乘以相同的值(例如10 ^ 6)，这等效于单位的更改，但将避免舍入错误和机器精度问题.

However, the minimum value for the step size will be the machine precision, also called machine epsilon. For your problem, you have very small inputs values which can be a problem. I would suggest multiply everyone of them by the same value (like 10^6), it is equivalent to a change of the units but will avoid rounding up errors and problems with machine precision.

当您查看提供的参数xtol=1e-15时，此问题也很重要.在您的错误消息中，它给出xtol=0.000000，因为它低于机器精度，因此无法考虑.同样，如果您查看行F[2] = K_W - OH*H3O，在给定机器精度的情况下，K_W是1e-15还是1e-30都没有关系.相对于机器精度，0是这两种情况的解决方案.为避免此问题，只需将所有值乘以更大的值即可.

This problem is also important when you look at the parameter xtol=1e-15 you provided. In your error message, it gives xtol=0.000000, as it is below machine precision and cannot be taken into account. Also, if you look at your line F[2] = K_W - OH*H3O, given the machine precision, it does not matter if K_W is 1e-15or 1e-30. 0 is a solution for both of this case compare to the machine precision. To avoid this problem, just multiply everything by a bigger value.

所以总结一下:

1. 对于Newton-Raphson算法，初始化点很重要！
2. 对于此算法，您应指定如何计算jacobian！
3. 在数值计算中，切勿使用较小的值.您可以轻松地将尺寸更改为其他尺寸:这是基本单位转换，例如以克而不是千克为单位.

这篇关于带解的方程组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！