Use NumPy for summation where j != i

Question

I have a function where a summation is calculated such that j != i as shown below:

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

I'm attempting to remove the for-loops using NumPy features but I'm not sure how to account for the j != i operation. Here's my attempt at using 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

An example of using the functions is below:

# 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

As expected, the NumPy version does not give the correct answer. Is it even possible to accomplish this with NumPy?

user_4975981 · Answer

Here is vectorize version of your function (it returns the same output as your function minus floating point error):

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

test:

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

Use NumPy for summation where j != i

1 Answers