八度的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文件中，它提到了以下内容：

这实际上是ACM "Algorithm 748: Enclosing Zeros of Continuous Functions" due to Alefeld,Potra and Shi,ACM Transactions on Mathematical Software,Vol. 21,No. 3,September 1995。

尽管工作流程应该相同，但是算法的结构已经进行了简单的转换；代替了作者按顺序调用构建块子程序的方法，我们在这里实现了一个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。