我必須找到以下等式的零：這是一個狀態方程，如果您不確切知道 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。你對此有何看法？有什么建議嗎？
2 回答
繁星coding
TA貢獻1797條經驗 獲得超4個贊
要事第一。您的兩個文件不相同，因此很難直接比較底層算法。我在這里附上一個八度和一個 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() 中的某些計算會導致被零除或不切實際的零等，您可以將其檢測并視為特殊情況？
梵蒂岡之花
TA貢獻1900條經驗 獲得超5個贊
（請將代碼修剪為最小的示例，該示例僅顯示找到不需要的根的根查找部分和參數。）
然后程序是手動檢查方程以找到您想要的根的定位間隔并使用它。我通常使用brentq.
添加回答
舉報
0/150
提交
取消