如何解决八度的fzero和Scipy的root函数不会产生相同的结果
我必须找到以下等式的零:
这是一个状态方程,如果您不确切知道什么是EoS,那么它并不重要。通过上述方程式的根,我可以计算(其中包括)不同压力和温度下气态物质 Z 的压缩系数。通过这些解决方案,我可以绘制压力为横坐标, Z s为纵坐标,温度为参数的曲线系列。 Beta,delta,eta和phi以及pr和Tr是常数。
在对Newton-Raphson方法(与其他几种EoS很好地配合)失败之后,我决定尝试Scipy的role
函数。令我不满的是,我获得了这张图表:
正如一个人容易理解的那样,这张锯齿形的图表是完全有缺陷的。我应该得到平滑的曲线。另外, Z 的范围通常在0.25到2.0之间。因此,等于或大于3的 Z 完全不符合要求。然而, Z
然后我尝试使用Octave的select `person`.*
from `people` `person`
where `person`.`id` in (
select `credit`.`person_id`
from `roles` `role`
join`credits` `credit`
on `role`.`id` = `credit`.`role_id`
where `role` like "Director"
) and `person`.`id` in (
select `credit`.`person_id`
from `roles` `role`
join`credits` `credit`
on `role`.`id` = `credit`.`role_id`
where `role` like "Actor"
);
求解器,并获得了此信息:
这正是我应该得到的,因为那是具有正确/预期形状的曲线!
这是我的问题。显然,Scipy的root()
和Octave的fzero()
基于MINPACK的同一算法 hybrid 。尽管如此,结果显然还是不一样。你们当中有人知道为什么吗?
我绘制了由Octave(横坐标)获得的Zs与由Scipy获得的Zs的曲线,并得出了这一点:
底部提示直线的点代表root()
,即Octave和Scipy在他们提出的解决方案中达成一致的点。其他要点完全不同,不幸的是,它们太多了,不能被简单地忽略。
从现在开始,我可能会一直使用Octave,但我想继续使用Python。
您对此有何看法?有什么建议吗?
PS:这是原始的Python代码。它会生成此处显示的第一个图表。
fzero()
现在是八度脚本:
y = x
解决方法
第一件事。您的两个文件不相等,因此很难直接比较基础算法。我在这里附上一个八度音阶和一个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",15 );
title( ["Soave-Benedict-Webb-Rubin para\t",nome(ic)],12 );
endfunction % main
function Z = zSBWR( Tr,pr,Zc,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,Tr),y0,optimset( 'TolX',1.0e-06 ) );
Z = pr / yi / Tr;
endfunction % zSBWR
function Out = fx( y,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,44.0100 ],[ 126.19,28.0134 ],[ 647.14,18.0153 ],[ 190.56,16.0430 ],[ 305.33,30.0700 ],[ 369.83,44.0970 ] ] )
nTr = 11; Tr = numpy.linspace( 0.8,nTr )
npr = 43; pr = numpy.linspace( 0.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],0)
# endfor
# endfor
for i in range(nTr):
plt.plot(pr,zvalues[i,:],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,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,args = (beta,method = 'hybr',tol = 1.0e-06 ).x
return pr / yi / Tr
# endfunction zsbwr
def fx(y,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文件中,它提到了以下内容:
尽管工作流程应该相同,但是算法的结构已经进行了简单的转换;代替了作者按顺序调用构建块子程序的方法,我们在这里实现了一个FSM版本,该版本使用一个内部点确定和每个迭代一个括号,从而减少了临时变量的数量并简化了算法结构。此外,这种方法减少了对外部功能和错误处理的需求。该算法也做了一些修改。
根据help(root)
:
注释
本节介绍可通过“方法”参数选择的可用求解器。默认方法是 hybr 。方法 hybr 对Powell混合方法进行了修改, 在MINPACK [1]中实现。
参考
[1] More,Jorge J.,Burton S. Garbow,and Kenneth E. Hillstrom. 1980. User Guide for MINPACK-1.
我尝试了help(root)
中提到的几种替代方法。 df-sane
似乎已针对“标量”值(即“ fzero”)进行了优化。确实,虽然不如八度音阶的实现好,但这确实会产生稍微“更刺眼”(双关语意)的结果:
话虽如此,混合方法不会转储任何警告,但是如果您使用其他一些替代方法,则许多替代方法都将告知您有效除以零,nans和infs,在某些地方您不应该这样做,这大概就是为什么您得到如此奇怪的结果的原因。因此,也许并不是八度音阶算法本身就“更好”了,但它在这个问题上的处理方式为“被零除”的实例稍微更优雅一些。
我不知道您问题的确切性质,但可能是python方面的算法只是希望您提供条件良好的问题。也许您在zsbwr()中进行的某些计算会导致除以零的情况或不切实际的零等,您可以将其检测并视为特殊情况?
,(请将代码修剪为一个最小的示例,该示例仅显示找到不需要的根的根查找部分和参数。)
然后,该过程是手动检查方程式,以找到所需根的本地化间隔并使用它。我通常使用brentq
。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。