时间戳

首页 > 解决方案 > Sympy nsolve 用于联立方程化学平衡

python - Sympy nsolve 用于联立方程化学平衡

问题描述

我正在尝试使用 Sympy 的 nsolve 作为一种数值方法来解决同时化学平衡中的反应进程变量。不幸的是，对于目前的情况，似乎没有多少初始猜测可以得出令人满意的解决方案。我相信这个系统确实有一个合适的解决方案，其中所有四个进度变量（未知数）都有介于 0 和 1 之间的真实解决方案。

这是简化的代码：

from sympy import *
var('x1,x2,x3,x4')
P = 1.          # Pressure [bar]
n1 = 2.*x1      # x1 is progress variable for reaction 1
n2 = x2         # x2 is progress variable for reaction 2
n3 = 2.*x3      # x3 is progress variable for reaction 3
n4 = x4         # x4 is progress variable for reaction 4
n5 = 1.-x1-x2-3.*x3-2.*x4
nTot = n1+n2+n3+n4+n5

rxn1 = Eq(P*n1**2/(nTot*n5),3.06e-3)
rxn2 = Eq(n2/n5,8.02e5)
rxn3 = Eq(nTot*n3**2/(P*n5**3),1.07e21)
rxn4 = Eq(nTot*n4/(P*n5**2),8.25e14)

reactions = (rxn1,rxn2,rxn3,rxn4)
unknowns = (x1,x2,x3,x4)
answer = nsolve(reactions,unknowns,(1e-4,.1,.3,.2))
print (answer)

以下是 Sympy 的输出（我在具有 64 位 Python 2.7 的 Spyder3 环境中使用 Sympy v1.1.1）：

Traceback (most recent call last):

  File "<ipython-input-9-43632f4fa60f>", line 18, in <module>
    answer = nsolve(reactions,unknowns,(1e-4,.1,.3,.2))

  File "C:\Users\<user>\AppData\Local\Continuum\anaconda2\lib\site-packages\sympy\utilities\decorator.py", line 91, in func_wrapper
    return func(*args, **kwargs)

  File "C:\Users\<user>\AppData\Local\Continuum\anaconda2\lib\site-packages\sympy\solvers\solvers.py", line 2847, in nsolve
    x = findroot(f, x0, J=J, **kwargs)

  File "C:\Users\<user>\AppData\Local\Continuum\anaconda2\lib\site-packages\mpmath\calculus\optimization.py", line 975, in findroot
    % (norm(f(*xl))**2, tol))

ValueError: Could not find root within given tolerance. (1.1449e+42 > 2.1684e-19)
Try another starting point or tweak arguments.

我尝试使用各种求解（不出所料，没有可用的解析解）和 nsolve 方法，到目前为止，对于这个非线性方程组，这些方法都没有成功。我也尝试过以对数形式重铸方程，得到类似的结果。

您是否建议我尝试使用其他求解器设置？

感谢您的任何帮助！戴夫

标签： pythonsympy

解决方案

找到解决方案的部分困难是由于常数大小的广泛变化。例如，rxn3右侧的值为 1070000000000000000000，而rxn1右侧的值为 0.00306。
数值算法通常在处理具有极端浮点值的方程时遇到问题。

尽管如此，我们可以通过符号求解n1, n2, n3,的 4 个方程来获得一些接近n4。（根据 - 变量找到解决方案n与求解 - 变量基本相同x，但允许我们用更少的变量/方程解决主要问题）。

soln = sym.solve([rxn1, rxn2, rxn3, rxn4], [n1, n2, n3, n4], dict=True)[0]
# {n1: -0.0553172667437573*sqrt(n5*nTot),
#  n2: 802000.0*n5,
#  n3: -32710854467.5922*sqrt(n5**3/nTot),
#  n4: 825000000000000.0*n5**2/nTot}

这给了我们n1,n2和n3和。我们现在可以尝试求解剩下的两个方程，它们定义了和-和，（见下文） -和。n4n5nTotn5nTotrxn5rxn6n5nTot

不幸的是，sym.solve未能找到一个象征性的解决方案：

In [114]: sym.solve([rxn5, rxn6], [n5, nTot])
Out[114]: []

请注意，当sym.solve返回一个空列表时，它并不一定意味着没有解决方案，只是sym.solve没有找到一个。

所以让我们试着找到一个数值解。根据文档sym.nsolve，

mpmath.findroot 被使用，您可以在那里找到更广泛的文档，特别是关于关键字参数和可用求解器的文档。

以及描述参数的文档。mpmath.findrootmaxstep通过增加maxstep参数，事实证明sym.nsolve可以找到n5和的解决方案nTot：

In [148]: sym.nsolve([rxn5, rxn6], [n5, nTot], [0.5, 0.5])
ValueError: Could not find root within given tolerance. (19920803.9870484600143 > 2.16840434497100886801e-19)
Try another starting point or tweak arguments.

In [149]: sym.nsolve([rxn5, rxn6], [n5, nTot], [0.5, 0.5], maxsteps=50)
Out[149]: 
Matrix([
[1.83437114233209e-8],
[  0.477964294180007]])

现在我们可以检查解决方案是否有意义：

import sympy as sym

P = 1.0          # Pressure [bar]
n1, n2, n3, n4, n5, nTot = sym.symbols('n1,n2,n3,n4,n5,nTot', positive=True, real=True)
rxn1 = sym.Eq(P*n1**2/(nTot*n5),3.06e-3)
rxn2 = sym.Eq(n2/n5,8.02e5)
rxn3 = sym.Eq(nTot*n3**2/(P*n5**3),1.07e21)
rxn4 = sym.Eq(nTot*n4/(P*n5**2),8.25e14)
rxn5 = sym.Eq(n5, 1.0-(n1/2.0)-n2-3.0*(n3/2.0)-2.0*n4)
rxn6 = sym.Eq(nTot, n1+n2+n3+n4+n5)

soln = sym.solve([rxn1, rxn2, rxn3, rxn4], [n1, n2, n3, n4], dict=True)[0]
rxn5 = rxn5.subs(soln)
rxn6 = rxn6.subs(soln)

expr5 = rxn5.rhs - rxn5.lhs
expr6 = rxn6.rhs - rxn6.lhs
print(expr5)
print(expr6)
# -1.65e+15*n5**2/nTot - 802001.0*n5 + 0.0276586333718787*sqrt(n5*nTot) + 49066281701.3884*sqrt(n5**3/nTot) + 1.0
# 825000000000000.0*n5**2/nTot + 802001.0*n5 - nTot - 0.0553172667437573*sqrt(n5*nTot) - 32710854467.5922*sqrt(n5**3/nTot)

root = sym.nsolve([rxn5, rxn6], [n5, nTot], [0.5, 0.5], maxsteps=100)

# or using scipy.optimize.fsolve
# import scipy.optimize as optimize
# f = sym.lambdify([n5, nTot], (expr5, expr6))
# root = optimize.fsolve(lambda x: f(*x), [[0.5, 1e7]], xtol=1e-16)
# print('f({}) = {}'.format(root, f(*root)))

# or using bisect with verify=False -- mentioned in the `sym.nsolve` docs
# root = sym.nsolve([rxn5, rxn6], [n5, nTot], [0.5, 0.5], solver='bisect', verify=False)

print(root)
# Matrix([[1.64738651107199e-8], [0.530353285064842]])

root = {n5:root[0], nTot:root[1]}
nsoln = {var:val.subs(root) for var,val in soln.items()}
nsoln.update(root)
for var in [n1, n2, n3, n4, n5, nTot]:
    print('{} = {}'.format(var, nsoln[var]))
for i, rxn in enumerate((rxn1, rxn2, rxn3, rxn4, rxn5, rxn6), 1):
    print('rxn{}: {}'.format(i, (rxn.lhs-rxn.rhs).subs(nsoln)))

这产生了这个解决方案：

n1 = 0.00000517060185532663
n2 = 0.0132120398187973
n3 = 0.0949735158684791
n4 = 0.422162542301846
n5 = 1.64738651107199E-8
nTot = 0.530353285064842

以及这些误差（各个方程左右两侧的差异）：

rxn1: -8.67361737988404E-19
rxn2: 0
rxn3: 131072.000000000
rxn4: -0.125000000000000
rxn5: -2.22044604925031E-16
rxn6: 5.55111512312578E-17

请注意，虽然 的绝对误差rxn3相当大，但相对误差 < 2e-16：

In [13]: 131072/1.07e21
Out[13]: 1.2249719626168225e-16

顺便说一句，文档sym.nsolve还说：

但是请注意，在根附近非常陡峭的函数可能会导致解决方案的验证失败。在这种情况下，您应该使用标志 verify=False 并独立验证解决方案。

和

如果根的边界已知并且使用二等分方法，则可以安全地跳过验证：

In [155]: sym.nsolve([rxn5, rxn6], [n5, nTot], [0.5, 0.5], verify=False, solver='bisect')
Out[155]: 
Matrix([
[2.77113579502186e-9],
[2.83553974806881e-6]])

但是这个解决方案并不像使用maxsteps=50（或更高）找到的解决方案一样好。

推荐阅读