我在执行此代码时遇到问题,我怀疑这是由于在贝塞尔函数 J1 内使用带有多个变量的 lambda 函数而导致的问题。函数 q(theta) 给出错误。 INTEGRATEZI只是一种从0到10**10的积分方法(Gauss-Kronrod方法)。
from scipy.special import j1 as J1
from cmath import sin, cos
import scipy
from scipy import array
def cml(function):
return (scipy.real(function)**2 - scipy.imag(function)**2)
def quad_routine(func, a, b, x_list, w_list):
c_1 = (b-a)/2.0
c_2 = (b+a)/2.0
eval_points = map(lambda x: c_1*x+c_2, x_list)
func_evals = map(func, eval_points)
return c_1 * sum(array(func_evals) * array(w_list))
def quad_kronrod_15(func, a, b):
x_kr = [-0.991455371120813,-0.949107912342759, -0.864864423359769, -0.741531185599394, -0.586087235467691,-0.405845151377397, -0.207784955007898, 0.0, 0.207784955007898,0.405845151377397, 0.586087235467691, 0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813]
w_kr = [0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525, 0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728, 0.204432940075298, 0.190350578064785, 0.169004726639267, 0.140653259715525, 0.104790010322250, 0.063092092629979, 0.022935322010529]
return quad_routine(func,a,b,x_kr, w_kr)
class Memoize(object):
def __init__(self, func):
self.func = func
self.eval_points = {}
def __call__(self, *args):
if args not in self.eval_points:
self.eval_points[args] = self.func(*args)
return self.eval_points[args]
def quadt(func,a,b):
func = Memoize(func) # Memoize function to skip repeated function calls.
k15 = quad_kronrod_15(func,a,b)
return k15
def INTEGRATEZI(func, **kwargs): # INTEGRATEZI short for Integrating from zero to infinity
def real_func(x):
return scipy.real(func(x))
def imag_func(x):
return scipy.imag(func(x))
real_integral = quadt(real_func, 0, 10**10)
imag_integral = quadt(imag_func, 0, 10**10)
return (real_integral + 1j*imag_integral)
r = 2.0
p = lambda theta: r*sin(theta)
z = lambda theta: r*cos(theta)
q = lambda theta: (p(theta) + z(theta) + INTEGRATEZI(lambda kp: (J1(kp * p(theta))) ))
print q(0.5)
这是错误报告:
%run "C:/Users/ENVY14-i7-SPECTRE/Documents/University/Year 3/project/123123.py"
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
C:\Users\ENVY14-i7-SPECTRE\AppData\Local\Enthought\Canopy\App\appdata\canopy-1.4.1.1975.win-x86_64\lib\site-packages\IPython\utils\py3compat.pyc in execfile(fname, glob, loc)
195 else:
196 filename = fname
--> 197 exec compile(scripttext, filename, 'exec') in glob, loc
198 else:
199 def execfile(fname, *where):
C:\Users\ENVY14-i7-SPECTRE\Documents\University\Year 3\project\123123.py in <module>()
45 z = lambda theta: r*cos(theta)
46 q = lambda theta: (p(theta) + z(theta) + INTEGRATEZI(lambda kp: (J1(kp * p(theta))) ))
---> 47 print q(0.5)
C:\Users\ENVY14-i7-SPECTRE\Documents\University\Year 3\project\123123.py in <lambda>(theta)
44 p = lambda theta: r*sin(theta)
45 z = lambda theta: r*cos(theta)
---> 46 q = lambda theta: (p(theta) + z(theta) + INTEGRATEZI(lambda kp: (J1(kp * p(theta))) ))
47 print q(0.5)
C:\Users\ENVY14-i7-SPECTRE\Documents\University\Year 3\project\123123.py in INTEGRATEZI(func, **kwargs)
38 def imag_func(x):
39 return scipy.imag(func(x))
---> 40 real_integral = quadt(real_func, 0, 10**10)
41 imag_integral = quadt(imag_func, 0, 10**10)
42 return (real_integral + 1j*imag_integral)
C:\Users\ENVY14-i7-SPECTRE\Documents\University\Year 3\project\123123.py in quadt(func, a, b)
30 def quadt(func,a,b):
31 func = Memoize(func) # Memoize function to skip repeated function calls.
---> 32 k15 = quad_kronrod_15(func,a,b)
33 return k15
34
C:\Users\ENVY14-i7-SPECTRE\Documents\University\Year 3\project\123123.py in quad_kronrod_15(func, a, b)
17 x_kr = [-0.991455371120813,-0.949107912342759, -0.864864423359769, -0.741531185599394, -0.586087235467691,-0.405845151377397, -0.207784955007898, 0.0, 0.207784955007898,0.405845151377397, 0.586087235467691, 0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813]
18 w_kr = [0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525, 0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728, 0.204432940075298, 0.190350578064785, 0.169004726639267, 0.140653259715525, 0.104790010322250, 0.063092092629979, 0.022935322010529]
---> 19 return quad_routine(func,a,b,x_kr, w_kr)
20
21 class Memoize(object):
C:\Users\ENVY14-i7-SPECTRE\Documents\University\Year 3\project\123123.py in quad_routine(func, a, b, x_list, w_list)
11 c_2 = (b+a)/2.0
12 eval_points = map(lambda x: c_1*x+c_2, x_list)
---> 13 func_evals = map(func, eval_points)
14 return c_1 * sum(array(func_evals) * array(w_list))
15
C:\Users\ENVY14-i7-SPECTRE\Documents\University\Year 3\project\123123.py in __call__(self, *args)
25 def __call__(self, *args):
26 if args not in self.eval_points:
---> 27 self.eval_points[args] = self.func(*args)
28 return self.eval_points[args]
29
C:\Users\ENVY14-i7-SPECTRE\Documents\University\Year 3\project\123123.py in real_func(x)
35 def INTEGRATEZI(func, **kwargs): # INTEGRATEZI short for Integrating from zero to infinity
36 def real_func(x):
---> 37 return scipy.real(func(x))
38 def imag_func(x):
39 return scipy.imag(func(x))
C:\Users\ENVY14-i7-SPECTRE\Documents\University\Year 3\project\123123.py in <lambda>(kp)
44 p = lambda theta: r*sin(theta)
45 z = lambda theta: r*cos(theta)
---> 46 q = lambda theta: (p(theta) + z(theta) + INTEGRATEZI(lambda kp: (J1(kp * p(theta))) ))
47 print q(0.5)
TypeError: ufunc 'j1' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
最佳答案
发生错误的原因是您向 Bessel 函数提供的参数:
kp * p(theta)
不是实数。然后 Scipy 会抛出一个类型错误,因为它需要一个实数。
为什么不是真的?您正在使用 Python 的复数数学库 cmath
,其 sin
函数以复数形式返回结果。因此,例如,
In [9]: cmath.sin(.5)
Out[9]: (0.479425538604203+0j)
有一个(零)虚部。但是 Scipy 的 Bessel 函数接收到一个复杂类型的变量,并且不知道如何处理它。快速解决方法是使用标准 math
库,或将 r*sin(theta)
替换为 r*sin(theta).real
并p
和 q
定义中的 r*cos(theta)
通过 r*cos(theta).real
.
要查看此内容,您可以使用类似 python debugger 的工具找出引发异常的确切原因。看起来您正在使用 IPython,其中包含一些易于使用的调试工具。
在您的情况下,在 IPython 终端运行代码,您将得到一个异常:
In [4]: run code.py
... (cut for space) ...
45 p = lambda theta: r*sin(theta)
46 z = lambda theta: r*cos(theta)
---> 47 q = lambda theta: (p(theta) + z(theta) + INTEGRATEZI(lambda kp: (J1(kp * p(theta))) ))
48 print q(0.5)
TypeError: ufunc 'j1' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
您将看到在使用 --->
评估该行时引发了 TypeError
。但这一行发生了很多事情。究竟是什么导致了这个问题?
一种简单的检查方法是“事后分析”启动调试器。每个变量的值都会像异常发生时一样及时“卡住”,因此我们可以四处查看问题可能出在哪里。
运行代码并看到异常后,立即在 IPython 提示符处输入 %debug
:
In [5]: %debug
> /home/eldridge/sandbox/stack/stack.py(47)<lambda>()
46 z = lambda theta: r*cos(theta)
---> 47 q = lambda theta: (p(theta) + z(theta) + INTEGRATEZI(lambda kp: (J1(kp * p(theta))) ))
48 print q(0.5)
我们处于导致错误的框架中。您认为问题可能出在贝塞尔函数上,所以我决定询问争论是什么。您可以使用 p
命令打印变量:
ipdb> p kp
42723144.39593506
看起来不错。 p(theta)
怎么样?
ipdb> p p(theta)
(0.958851077208406+0j)
嗯。这是想象的。贝塞尔函数可以处理复杂的输入吗?我不记得我的数学类(class),但我可以尝试一下看看:
ipdb> p J1(1j)
*** TypeError: TypeError("ufunc 'j1' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''",)
显然不是。
当然,调试器还有更多功能,值得查看 tutorial当你有时间的时候。
关于python - 将贝塞尔函数与 2 个变量积分,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/27360726/