Octave 的 fzero() 和 Scipy 的 root() 函数不会产生相同的结果

Question

我必须找到以下等式的零：这是一个状态方程，如果您不确切知道 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。你对此有何看法？有什么建议吗？

Answer 1

要事第一。您的两个文件不相同，因此很难直接比较底层算法。我在这里附上一个八度和一个 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() 中的某些计算会导致被零除或不切实际的零等，您可以将其检测并视为特殊情况？

Answer 2

（请将代码修剪为最小的示例，该示例仅显示找到不需要的根的根查找部分和参数。）

然后程序是手动检查方程以找到您想要的根的定位间隔并使用它。我通常使用brentq.