首页 文章

使用Sympy lambdify函数实例化bessel函数时出错

提问于
浏览
1

我目前正在开发一个python程序,将 sympy 计算的符号表达式转换为包含所有数值的 numpy 数组 . 我使用 sympy.lambdify 函数实例化符号表达式 .

一些符号表达式包含贝塞尔函数,我将 scipy.special.jv/jy 等作为参数传递给 lambdify 函数 . 这是单个代码(来自程序的appart):

m = Symbol('m')
expr= -1.06048319593874*(-(3.14159265358979*m*besselj(27.9937791601866*m, 4.46681007624482*sqrt(-I)) - 0.501286290793831*sqrt(-I)*besselj(27.9937791601866*m + 1, 4.46681007624482*sqrt(-I)))*bessely(27.9937791601866*m, 6.81776274795262*sqrt(-I))/(3.14159265358979*m*bessely(27.9937791601866*m, 4.46681007624482*sqrt(-I)) - 0.501286290793831*sqrt(-I)*bessely(27.9937791601866*m + 1, 4.46681007624482*sqrt(-I))) + besselj(27.9937791601866*m, 6.81776274795262*sqrt(-I)))*sin(0.942966379693359*m)*cos(1.5707963267949*m)/m
nm = 5
vm = arange(nm) +1 
bessel = {'besselj':jv,'besselk':kv,'besseli':iv,'bessely':yv}
libraries = [bessel, "numpy"]     
result = lambdify(m, expr, modules=libraries)(vm)

In[1] : result
Out[1]: 
array([ -7.51212638e-030 -3.22606326e-030j,
         4.81143544e-046 +1.04405860e-046j,
         1.97977798e-097 +3.02047228e-098j,
         3.84986092e-124 +4.73598141e-125j,
         1.12934434e-181 +1.21145535e-182j])

结果如我所料:一个5行1-d数组,每个整数值为符号 m .

问题是当我尝试在程序中实现它时 . 这是实施:

expr = list_symbolic_expr[index]
vcm = arange(Dim_col[0]) +1   #Dim_col = list of integer range values to instantiate                     
m = Symbol(str(Var_col[0])) #Var_col = list of symbolic integer parameters
print expr, m
smat = lambdify(m, expr, modules=libraries)(vcm)

使用 scipy.special bessel函数时出现以下错误:

In [2]: run Get_System_Num
Out [2]:
-1.06048319593874*(-(3.14159265358979*m*besselj(27.9937791601866*m, 4.46681007624482*sqrt(-I)) - 0.501286290793831*sqrt(-I)*besselj(27.9937791601866*m + 1, 4.46681007624482*sqrt(-I)))*bessely(27.9937791601866*m, 6.81776274795262*sqrt(-I))/(3.14159265358979*m*bessely(27.9937791601866*m, 4.46681007624482*sqrt(-I)) - 0.501286290793831*sqrt(-I)*bessely(27.9937791601866*m + 1, 4.46681007624482*sqrt(-I))) + besselj(27.9937791601866*m, 6.81776274795262*sqrt(-I)))*sin(0.942966379693359*m)*cos(1.5707963267949*m)/m m
  File "Numeric\Get_System_Num.py", line 183, in get_Mat
    smat = lambdify(m, expr, modules=libraries)(vcm)

  File "<string>", line 1, in <lambda>

TypeError: ufunc 'jv' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

使用 sympy.special bessel函数时出现以下错误:

File "Numeric\Get_System_Num.py", line 183, in get_Mat
    smat = lambdify(m, expr, modules=libraries)(vcm)

  File "<string>", line 1, in <lambda>

  File "C:\Users\Emile\Anaconda2\lib\site-packages\sympy\core\function.py", line 375, in __new__
    result = super(Function, cls).__new__(cls, *args, **options)

  File "C:\Users\Emile\Anaconda2\lib\site-packages\sympy\core\function.py", line 199, in __new__
    evaluated = cls.eval(*args)

  File "C:\Users\Emile\Anaconda2\lib\site-packages\sympy\functions\special\bessel.py", line 161, in eval
    if nu.is_integer:

AttributeError: 'numpy.ndarray' object has no attribute 'is_integer'

似乎 lambdify 无法正确解释程序代码中bessel函数中的输入值(因为它是scipy或sympy bessel函数),因为输入值似乎是一个数组对象,而它不是单个问题码 . 但我没有看到它们之间的区别 .

感谢您阅读本文,我希望有人可以帮助我解决这个问题 . 请告诉我是否需要更多信息来说明这个问题 .

  • 编辑1--

我一直在搜索问题,甚至当我尝试lambdify这个简单的表达式: 17469169.5935065*sin(0.942966379693359*m)*cos(1.5707963267949*m)/m ,以"numpy"作为模块参数时,它给出: 'float' object has no attribute 'sin'

我也尝试了打包的模块'numexpr'没有成功,因为它在另一个表达式上停止: 0.058**(46.6321243523316*k) ,说:

File "C:\Users\Emile\Anaconda2\lib\site-packages\numexpr\necompiler.py", line 756, in evaluate
    signature = [(name, getType(arg)) for (name, arg) in zip(names, arguments)]

  File "C:\Users\Emile\Anaconda2\lib\site-packages\numexpr\necompiler.py", line 654, in getType
    raise ValueError("unknown type %s" % a.dtype.name)

ValueError: unknown type object

所以我真的不知道这个问题的本质,我甚至无法跟踪它,因为提出的异常来自内心的同情,我不知道如何在sympy函数中放置标志 .

  • 编辑2--

这是Get_System_Num类的主体:首先是导入,然后我也从符号项目中导入数据,这意味着:我要实例化的每个符号系数(在Mat_Num中列出);每个系数实例化(在“list_Var”中列出)的符号,取决于行和列;真实大矩阵中实例化矩阵的大小和位置(在list_Dim中列出) . 因为这些列表在开头是符号的我需要用我在Pr_Num <=>数值项目中放入的实数值替换符号,并在sub_system_num方法中执行 .

from numpy import diag, zeros, concatenate, arange, meshgrid, array
from scipy.special import jv,kv,iv,yv
from sympy import srepr, Matrix, lambdify, Symbol
from Numeric.Sub_System_Num import Sub_System_Num

class Get_System_Num :

    #Calling constructor
    def __init__(self, Pr_Num) :

        #Assigning inputs values to local values
        self.Pr_Num = Pr_Num
        self.Pr_Symb = Pr_Num.Pr_Symb

        #Load data from Pr_Symb              
        self.Mat_Num_Index = self.Pr_Symb.Mat_Num_Index        

        #Subs symbols with numeric values
        print "Substitute symbols with real values..."
        self.Sub_System_Num = Sub_System_Num(self)

        #Gathering the results        
        self.Mat_Num = self.Sub_System_Num.Mat_Num
        self.Mat_Size = self.Sub_System_Num.Mat_Size
        self.list_Index = self.Sub_System_Num.list_Index        
        self.list_Var_row = self.Sub_System_Num.list_Var.col(0)  
        self.list_Dim_row = self.Sub_System_Num.list_Dim.col(0)
        self.list_Var_col = self.Sub_System_Num.list_Var.col(1)  
        self.list_Dim_col = self.Sub_System_Num.list_Dim.col(1)

        self.count = 0

        print "Compute numerical matrix..."
        self.Mat = self.get_Mat()

这是填充数字矩阵的方法 . 我只假设只有一个变量来实例化行和列 . 它可能是9(3²)中的一种情况,因为每行和每列都有0,1或2个变量来实例化 . 关于 if srepr(coeff).__contains__(srepr(Var_row[0])) :这检查变量是否真的包含在expr中,如果不是,我手动复制coeff . 由于表达式是先前以符号方式计算的,因此有时候可以简化它们,以便变量不再在这里 .

def get_Mat(self) : 

    bessel = {'besselj':jv,'besselk':kv,'besseli':iv,'bessely':yv}
    libraries = [bessel, 'numpy']

    Mat = zeros((self.Mat_Size[0], self.Mat_Size[0]), dtype=complex)

    for index in range(0, self.Mat_Num_Index.__len__()) :

        Nb_row = self.list_Index[self.Mat_Num_Index[index][0], 0]
        Nb_col = self.list_Index[self.Mat_Num_Index[index][1], 1]

        Pos_row = self.list_Index[self.Mat_Num_Index[index][0], 2]
        Pos_col = self.list_Index[self.Mat_Num_Index[index][1], 3]

        Var_row = self.list_Var_row[self.Mat_Num_Index[index][0]]
        Var_col = self.list_Var_col[self.Mat_Num_Index[index][1]]

        Dim_row = self.list_Dim_row[self.Mat_Num_Index[index][0]]
        Dim_col = self.list_Dim_col[self.Mat_Num_Index[index][1]]

        coeff = self.Mat_Num[index]           

        #M(K or Z, M or Z)
        if Var_row.__len__() == 1 :    
            if Var_col.__len__() == 1 :                    
                if Var_row[0] == Var_col[0] :  #M(K,M=K) or M(Z,Z)  
                    vrk = arange(Dim_row[0]) +1  
                    k = Symbol(str(Var_row[0]))
                    smat = lambdify(k, coeff, modules=libraries)(vrk)                                           
                    if srepr(coeff).__contains__(srepr(Var_row[0])) :
                        smat = diag(smat)
                        print 'coeff n°', index, ' Case 1/1 M(K,M=K) or M(Z,Z) : i dependent '
                        print coeff, Var_col, Var_row

                    else :                           
                        smat = diag([smat]*Dim_row[0])                            
                        print 'coeff n°', index, ' Case 1/1 M(K,M=K) or M(Z,Z) : i non dependent'
                        print coeff, Var_col, Var_row                                  

                else :
                    if srepr(coeff).__contains__(srepr(Var_col[0]) and srepr(Var_row[0])) :  #M(K,Z) or M(Z,M) or M(K, M) 
                        print 'coeff n°', index, ' Case 1/1 M(K,Z) or M(Z,M) : i dependent ' 
                        print coeff, Var_col, Var_row, index                                                    
                        vrk = arange(Dim_row[0]) +1
                        vcm = arange(Dim_col[0]) +1
                        mcm, mrk = meshgrid(vcm, vrk)
                        k = Symbol(str(Var_row[0]))
                        i = Symbol(str(Var_col[0]))
                        smat = lambdify((k, i), coeff, modules=libraries)(mrk, mcm)

                    elif not(srepr(coeff).__contains__(srepr(Var_row[0]))) : #M(Z,M)
                        print 'coeff n°', index, ' Case 1/1 M(Z,M) : i non dependent '
                        print coeff, Var_col, Var_row, index                        
                        vcm = arange(Dim_col[0]) +1
                        m = Symbol(str(Var_col[0]))
                        smat = lambdify(m, coeff, modules=libraries)(vcm)
                        smat = [smat]*Dim_row[0]
                        smat = concatenate(list(smat), axis=0)

                    elif not(srepr(coeff).__contains__(srepr(Var_col[0]))) : #M(K,Z)
                        print 'coeff n°', index, ' Case 1/1 M(K,Z) : i non dependent'
                        print coeff, Var_col, Var_row, index                        
                        vrk = arange(Dim_row[0]) +1
                        k = Symbol(str(Var_row[0]))
                        smat = lambdify(k, coeff, modules=libraries)(vrk)
                        smat = [smat]*Dim_col[0]
                        smat = concatenate(list(smat), axis=1)                            

        self.count = self.count +1
        Mat[Pos_row:Pos_row+Nb_row, Pos_col:Pos_col+Nb_col] = smat

    return Mat

最后,我用lambdified coeff填充矩阵,在这种情况下,lambdified coeff是一个较小的矩阵(Nb_row,Nb_col) .

我将继续自己搜索,不要犹豫,要求更多细节!

  • 编辑3 -

我发现如果我这样做:

m = Symbol(str(Var_col[0]))
coeff = lambdify(m, coeff, modules=libraries)
for nm in range(0, Dim_col[0]) :
     smat.append(coeff(nm+1))

它可以工作,但它耗时更多(虽然远远少于使用subs和evalf) . 当我使用数组对象调用lambdify创建的lambda函数时,会出现错误(无论是numpy还是sympy数组类型) .

1 回答

  • 0

    我可以重现你的错误 . 基于这些输入:

    m = Symbol('m')
    expr = -1.06048319593874*(-(3.14159265358979*m*besselj(27.9937791601866*m, 4.46681007624482*sqrt(-I)) - 0.501286290793831*sqrt(-I)*besselj(27.9937791601866*m + 1, 4.46681007624482*sqrt(-I)))*bessely(27.9937791601866*m, 6.81776274795262*sqrt(-I))/(3.14159265358979*m*bessely(27.9937791601866*m, 4.46681007624482*sqrt(-I)) - 0.501286290793831*sqrt(-I)*bessely(27.9937791601866*m + 1, 4.46681007624482*sqrt(-I))) + besselj(27.9937791601866*m, 6.81776274795262*sqrt(-I)))*sin(0.942966379693359*m)*cos(1.5707963267949*m)/m
    nm = 2
    bessel = {'besselj': jv,'besselk':kv,'besseli':iv,'bessely':yv}
    libraries = [bessel, "numpy"]
    

    它在 vm 中运行整数:

    vm = np.arange(nm, dtype=np.int) +1 
    result = lambdify(m, expr, modules=libraries)(vm)
    

    但复数会产生与您相同的错误消息:

    vm = np.arange(nm, dtype=np.complex) +1 
    result = lambdify(m, expr, modules=libraries)(vm)
    

    错误消息:

    TypeError                                 Traceback (most recent call last)
    <ipython-input-30-f0e31009275a> in <module>()
          1 vm = np.arange(nm, dtype=np.complex) +1
    ----> 2 result = lambdify(m, expr, modules=libraries)(vm)
    
    /Users/mike/anaconda/envs/py34/lib/python3.4/site-packages/numpy/__init__.py in <lambda>(_Dummy_27)
    
    TypeError: ufunc 'jv' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
    

    因此,检查 vcm 中的数据类型,即 arange(Dim_col[0]) +1 返回的内容 .

    您只能使用属性 real 获取复数的实部 . 例如, vm.real 为您提供 vm 的真实部分 . 如果这是你想要的,你应该得到你的结果:

    result = lambdify(m, expr, modules=libraries).(vm.real)
    

相关问题