【问题标题】：Defining a integral of multi variable function as a second function python将多变量函数的积分定义为第二个函数python
【发布时间】：2020-04-14 03:06:37
【问题描述】：

将多变量函数的积分定义为第二个函数python

我正在使用 python 仅对其中一个变量积分一个多变量函数（它是 x 和 theta 的函数，我正在对从 0 到 2*pi 的 theta 积分，因此结果是 x 的函数）。我尝试了以下方法：

import numpy as np
import scipy.integrate as inte
d=10.0
xvals=np.linspace(-d,d,1000)
def aIntegrand(theta,x):
    return 1/(2*np.pi)*np.sin(2*np.pi*x*np.sin(theta)/d)**2
def A(x):
    return (inte.quad(aIntegrand,0,2*np.pi,args=(x,))[0])**(1/2)

plt.plot(xvals,A(xvals))
plt.xlabel("x")
plt.ylabel("A(x)")
plt.show()

我收到以下错误：

TypeError: only size-1 arrays can be converted to Python scalars

我认为这是因为四元积分器的结果是一个包含两个元素的数组，而 python 不喜欢基于索引数组定义函数？不过，这是对问题的完整猜测。如果有人知道我该如何解决这个问题并可以让我知道，那就太好了:)

第二次尝试

我已经成功地使用以下代码获得了积分图：

import numpy as np
import scipy.integrate as inte
import matplotlib.pyplot as plt
d=10.0
xvals=np.linspace(-d,d,1000)
thetavals=np.linspace(0.0,2*np.pi,1000)
def aIntegrand(theta,x):
    return 1/(2*np.pi)*np.sin(2*np.pi*x*np.sin(theta)/d)**2
def A(x):
    result=np.zeros(len(x))
    for i in range(len(x)):
        result[i]=(inte.quad(aIntegrand,0,2*np.pi,args=(x[i],))[0])**(1/2)
    return result
def f(x,theta):
    return x**2* np.sin(theta)

plt.plot(xvals,A(xvals))
plt.xlabel("x")
plt.ylabel("A(x)")
plt.show()

但这并没有将 A(x) 作为函数给出，由于我定义它的方式，它需要一个数组形式的输入。我需要该函数与 aIntegrand 具有相同的形式，其中当给定参数返回单个值时，该函数可以重复集成。

【问题讨论】：

标签： python function scipy integration quad

【解决方案1】：

我不认为你所寻求的存在于 Scipy 中。但是，您至少有两种选择。

首先，您可以使用interp1d 创建一个插值器。这意味着您最初给出的 x 值的范围决定了您将能够调用插值器的范围。然后可以针对此区间内的任何 x 值调用此插值器。
即使您没有可调用对象，也可以执行第二个积分。函数simps 只需要值及其位置来提供对集成过程的估计。否则，您可以通过dblquad 一次通话执行双重积分。

【讨论】：

【解决方案2】：

首先请注意，您的积分可以通过解析计算。这是

0.5 * (1 - J0(4 * pi * x / d))

其中J0 是第一类贝塞尔函数。

其次，你可以使用quadpy（我的一个项目）；它具有完全矢量化的计算。

import numpy as np
import quadpy
import matplotlib.pyplot as plt
import scipy.special

d = 10.0
x = np.linspace(-d, d, 1000)


def aIntegrand(theta):
    return (
        1
        / (2 * np.pi)
        * np.sin(2 * np.pi * np.multiply.outer(x, np.sin(theta)) / d) ** 2
    )


Ax2, err = quadpy.quad(aIntegrand, 0, 2 * np.pi)
Ax = np.sqrt(Ax2)

plt.plot(x, Ax, label="quadpy")
plt.xlabel("x")
plt.ylabel("A(x)")
plt.plot(x, np.sqrt(0.5 * (1 - scipy.special.jv(0, 4*np.pi*x/d))), label="bessel")
# plt.show()
plt.savefig("out.png", transparent=True, bbox_inches="tight")

【讨论】：