将多变量函数的积分定义为第二函数python
Defining a integral of multi variable function as a second function python
将多变量函数的积分定义为第二函数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) 作为函数给出,由于我定义它的方式,它需要一个数组形式的输入。我需要函数与被积函数具有相同的形式,其中当给定参数时 returns 单个值 & 因此函数可以重复积分。
我认为Scipy中不存在您要寻找的东西。但是,您至少有两个选择。
首先请注意,您的积分可以通过分析计算。是
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")
将多变量函数的积分定义为第二函数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) 作为函数给出,由于我定义它的方式,它需要一个数组形式的输入。我需要函数与被积函数具有相同的形式,其中当给定参数时 returns 单个值 & 因此函数可以重复积分。
我认为Scipy中不存在您要寻找的东西。但是,您至少有两个选择。
首先请注意,您的积分可以通过分析计算。是
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")