以下に示すように、 j != i となるように合計が計算される関数があります。
def mu_brokaw(mus, mws, xs):
n = len(mus)
mu_mix = 0.0
for i in range(n):
d = 0.0
for j in range(n):
if j != i:
mij = ((4 * mws[i] * mws[j]) / ((mws[i] + mws[j])**2))**0.25
num = (mws[i] / mws[j]) - (mws[i] / mws[j])**0.45
den = 2 * (1 + mws[i] / mws[j]) + ((1 + (mws[i] / mws[j])**0.45) / (1 + mij)) * mij
aij = mij * ((mws[j] / mws[i])**0.5) * (1 + num / den)
sij = 1.0
d = d + sij * aij * xs[j] / np.sqrt(mus[j])
mu_mix = mu_mix + (xs[i] * np.sqrt(mus[i])) / (xs[i] / np.sqrt(mus[i]) + d)
return mu_mix
NumPy 機能を使用して for ループを削除しようとしていますが、 j != i 操作をどのように考慮すればよいかわかりません。 NumPy を使用してみた例は次のとおりです。
def mu_brokaw2(mus, mws, xs):
mij = ((4 * np.outer(mws, mws)) / ((np.add.outer(mws, mws))**2))**0.25
num = np.divide.outer(mws, mws) - np.divide.outer(mws, mws)**0.45
den = 2 * (1 + np.divide.outer(mws, mws)) + (1 + np.divide.outer(mws, mws)**0.45) / (1 + mij) * mij
aij = mij * (np.divide.outer(mws, mws)**0.5) * (1 + num / den)
sij = 1.0
d = np.sum(sij * aij * xs / np.sqrt(mus))
mu_mix = np.sum((xs * np.sqrt(mus)) / (xs / np.sqrt(mus) + d))
return mu_mix
関数の使用例を以下に示します。
# dynamic gas viscosity in µP
mu_h2 = 179.75
mu_n2 = 363.87
# molecular weight in g/mol
mw_h2 = 2.016
mw_n2 = 28.014
# mole fraction
x_h2 = 0.85
x_n2 = 0.15
mu_mix = mu_brokaw([mu_h2, mu_n2], [mw_h2, mw_n2], [x_h2, x_n2])
print(f'mu_mix = {mu_mix:.4f}')
mu_mix2 = mu_brokaw2([mu_h2, mu_n2], [mw_h2, mw_n2], [x_h2, x_n2])
print(f'mu_mix2 = {mu_mix2:.4f}')
# mu_mix = 257.9015
# mu_mix2 = 41.1099
予想どおり、NumPy バージョンでは正しい答えが得られません。 NumPy を使用してこれを実現することは可能ですか?
1
サンプルの入力と出力を追加してください
– アクシャイ・セーガル
2020 年 9 月 3 日 2:44
d = sij * aij * xs / を置き換えてみてください。 np.sqrt(mus); d = d.sum(1)-np.diagonal(d)。加えてブロードキャストの修正も加えました。このコードをベクトル化してループを取り除くことができます。
– イーサン
2020 年 9 月 3 日 2:56
@Ehsan まだ正しい返答が得られません恨む。 NumPy 関数に別の問題がある可能性があります。もう一度見てみます。コメントの代わりに回答を送信していただけますか?
– かつら
2020 年 9 月 3 日 3:18
新しい i ループごとに d をゼロに設定したため、d は行に対して計算されます。これは、mu_mix で行を列計算から分離する必要があることを意味します。余談ですが、aij の外側の除算は、必要なものの逆になります。そこでtranspose np.divide.outer(mws, mws).T**0.5を使用してみてください。
– アガイ
2020 年 9 月 3 日4:20
------------------------
これは関数のベクトル化バージョンです (関数と同じ出力から浮動小数点エラーを除いたものを返します):
def mu_brokaw2(mus, mws, xs):
mij = ((4 * np.outer(mws, mws)) / ((np.add.outer(mws, mws))**2))**0.25
num = np.divide.outer(mws, mws) - np.divide.outer(mws, mws)**0.45
den = 2 * (1 + np.divide.outer(mws, mws)) + ((1 + np.divide.outer(mws, mws)**0.45) / (1 + mij)) * mij
aij = mij * (np.divide.outer(mws, mws).T**0.5) * (1 + num / den)
sij = 1.0
tij = sij * aij * xs / np.sqrt(mus)
d = tij.sum(1) - np.diagonal(tij)
mu_mix = np.sum((xs * np.sqrt(mus)) / (xs / np.sqrt(mus) + d))
return mu_mix
テスト:
mus = np.random.rand(100)
mws = np.random.rand(100)
xs = np.random.rand(100)
print(mu_brokaw(mus,mws,xs)-mu_brokaw2(mus,mws,xs))
#0.0
3
1 / を使用する理由はありますか? np.divide.outer(mws, mws) を np.divide.outer(mws, mws).T であるアプローチと比較しますか?どちらも同じ結果になります。
– かつら
2020 年 9 月 3 日 16:33
1
@かつらの配合に関してはノーです。パフォーマンス的にはそうです。ここでの .T は単に np.divide.outer(mws, mws) へのビューですが、1/np.divide.outer(mws, mws) は追加の除算ステップを実行します。問題が解決する場合は、同意して質問を閉じてください。ありがとう
– イーサン
2020 年 9 月 3 日 20:31
また、私があなただったら、パフォーマンスが重要な場合に備えて、外側の計算をすべて保存し、数式で使用するでしょう。方程式に二重計算がたくさんあるようです。
– イーサン
2020 年 9 月 3 日 20:34
