时间戳

首页 > 解决方案 > Octave 的 fzero() 和 Scipy 的 root() 函数不会产生相同的结果

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

问题描述

我必须找到以下等式的零：

在此处输入图像描述

这是一个状态方程，如果您不确切知道 EoS 是什么，这并不重要。根据上述等式的根，我计算（除其他外）气态物质Z在不同压力和温度下的压缩系数。通过这些解决方案，我可以绘制以压力为横坐标、 Z为纵坐标、温度为参数的曲线族。Beta、delta、eta 和 phi 是常数，pr 和 Tr 也是如此。

在用 Newton-Raphson 方法（它与其他几个 EoS 工作得很好）失败后，我决定尝试 Scipy 的root()功能。令我不满的是，我得到了这张图表：

在此处输入图像描述

可以很容易地看出，这张锯齿状的图表是完全有缺陷的。我应该得到平滑的曲线。此外，Z通常介于 0.25 和 2.0 之间。因此，Z s 等于或大于 3 是完全不符合要求的。然而Z < 2的曲线看起来还不错，尽管由于比例而被高度压缩。

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

在此处输入图像描述

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

我的问题来了。显然 Scipy'sroot()和 Octave'sfzero()基于 MINPACK 的相同算法混合。尽管如此，结果显然并不相同。你们中有人知道为什么吗？

我绘制了一条由 Octave（横坐标）获得的 Z 曲线与使用 Scipy 获得的 Z 曲线，并得到了这个：

在此处输入图像描述

底部暗示直线y = x的点代表 ，即 Octave 和 Scipy 在他们提出的解决方案中同意的点。其他观点完全不同，不幸的是，它们太多了，不能简单地忽略。

从现在开始我可能会一直使用 Octave，因为它可以工作，但我想继续使用 Python。

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

PS：这是原始的 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

标签： pythonscipyoctave

解决方案

第一件事。您的两个文件不等价，因此很难直接比较底层算法。我在这里附上一个八度音阶和一个 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 中，就在 docstring 之后，它提到了以下内容：

这实质上是Alefeld、Potra 和 Shi 的 ACM“算法 748：封闭连续函数的零点”，ACM Transactions on Mathematical Software，Vol。1995 年 9 月 21 日第 3 期。

尽管工作流程应该是相同的，但算法的结构已经发生了重大变化；代替作者顺序调用构建块子程序的方法，我们在这里实现了一个 FSM 版本，每次迭代使用一个内点确定和一个括号，从而减少临时变量的数量并简化算法结构。此外，这种方法减少了对外部功能和错误处理的需求。该算法也略有修改。

而根据help(root)：

注释
本节描述了可以通过 'method' 参数选择的可用求解器。默认方法是hybr。

方法hybr使用了在 MinPACK [1] 中实现的 Powell 混合方法的修改。

参考文献
[1]更多、Jorge J.、Burton S. Garbow 和 Kenneth E. Hillstrom。1980. MINPACK-1 用户指南。

我尝试了help(root). 这个df-sane似乎针对“标量”值（即像“fzero”）进行了优化。事实上，虽然不如 octave 的实现好，但这确实给出了一个稍微“理智”（双关语）的结果：

说了这么多，混合方法不会转储任何警告，但是如果您使用其他一些替代方法，它们中的许多会告诉您，您有很多有效的除以零、nans 和 infs 的地方是你不应该，这大概就是为什么你会得到如此奇怪的结果。所以，也许不是 octave 的算法本身“更好”，而是它在这个问题中处理“被零除”的实例稍微优雅一些​​。

我不知道你的问题的确切性质，但可能是 python 方面的算法只是希望你给它提供条件良好的问题。也许您在 zsbwr() 中的某些计算导致除以零次或不切实际的零等，您可以将其检测并视为特殊情况？

推荐阅读