numpy 球面上的快速二重积分

wqsoz72f 于 12个月前发布在其他

关注(0)|答案(1)|浏览(113)

我需要在一个单位球面上进行5次数值积分。
我试过dblquad，但这需要相当长的时间。
下面是一个MRE：

import numpy as np
from scipy import integrate

def foo():

    f1 = lambda theta, phi: np.cos(theta)**2
    f2 = lambda theta, phi: f1(theta, phi)*np.cos(theta)**2*np.sin(theta)  
    
    ret = np.zeros(5)
    
    ret[0] = integrate.dblquad(f2, 0, 2*np.pi, 0, np.pi)[0]
    ret[1] = integrate.dblquad(f2, 0, 2*np.pi, 0, np.pi)[0]
    ret[2] = integrate.dblquad(f2, 0, 2*np.pi, 0, np.pi)[0]
    ret[3] = integrate.dblquad(f2, 0, 2*np.pi, 0, np.pi)[0]
    ret[4] = integrate.dblquad(f2, 0, 2*np.pi, 0, np.pi)[0]

    return ret
    
def main():
    tst = np.zeros((1000, 5))
    for i in range(1000):
        tst[i] = foo()
        
    
main()

在我的机器中，当我执行%timeit main()时，我得到21.4 s ± 126 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
这个代码结构需要很长的时间来执行计算，有什么方法可以改进吗？

numpy

来源：https://stackoverflow.com/questions/76965342/fast-double-integration-over-a-sphere

1条答案

按热度按时间

ie3xauqp1#

你可以通过使用numba和jit f1/f2函数来加快计算速度：

import numpy as np
from numba import njit
from scipy import integrate

@njit
def f1(theta, phi):
    return np.cos(theta) ** 2

@njit
def f2(theta, phi):
    return f1(theta, phi) * np.cos(theta) ** 2 * np.sin(theta)

def foo():
    ret = np.zeros(5)

    ret[0] = integrate.dblquad(f2, 0, 2 * np.pi, 0, np.pi)[0]
    ret[1] = integrate.dblquad(f2, 0, 2 * np.pi, 0, np.pi)[0]
    ret[2] = integrate.dblquad(f2, 0, 2 * np.pi, 0, np.pi)[0]
    ret[3] = integrate.dblquad(f2, 0, 2 * np.pi, 0, np.pi)[0]
    ret[4] = integrate.dblquad(f2, 0, 2 * np.pi, 0, np.pi)[0]

    return ret

def main():
    tst = np.zeros((1000, 5))
    for i in range(1000):
        tst[i] = foo()

main()

在我的机器（AMD 5700 x/Python 3.11）上使用time运行它：

andrej@MyPC:/app$ time python3 script.py

real    0m2,204s
user    0m0,013s
sys     0m0,004s

作为比较，如果没有@njit装饰器：

andrej@MyPC:/app$ time python3 script.py

real    0m13,505s
user    0m0,014s
sys     0m0,004s

赞(0）回复(0）举报 12个月前

我来回答

numpy 球面上的快速二重积分

1条答案

相关问题

热门标签

最新问答