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

Question

我需要找到以下方程的零点：

这是一个状态方程，如果您不知道什么是EoS也没关系。通过上述方程的根，我计算出了气态物质在不同压力和温度下的压缩因子Z（以及其他一些内容）。利用这些解，我可以绘制一族曲线，其中压力为横坐标，Z为纵坐标，温度为参数。Beta、delta、eta和phi是常数，pr和Tr也是。

在试过牛顿-拉夫逊方法无果后（该方法在其他几种EoS中效果良好），我决定尝试Scipy的root()函数。但令我不满的是，我得到了这张图表：

作为一个很容易察觉到的问题，这个锯齿状的图表是完全有缺陷的。我本应该得到平滑的曲线。此外，Z通常在0.25到2.0之间变化。因此，例如Z大于或等于3就完全不准确了。然而，Z小于2的曲线看起来还可以，尽管由于比例尺的缘故高度压缩。

然后我试了Octave的fzero()求解器，并得到了这个结果：

这正是我应该得到的，因为这些曲线具有正确/预期形状！

接下来是我的问题。显然，Scipy的root()和Octave的fzero()都基于MINPACK的同一算法hybrid。然而，结果明显不同。你们中的任何人知道原因吗？

我绘制了一个由Octave获得的Zs曲线（横坐标）与使用Scipy获得的曲线对比：

底部提示一条直线的点表示y = x，即Octave和Scipy在它们提供的解决方案中达成一致的点。其他点完全不同意，而且不幸的是，它们太多了，不能简单地忽略。

我可能会一直使用Octave，因为它有效，但我想继续使用Python。

你对此有什么看法？有什么建议吗？

附注：这是原始的Python代码。它产生了第一个图表。

import numpy
from scipy.optimize import root
import matplotlib.pyplot as plt

def fx(x, beta, delta, eta, phi, pr_, Tr_):
    tmp = phi*x**2
    etmp = numpy.exp(-tmp)
    f = x*(1.0 + beta*x + delta*x**4 + eta*x**2*(1.0 + tmp)*etmp) - pr_/Tr_
    return f

def zsbwr(pr_, Tr_, pc_, Tc_, zc_, w_, MW_, phase=0):

    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))

    raiz = root(fx,y0,args=(beta, delta, eta, phi, pr_, Tr_),method='hybr',tol=1.0e-06)

    return pr_/raiz.x[0]/Tr_


if __name__ == "__main__":

    Tc = 304.13
    pc = 73.773
    omega = 0.22394
    zc = 0.2746
    MW = 44.01

    Tr = numpy.array([0.8, 0.93793103])
    pr = numpy.linspace(0.5, 14.5, 25)

    zfactor = numpy.zeros((2, 25))

    for redT in Tr:
        j = numpy.where(Tr == redT)[0][0]
        for redp in pr:
            indp = numpy.where(pr == redp)[0][0]
            zfactor[j][indp] = zsbwr(redp, redT, pc, Tc, zc, omega, MW, 0)

    for key, value in enumerate(zfactor):
        plt.plot(pr, value, '.-', linewidth=1, color='#ef082a')

    plt.figure(1, figsize=(7, 6))
    plt.xlabel('$p_R$', fontsize=16)
    plt.ylabel('$Z$', fontsize=16)
    plt.grid(color='#aaaaaa', linestyle='--', linewidth=1)
    plt.show()

现在是Octave脚本：

function SoaveBenedictWebbRubin

    format long;

    nTr = 11;
    npr = 43;

    ic = 1;

    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];

    Tc = comp(ic,1);
    pc = comp(ic,2);
    w = comp(ic,3);
    zc = comp(ic,4);
    MW = comp(ic,5);

    Tr = linspace(0.8, 2.8, nTr);
    pr = linspace(0.2, 7.2, npr);

    figure(1, 'position',[300,150,600,500])

    for i=1:size(Tr, 2)
        icont = 1;
        zval = zeros(1, npr);
        for j=1:size(pr, 2)
            [Z, phi, density] = SBWR(Tr(i), pr(j), Tc, pc, zc, w, MW, 0);
            zval(icont) = Z;
            icont = icont + 1;
        endfor
        plot(pr,zval,'o','markerfacecolor','white','linestyle','-','markersize',3);
        hold on;
    endfor

    str = strcat("Soave-Benedict-Webb-Rubin para","\t",nome(ic));
    xlabel("p_r",'fontsize',15);
    ylabel("Z",'fontsize',15);
    title(str,'fontsize',12);
end

function [Z,phi,density] = SBWR(Tr, pr, Tc, pc, Zc, w, MW, phase)
    R = 8.3144E-5; % universal gas constant (bar·m3/(mol·K))

    % 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;
    ff = 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));
        end
    end

    fun = @(y)y*(1.0 + beta*y + delta*y**4 + eta*y**2*(1.0 + ff*y**2)*exp(-ff*y**2)) - pr/Tr;

    options = optimset('TolX',1.0e-06);
    yi = fzero(fun,y0,options);

    Z = pr/yi/Tr;
    density = yi*pc*MW/(1000.0*R*Tc);
    phi = exp(Z - 1.0 - log(Z) + beta*yi + 0.25*delta*yi**4 - eta/ff*(exp(-ff*yi**2)*(1.0 + 0.5*ff*yi**2) - 1.0));
end

- Carlos Gouveia

---

展示一下Python代码。可能有几个根需要收敛。 - ev-br

@ev-br 当然有多个根。让我惊讶的是Octave可以从提供的初始猜测中找到正确的根。而且这种方法很可能与scipy.optimize.root()使用的方法相同。 - Carlos Gouveia

@ev-br 另外，函数 zsbwr() 返回的 phi 不是上面 f(x) 方程中显示的希腊字母 phi。后者应该是 ff。我的代码有点混乱。对此我感到抱歉。 - Carlos Gouveia

@TasosPapastylianou，完成。 - Carlos Gouveia

我注意到你的参数（pr 和 Tr）在两个版本中是不同的。也许这是个问题？（我也没有在 Python 代码中看到普适气体常数，但你没有绘制涉及那里的变量，所以我认为这与之无关） - Tasos Papastylianou

Answer 1

首先，你的两个文件并不相等，因此直接比较底层算法是困难的。我在这里附上一个Octave版本和一个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 Transactions on Mathematical Software，Vol. 21，No. 3，September 1995上发表的ACM“Algorithm 748: Enclosing Zeros of Continuous Functions”。”

尽管工作流程应该是相同的，但算法的结构已经被非常复杂地转换了；我们在这里使用一种FSM版本，使用一个内部点确定和一个括号来迭代，从而减少了临时变量的数量并简化了算法结构，而不是作者顺序调用构建块子程序的方法。此外，这种方法还减少了对外部函数和错误处理的需求。该算法也稍作修改。

而根据help(root)：

注
本节介绍了可以通过“method”参数选择的可用求解器。默认方法是hybr。

方法hybr使用了Powell混合方法的修改版本，该方法在MINPACK [1]中实现。

参考文献
[1] More, Jorge J., Burton S. Garbow, and Kenneth E. Hillstrom. 1980. MINPACK-1用户指南。

我尝试了一些在help(root)中提到的替代方案。 df-sane似乎针对“标量”值进行了优化（即像“fzero”一样）。实际上，虽然不如Octave的实现好，但这确实给出了一个稍微“更加理智”的结果（双关语）：

话虽如此，混合方法不会抛出任何警告，但如果您使用其他一些替代方法，则其中许多方法将通知您在不应该的地方存在大量有效的零除法、无穷大和未定义数，这可能是为什么您会得到如此奇怪的结果。因此，也许Octave算法本身并不是“更好”，而是在这个问题中稍微更优雅地处理了“除以零”的情况。

我不知道您的问题的确切性质，但可能是Python端的算法只是希望您提供条件良好的问题。也许zsbwr()中的某些计算会导致除零或不现实的零等情况，您可以检测并将其视为特殊情况处理？

- Tasos Papastylianou

---

首先，非常感谢您详细和友善的回复。确实，您的第二张图比我的好得多，尽管蓝色曲线完全是弹道式的，并显示出不现实的行为（通常，我不会期望Z>=2，尽管这是一个理论上可能的值）。我敢打赌这个蓝色曲线是Tr = 0.8的曲线。我认为（但无法证明）这与数值异常有关，并且Octave更优雅地处理它们。在我看来，这是一个实现缺陷，因为返回的值虽然异常高，但并没有像O(10 ** 10)之类的疯狂值。 - Carlos Gouveia

（续）因此，Z = 7很奇怪，但对于不熟悉的人来说，它看起来像是一个有教养的根。我认为我们永远不会找到这个谜题的最终答案。我能做的就是在下一次解决类似复杂方程的机会出现时，始终关注这样的解决方案。这是一个复杂的问题。我有机会研究和使用的EOS比这个SBWR简单得多。再次感谢您的回答；如果可以的话，我会投两票。;-) - Carlos Gouveia

@CarlosGouveia 嗯，如果您认为它回答了您的问题，那么您可以接受这个答案:p（我想它回答了“它们不是相同的算法”部分，基本上的答案是“不是”）。至于图表是否有意义，我想要做的合理的事情（如果您还没有）就是在各种pr值的x域上简单地绘制EOS方程（以及固定的Tr，例如您认为有问题的那个），看看您得到的图表是否有意义或者从视觉上解释为什么“根查找”算法在那种情况下难以找到零点！ - Tasos Papastylianou

我必须“正式”接受一个答案吗？比如，关闭这个话题？我该怎么做？我很乐意接受你的答案，Tasos。 - Carlos Gouveia

@CarlosGouveia 如果没有一个答案适合你，你不必“强制”接受一个答案，但如果有一个答案回答了你的问题，那么最好接受它，主要是为了让问题在问题列表中显示为已关闭，这有助于其他搜索类似问题的人。你可以通过点击答案旁边的“打勾”（就在投票箭头下方）来“接受”一个答案。请注意，接受一个答案并不会在某种意义上“关闭”一个话题。你可以随时改变主意，选择另一个答案作为被接受的答案。 - Tasos Papastylianou

Answer 2

（请将代码裁剪为仅显示寻找不需要的根和参数的最小示例）

然后，手动检查方程以找到所需根的定位间隔并使用它。我通常使用 brentq。

- ev-br

---

已完成。我将其尽可能地简化了。我还发现问题发生在Tr = 0.8和Tr = 0.93793103时，因此修剪后的代码仅绘制这两个降温度数的图表。出于某种难以理解的原因，SciPy无法很好地处理这两个Trs。然而，Octave可以很好地处理它们。 - Carlos Gouveia