numpy 函数给出不正确的结果 - 手动检查 excel

Question

我正在用 numpy 编写一些用于岩石物理建模的函数，我注意到我的一个函数给出了错误的结果。该函数是我对 Hertz-Mindlin 球体建模的实现：

Summary of the Hertz-Mindlin model

这是我目前的功能：

# Hertz-Mindlin sphere pack model: 

import numpy as np 

def hertzmindlin(K0, G0, PHIC, P, f=1.0):
'''
Hertz-Mindlin sphere-pack model, adapted from:
'Dvorkin, J. and Nur, A., 1996. Elasticity of high-porosity sandstones: 
Theory for two North Sea data sets. Geophysics, 61(5), pp.1363-1370."

Arguments:
K0 = Bulk modulus of mineral in GPa
G0 = Shear modulus of mineral in GPa
PHIC = Critical porosity for mineral-fluid mixture. Calculate using Dvorkin-Nuir (1995) or use literature
P = Confining pressure in GPa
f = Shear modulus correction factor. Default = 1

Results:
V0 = Theoretical poissons ratio of mineral
n = Coordination number of sphere-pack, calculated from Murphy's (1982) empirical relation
K_HM = Hertz-Mindlin effective dry Bulk modulus at pressure, P, in GPa
G_HM = Hertz-Mindlin effective dry Shear modulus at pressure, P, in GPa

'''
V0 = (3*K0-2*G0)/(6*K0+2*G0) # Calculated theoretical poissons ratio of bulk rock
n = 20-(34*PHIC)+(14*(PHIC**2)) # Coordination number at critical porosity (Murphy 1982)
K_HM = (P*(n**2*(1-PHIC)**2*G0**2) / (18*np.pi**2*(1-V0)**2))**(1/3)
G_HM = ((2+3*f-V0*(1+3*f))/(5*(2-V0))) * ((P*(3*n**2*(1-PHIC)**2*G0**2)/(2*np.pi**2*(1-V0)**2)))**(1/3)
return K_HM, G_HM

问题是当我运行这个函数用于输入时：

K, G, = 36, 45

PHIC = 0.4

P = 0.001

我得到的结果是 K_HM = 1.0，G_HM = 0.49009009009009

手算和excel计算值显示这是不正确的，我应该输出K_HM = 0.763265313，G_HM = 1.081083984

基于输入 K、G 的输出 G 应该大于 K（目前较小）这一事实，我相当确定函数中出现了问题

如有任何帮助，我们将不胜感激！我可以在 excel 中执行此操作，但理想情况下希望运行中的所有内容都在 python 中。

Answer 1

在Python2中，整数除法（使用/）return是一个整数。例如，1/3 = 0。在 Python3 中，整数除法（使用 /）可能 return 一个浮点数。

您似乎正在使用 Python2。要获得浮点数除法（在 Python2 和 Python3 中），确保每个除法运算至少涉及一个浮点数：例如，将 1/3 更改为 1.0/3 或 1/3.0 或（可接受但可能不太可读，1/3.）：

import numpy as np
def hertzmindlin(K0, G0, PHIC, P, f=1.0):
    K0, G0 = map(float, (K0, G0))
    V0 = (3*K0-2*G0)/(6*K0+2*G0) # Calculated theoretical poissons ratio of bulk rock
    n = 20-(34*PHIC)+(14*(PHIC**2)) # Coordination number at critical porosity (Murphy 1982)
    K_HM = (P*(n**2*(1-PHIC)**2*G0**2) / (18*np.pi**2*(1-V0)**2))**(1/3.0)
    G_HM = ((2+3*f-V0*(1+3*f))/(5*(2-V0))) * ((P*(3*n**2*(1-PHIC)**2*G0**2)/(2*np.pi**2*(1-V0)**2)))**(1/3.0)
    return K_HM, G_HM

K, G, = 36, 45
PHIC = 0.4
P = 0.001

print(hertzmindlin(K, G, PHIC, P))

或者，在 Python2 的更高版本中（例如 Python2.7），您可以放置

from __future__ import division

在脚本的顶部（在所有其他导入语句之前）到 activate Python3-style floating-point division。

numpy 函数给出不正确的结果 - 手动检查 excel

numpy function giving incorrect results - checked by hand and excel

math

numpy

geo