scipy Python中矩阵值的数值积分

xqkwcwgp  于 2023-06-23  发布在  Python
关注(0)|答案(4)|浏览(155)

我试图在Python中集成一些矩阵条目。我希望避免循环,因为我的任务包括100万次模拟。我正在寻找一个规范,将有效地解决我的问题。
我得到以下错误:只能将size-1数组转换为Python标量

from scipy import integrate
import numpy.random as npr

n = 1000
m = 30

x = npr.standard_normal([n, m])

def integrand(k):
    return k * x ** 2

integrate.quad(integrand, 0, 100)

这是一个简单的例子,我的情况。我有多个嵌套函数,所以我不能简单地把x放在积分前面。

guz6ccqo

guz6ccqo1#

你可能想使用并行执行来实现这一点。只要你只想执行integrate.quad30000000次,这应该很容易。只需将您的工作负载拆分为小包并将其交给线程池即可。当然,加速是有限的核心数量,你在你的个人电脑。我不是Python程序员,但这是可能的。你也可以在quad函数中增加epsabs和epsrel参数,这取决于实现,这也会加快程序的速度。当然,你会得到一个不太精确的结果,但这可能是确定的,这取决于你的问题。

import threading
from scipy import integrate
import numpy.random as npr

n = 2
m = 3
x = npr.standard_normal([n,m])

def f(a):
    for j in range(m):
        integrand = lambda k: k * x[a,j]**2
        i =integrate.quad(integrand, 0, 100)
        print(i) ##write it to result array

for i in range(n):
    threading.Thread(target=f(i)).start();
##better split it up even more and give it to a threadpool to avoid
##overhead because of thread init
oalqel3c

oalqel3c2#

这可能不是理想的解决方案,但它应该有所帮助。可以使用numpy.vectorize。连医生都说 * 提供向量化功能主要是为了方便,而不是为了性能。实现本质上是一个for循环。* 但是,在您提供的简单示例中,%timeit显示了2.3倍的加速。
实现是

from scipy import integrate
from numpy import vectorize
import numpy.random as npr  

n = 1000
m = 30

x = npr.standard_normal([n,m])

def g(x):
    integrand = lambda k: k * x**2
    return integrate.quad(integrand, 0, 100)

vg = vectorize(g)
res = vg(x)
shstlldc

shstlldc3#

2003更新:四矢量

截至2023年,有一个Scipy函数integrate.quad_vec用于向量函数的有效正交。
这个问题的解决方案是以下高度矢量化的过程

from scipy import integrate
import numpy as np

x = np.random.standard_normal([1000, 30])

def integrand(k):
    return k * x**2

res = integrate.quad_vec(integrand, 0, 100)

输出res[0]包含一个1000x30矩阵,其中包含每个参数x的数值积分。

nx7onnlm

nx7onnlm4#

quadpy(一个矿山项目,商业项目)做矢量化正交:

import numpy
import numpy.random as npr
import quadpy

x = npr.standard_normal([1000, 30])

def integrand(k):
    return numpy.multiply.outer(x ** 2, k)

scheme = quadpy.line_segment.gauss_legendre(10)
val = scheme.integrate(integrand, [0, 100])

这比所有其他答案都要快得多。

相关问题