我必須找到以下等式的零：這是一個(gè)狀態(tài)方程，如果您不確切知道 EoS 是什么，那也沒(méi)關(guān)系。根據(jù)上述等式的根，我計(jì)算（除其他事項(xiàng)外）不同壓力和溫度下氣態(tài)物質(zhì)Z的壓縮系數(shù)。通過(guò)這些解決方案，我可以繪制以壓力為橫坐標(biāo)、以Z s 為縱坐標(biāo)、以溫度為參數(shù)的曲線族。Beta、delta、eta 和 phi 是常數(shù)，pr 和 Tr 也是。在與 Newton-Raphson 方法（它與其他幾個(gè) EoS 一起工作得很好）我的頭腦沒(méi)有成功之后，我決定嘗試 Scipy 的root()功能。令我不滿(mǎn)的是，我得到了這張圖表：人們可以很容易地看出，這張鋸齒形圖表是完全有缺陷的。我應(yīng)該得到平滑的曲線。此外，Z通常介于 0.25 和 2.0 之間。因此，Z s 等于，比方說(shuō)，3 或以上是完全不合時(shí)宜的。然而Z < 2的曲線看起來(lái)還不錯(cuò)，盡管由于比例而高度壓縮。然后我嘗試了 Octave 的fzero()求解器，得到了這個(gè)：這正是我應(yīng)該得到的，因?yàn)槟切┦蔷哂姓_/預(yù)期形狀的曲線！我的問(wèn)題來(lái)了。顯然，Scipyroot()和 Octavefzero()是基于MINPACK 的相同算法混合體。盡管如此，結(jié)果顯然并不相同。你們有人知道為什么嗎？我繪制了 Octave（橫坐標(biāo)）獲得的 Zs 與 Scipy 獲得的 Zs 的曲線并得到了這個(gè)：底部暗示直線的點(diǎn)代表y = x，即 Octave 和 Scipy 在他們提出的解決方案中達(dá)成一致的點(diǎn)。其他觀點(diǎn)完全不一致，不幸的是，它們太多了，不能簡(jiǎn)單地忽略。我可能從現(xiàn)在開(kāi)始一直使用 Octave，因?yàn)樗梢怨ぷ?，但我想繼續(xù)使用 Python。你對(duì)此有何看法？有什么建議嗎？
2 回答
繁星coding
TA貢獻(xiàn)1797條經(jīng)驗(yàn) 獲得超4個(gè)贊
要事第一。您的兩個(gè)文件不相同，因此很難直接比較底層算法。我在這里附上一個(gè)八度和一個(gè) 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()
這證實(shí)了兩個(gè)系統(tǒng)的輸出實(shí)際上確實(shí)存在差異，部分原因在于所使用的底層算法的輸出，而不是因?yàn)槌绦驅(qū)嶋H上并不相同。然而，現(xiàn)在的比較并沒(méi)有那么糟糕：
至于“算法是一樣的”，其實(shí)不然。Octave 通常在源代碼中隱藏了更多的技術(shù)實(shí)現(xiàn)細(xì)節(jié)，所以這總是值得檢查的。特別是，在文件 fzero.m 中，就在文檔字符串之后，它提到了以下內(nèi)容：
這本質(zhì)上是Alefeld、Potra 和 Shi 的 ACM“算法 748：封閉連續(xù)函數(shù)的零點(diǎn)”，ACM Transactions on Mathematical Software，Vol。21，第 3 期，1995 年 9 月。
雖然工作流程應(yīng)該是一樣的，但算法的結(jié)構(gòu)已經(jīng)進(jìn)行了不平凡的改造；與作者順序調(diào)用構(gòu)建塊子程序的方法不同，我們?cè)谶@里實(shí)現(xiàn)了一個(gè) FSM 版本，每次迭代使用一次內(nèi)點(diǎn)確定和一次包圍，從而減少了臨時(shí)變量的數(shù)量并簡(jiǎn)化了算法結(jié)構(gòu)。此外，這種方法減少了對(duì)外部函數(shù)和錯(cuò)誤處理的需要。該算法也略有修改。
而根據(jù)help(root)：
備注
本節(jié)介紹可通過(guò)“方法”參數(shù)選擇的可用求解器。默認(rèn)方法是hybr。方法hybr使用 MINPACK [1] 中實(shí)現(xiàn)的 Powell 混合方法的修改。
參考文獻(xiàn)
[1]?More、Jorge J.、Burton S. Garbow 和 Kenneth E. Hillstrom。1980. MINPACK-1 用戶(hù)指南。
我嘗試了 中提到的幾個(gè)備選方案help(root)。一個(gè)df-sane似乎針對(duì)“標(biāo)量”值（即像“fzero”）進(jìn)行了優(yōu)化。事實(shí)上，雖然不如 Octave 的實(shí)現(xiàn)那么好，但這確實(shí)給出了稍微“理智”（雙關(guān)語(yǔ)意）的結(jié)果：
話(huà)雖如此，混合方法不會(huì)轉(zhuǎn)儲(chǔ)任何警告，但如果您使用其他一些替代方法，它們中的許多方法會(huì)告訴您您有很多有效的除以零、nans 和 infs 的地方你不應(yīng)該，這大概就是為什么你會(huì)得到如此奇怪的結(jié)果。所以，也許不是八度的算法本身“更好”，而是它在這個(gè)問(wèn)題中處理“被零除”的實(shí)例稍微更優(yōu)雅一些。
我不知道你的問(wèn)題的確切性質(zhì)，但可能是 python 方面的算法只是希望你提供條件良好的問(wèn)題。也許您在 zsbwr() 中的某些計(jì)算會(huì)導(dǎo)致被零除或不切實(shí)際的零等，您可以將其檢測(cè)并視為特殊情況？
梵蒂岡之花
TA貢獻(xiàn)1900條經(jīng)驗(yàn) 獲得超5個(gè)贊
（請(qǐng)將代碼修剪為最小的示例，該示例僅顯示找到不需要的根的根查找部分和參數(shù)。）
然后程序是手動(dòng)檢查方程以找到您想要的根的定位間隔并使用它。我通常使用brentq.
添加回答
舉報(bào)
0/150
提交
取消