python - 将贝塞尔函数与 2 个变量积分

Question

我在执行此代码时遇到问题，我怀疑这是由于在贝塞尔函数 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'' 

Answer 1

发生错误的原因是您向 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当你有时间的时候。