1

我正在尝试使用quadpy，因为我想对 2D 积分进行 2D 数值矢量积分。为了了解工作速度有多快quadpy，我想对其进行测试并将其与scipy一维矢量积分进行比较。于是我写了一个简单的代码：

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import quadpy

q = np.linspace(0.03, 1.0, 500)


def f(t):
    return t * 0.5 * (erf((t - 40) / 3) - 1) * j0(q * t)


y, _ = integrate.quad_vec(f, 0, 50)
y1, _ = quadpy.quad(f, 0, 50)

print(y)
print(y1)

没有使用 quadpy 的经验，我收到以下错误：

Traceback (most recent call last):
  File "test.py", line 8, in <module>
    y1 = quadpy.quad(lambda t: t*0.5*(erf((t-40)/3)-1)*j0(q*t),0.0,50)
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/_scipy_compat.py", line 44, in quad
    g, [a, b], eps_abs=epsabs, eps_rel=epsrel, max_num_subintervals=limit
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/line_segment/_tools.py", line 42, in integrate_adaptive
    range_shape=range_shape,
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/line_segment/_gauss_kronrod.py", line 124, in _gauss_kronrod_integrate
    fx_gk = numpy.asarray(f(sp))
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/_scipy_compat.py", line 41, in g
    return f(x, *args)
  File "test.py", line 8, in <lambda>
    y1 = quadpy.quad(lambda t: t*0.5*(erf((t-40)/3)-1)*j0(q*t),0.0,50)
ValueError: operands could not be broadcast together with shapes (500,) (15,)
4

1 回答 1

2

你缺少一个multiply.outer：

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import quadpy

q = np.linspace(0.03, 1.0, 500)


def f(t):
    return t * 0.5 * (erf((t - 40) / 3) - 1) * j0(np.multiply.outer(q, t))


y, _ = integrate.quad_vec(f, 0, 50)
y1, _ = quadpy.quad(f, 0, 50)

print(y - y1)
于 2020-03-29T12:33:42.673 回答