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

要事第一。您的两个文件不相同，因此很难直接比较底层算法。我在这里附上一个八度和一个 python 版本，它们可以直接逐行比较，可以并排比较。%%% File: SoaveBenedictWebbRubin.m:% No package imports necessaryfunction 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 % mainfunction 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 % zSBWRfunction 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.pyimport numpy; from scipy.optimize import root; import matplotlib.pyplot as pltdef 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 maindef 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 zsbwrdef 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 fxif __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() 中的某些计算会导致被零除或不切实际的零等，您可以将其检测并视为特殊情况？

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

2回答