我必须找到以下等式的零：这是一个状态方程，如果您不确切知道 EoS 是什么，那也没关系。根据上述等式的根，我计算（除其他事项外）不同压力和温度下气态物质Z的压缩系数。通过这些解决方案，我可以绘制以压力为横坐标、以Z s 为纵坐标、以温度为参数的曲线族。Beta、delta、eta 和 phi 是常数，pr 和 Tr 也是。在与 Newton-Raphson 方法（它与其他几个 EoS 一起工作得很好）我的头脑没有成功之后，我决定尝试 Scipy 的root()功能。令我不满的是，我得到了这张图表：人们可以很容易地看出，这张锯齿形图表是完全有缺陷的。我应该得到平滑的曲线。此外，Z通常介于 0.25 和 2.0 之间。因此，Z s 等于，比方说，3 或以上是完全不合时宜的。然而Z < 2的曲线看起来还不错，尽管由于比例而高度压缩。然后我尝试了 Octave 的fzero()求解器，得到了这个：这正是我应该得到的，因为那些是具有正确/预期形状的曲线！我的问题来了。显然，Scipyroot()和 Octavefzero()是基于MINPACK 的相同算法混合体。尽管如此，结果显然并不相同。你们有人知道为什么吗？我绘制了 Octave（横坐标）获得的 Zs 与 Scipy 获得的 Zs 的曲线并得到了这个：底部暗示直线的点代表y = x，即 Octave 和 Scipy 在他们提出的解决方案中达成一致的点。其他观点完全不一致，不幸的是，它们太多了，不能简单地忽略。我可能从现在开始一直使用 Octave，因为它可以工作，但我想继续使用 Python。你对此有何看法？有什么建议吗？
2 回答
繁星coding
TA贡献1797条经验 获得超4个赞
要事第一。您的两个文件不相同，因此很难直接比较底层算法。我在这里附上一个八度和一个 python 版本，它们可以直接逐行比较，可以并排比较。
%%% File: SoaveBenedictWebbRubin.m:
% No package imports necessary
function SoaveBenedictWebbRubin()
nome = {"CO2"; "N2"; "H2O"; "CH4"; "C2H6"; "C3H8"};
comp = [ 304.13, 73.773, 0.22394, 0.2746, 44.0100;
126.19, 33.958, 0.03700, 0.2894, 28.0134;
647.14, 220.640, 0.34430, 0.2294, 18.0153;
190.56, 45.992, 0.01100, 0.2863, 16.0430;
305.33, 48.718, 0.09930, 0.2776, 30.0700;
369.83, 42.477, 0.15240, 0.2769, 44.0970 ];
nTr = 11; Tr = linspace( 0.8, 2.8, nTr );
npr = 43; pr = linspace( 0.2, 7.2, npr );
ic = 1;
Tc = comp(ic, 1);
pc = comp(ic, 2);
w = comp(ic, 3);
zc = comp(ic, 4);
MW = comp(ic, 5);
figure(1, 'position',[300,150,600,500])
zvalues = zeros( nTr, npr );
for i = 1 : nTr
for j = 1 : npr
zvalues(i,j) = zSBWR( Tr(i), pr(j), Tc, pc, zc, w, MW, 0 );
endfor
endfor
hold on
for i = 1 : nTr
plot( pr, zvalues(i,:), 'o-', 'markerfacecolor', 'white', 'markersize', 3);
endfor
hold off
xlabel( "p_r", 'fontsize', 15 );
ylabel( "Z" , 'fontsize', 15 );
title( ["Soave-Benedict-Webb-Rubin para\t", nome(ic)], 'fontsize', 12 );
endfunction % main
function Z = zSBWR( Tr, pr, Tc, pc, Zc, w, MW, phase )
% Definition of parameters
d1 = 0.4912 + 0.6478 * w;
d2 = 0.3 + 0.3619 * w;
e1 = 0.0841 + 0.1318 * w + 0.0018 * w ** 2;
e2 = 0.075 + 0.2408 * w - 0.014 * w ** 2;
e3 = -0.0065 + 0.1798 * w - 0.0078 * w ** 2;
f = 0.77;
ee = (2.0 - 5.0 * Zc) * exp( f ) / (1.0 + f + 3.0 * f ** 2 - 2.0 * f ** 3 );
d = (1.0 - 2.0 * Zc - ee * (1.0 + f - 2.0 * f ** 2) * exp( -f ) ) / 3.0;
b = Zc - 1.0 - d - ee * (1.0 + f) * exp( -f );
bc = b * Zc;
dc = d * Zc ** 4;
ec = ee * Zc ** 2;
phi = f * Zc ** 2;
beta = bc + 0.422 * (1.0 - 1.0 / Tr ** 1.6) + 0.234 * w * (1.0 - 1.0 / Tr ** 3);
delta = dc * (1.0 + d1 * (1.0 / Tr - 1.0) + d2 * (1.0 / Tr - 1.0) ** 2);
eta = ec + e1 * (1.0 / Tr - 1.0) + e2 * (1.0 / Tr - 1.0) ** 2 + e3 * (1.0 / Tr - 1.0) ** 3;
if Tr > 1
y0 = pr / Tr / (1.0 + beta * pr / Tr);
else
if phase == 0
y0 = pr / Tr / (1.0 + beta * pr / Tr);
else
y0 = 1.0 / Zc ** (1.0 + (1.0 - Tr) ** (2.0 / 7.0) );
endif
endif
yi = fzero( @(y) fx(y, beta, delta, eta, phi, pr, Tr), y0, optimset( 'TolX', 1.0e-06 ) );
Z = pr / yi / Tr;
endfunction % zSBWR
function Out = fx( y, beta, delta, eta, phi, pr, Tr )
Out = y * (1.0 + beta * y + delta * y ** 4 + eta * y ** 2 * (1.0 + phi * y ** 2) * exp( -phi * y ** 2 ) ) - pr / Tr;
endfunction
### File: SoaveBenedictWebbRubin.py
import numpy; from scipy.optimize import root; import matplotlib.pyplot as plt
def SoaveBenedictWebbRubin():
nome = ["CO2", "N2", "H2O", "CH4", "C2H6", "C3H8"]
comp = numpy.array( [ [ 304.13, 73.773, 0.22394, 0.2746, 44.0100 ],
[ 126.19, 33.958, 0.03700, 0.2894, 28.0134 ],
[ 647.14, 220.640, 0.34430, 0.2294, 18.0153 ],
[ 190.56, 45.992, 0.01100, 0.2863, 16.0430 ],
[ 305.33, 48.718, 0.09930, 0.2776, 30.0700 ],
[ 369.83, 42.477, 0.15240, 0.2769, 44.0970 ] ] )
nTr = 11; Tr = numpy.linspace( 0.8, 2.8, nTr )
npr = 43; pr = numpy.linspace( 0.2, 7.2, npr )
ic = 0
Tc = comp[ic, 0]
pc = comp[ic, 1]
w = comp[ic, 2]
zc = comp[ic, 3]
MW = comp[ic, 4]
plt.figure(1, figsize=(7, 6))
zvalues = numpy.zeros( (nTr, npr) )
for i in range( nTr ):
for j in range( npr ):
zvalues[i,j] = zsbwr( Tr[i], pr[j], pc, Tc, zc, w, MW, 0)
# endfor
# endfor
for i in range(nTr):
plt.plot(pr, zvalues[i, :], 'o-', markerfacecolor='white', markersize=3 )
plt.xlabel( '$p_r$', fontsize = 15 )
plt.ylabel( '$Z$' , fontsize = 15 )
plt.title( "Soave-Benedict-Webb-Rubin para\t" + nome[ic], fontsize = 12 );
plt.show()
# end function main
def zsbwr( Tr, pr, pc, Tc, zc, w, MW, phase=0):
# Definition of parameters
d1 = 0.4912 + 0.6478 * w
d2 = 0.3000 + 0.3619 * w
e1 = 0.0841 + 0.1318 * w + 0.0018 * w ** 2
e2 = 0.075 + 0.2408 * w - 0.014 * w ** 2
e3 = -0.0065 + 0.1798 * w - 0.0078 * w ** 2
f = 0.770
ee = (2.0 - 5.0 * zc) * numpy.exp( f ) / (1.0 + f + 3.0 * f ** 2 - 2 * f ** 3)
d = (1.0 - 2.0 * zc - ee * (1.0 + f - 2.0 * f ** 2) * numpy.exp( -f )) / 3.0
b = zc - 1.0 - d - ee * (1.0 + f) * numpy.exp( -f )
bc = b * zc
dc = d * zc ** 4
ec = ee * zc ** 2
phi = f * zc ** 2
beta = bc + 0.422 * (1.0 - 1.0 / Tr ** 1.6) + 0.234 * w * (1.0 - 1.0 / Tr ** 3)
delta = dc * (1.0 + d1 * (1.0 / Tr - 1.0) + d2 * (1.0 / Tr - 1.0) ** 2)
eta = ec + e1 * (1.0 / Tr - 1.0) + e2 * (1.0 / Tr - 1.0) ** 2 + e3 * (1.0 / Tr - 1.0) ** 3
if Tr > 1:
y0 = pr / Tr / (1.0 + beta * pr / Tr)
else:
if phase == 0:
y0 = pr / Tr / (1.0 + beta * pr / Tr)
else:
y0 = 1.0 / zc ** (1.0 + (1.0 - Tr) ** (2.0 / 7.0))
# endif
# endif
yi = root( fx, y0, args = (beta, delta, eta, phi, pr, Tr), method = 'hybr', tol = 1.0e-06 ).x
return pr / yi / Tr
# endfunction zsbwr
def fx(y, beta, delta, eta, phi, pr, Tr):
return y*(1.0 + beta*y + delta*y**4 + eta*y**2*(1.0 + phi*y**2)*numpy.exp(-phi*y**2)) - pr/Tr
# endfunction fx
if __name__ == "__main__": SoaveBenedictWebbRubin()
这证实了两个系统的输出实际上确实存在差异，部分原因在于所使用的底层算法的输出，而不是因为程序实际上并不相同。然而，现在的比较并没有那么糟糕：
至于“算法是一样的”，其实不然。Octave 通常在源代码中隐藏了更多的技术实现细节，所以这总是值得检查的。特别是，在文件 fzero.m 中，就在文档字符串之后，它提到了以下内容：
这本质上是Alefeld、Potra 和 Shi 的 ACM“算法 748：封闭连续函数的零点”，ACM Transactions on Mathematical Software，Vol。21，第 3 期，1995 年 9 月。
虽然工作流程应该是一样的，但算法的结构已经进行了不平凡的改造；与作者顺序调用构建块子程序的方法不同，我们在这里实现了一个 FSM 版本，每次迭代使用一次内点确定和一次包围，从而减少了临时变量的数量并简化了算法结构。此外，这种方法减少了对外部函数和错误处理的需要。该算法也略有修改。
而根据help(root)：
备注
本节介绍可通过“方法”参数选择的可用求解器。默认方法是hybr。方法hybr使用 MINPACK [1] 中实现的 Powell 混合方法的修改。
参考文献
[1] More、Jorge J.、Burton S. Garbow 和 Kenneth E. Hillstrom。1980. MINPACK-1 用户指南。
我尝试了 中提到的几个备选方案help(root)。一个df-sane似乎针对“标量”值（即像“fzero”）进行了优化。事实上，虽然不如 Octave 的实现那么好，但这确实给出了稍微“理智”（双关语意）的结果：
话虽如此，混合方法不会转储任何警告，但如果您使用其他一些替代方法，它们中的许多方法会告诉您您有很多有效的除以零、nans 和 infs 的地方你不应该，这大概就是为什么你会得到如此奇怪的结果。所以，也许不是八度的算法本身“更好”，而是它在这个问题中处理“被零除”的实例稍微更优雅一些。
我不知道你的问题的确切性质，但可能是 python 方面的算法只是希望你提供条件良好的问题。也许您在 zsbwr() 中的某些计算会导致被零除或不切实际的零等，您可以将其检测并视为特殊情况？
梵蒂冈之花
TA贡献1900条经验 获得超5个赞
（请将代码修剪为最小的示例，该示例仅显示找到不需要的根的根查找部分和参数。）
然后程序是手动检查方程以找到您想要的根的定位间隔并使用它。我通常使用brentq.
添加回答
举报
0/150
提交
取消