使用Quadpy进行集成时出现问题 [英] Problems integrating using quadpy

问题描述

我想使用quadpy，因为我想对2D积分进行2D数值矢量积分.为了了解quadpy的运行速度，我想对其进行测试并将其与scipy 1D矢量积分进行比较.因此，我编写了一个简单的代码:

I am trying to use quadpy as I want to do 2D numerical vector integration for 2D integrals. To see how fast quadpy works, I wanted to test it and compare it with scipy 1D vector integration. Thus I wrote a simple code:

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的经验，我遇到以下错误:

Having no experience with quadpy, I get the following errors:

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,) 

推荐答案

您缺少一个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)

这篇关于使用Quadpy进行集成时出现问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！