我目前正在开发python程序,以将由sympy计算的符号表达式转换为包含所有数值的numpy数组。我用sympy.lambdify函数实例化符号表达式。
一些符号表达式包含Bessel函数,我将scipy.special.jv/jy等作为参数传递给lambdify函数。这是单个代码(程序中的appart):
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15m = 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维数组,每个整数m值。
问题是当我尝试在程序中实现它时。这是实现:
1
2
3
4
5expr = 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贝塞尔函数时的错误:
1
2
3
4
5
6
7
8
9In [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"", line 1, in
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贝塞尔函数时的错误:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15File"Numeric\Get_System_Num.py", line 183, in get_Mat
smat = lambdify(m, expr, modules=libraries)(vcm)
File"", line 1, in
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--
我一直在寻找问题,甚至当我尝试对这个简单的表达式:17469169.5935065*sin(0.942966379693359*m)*cos(1.5707963267949*m)/m进行lambdify时,以" numpy"作为模块参数,它也会给出:'float' object has no attribute 'sin'
我还尝试了打包模块'numexpr',但没有成功,因为它停止在另一个表达式上:0.058**(46.6321243523316*k),说:
1
2
3
4
5
6
7
8
9
10
11File"C:\Users\Emile\Anaconda2\lib\site-packages
umexpr
ecompiler.py", line 756, in evaluate
signature = [(name, getType(arg)) for (name, arg) in zip(names, arguments)]
File"C:\Users\Emile\Anaconda2\lib\site-packages
umexpr
ecompiler.py", line 654, in getType
raise ValueError("unknown type %s" % a.dtype.name)
ValueError: unknown type object
因此,我真的不知道此问题的性质,甚至无法跟踪它,因为引发的异常来自sympy内部,而且我也不知道如何在sympy函数中放置标志。
-编辑2--
这是Get_System_Num类的主体:首先是导入,然后我还从符号项目中导入数据,这意味着:我要实例化的每个符号系数(列在Mat_Num中);每个系数要实例化的符号(列在" list_Var"中),具体取决于行和列;实例化矩阵在实大矩阵中的大小和位置(列在list_Dim中),因为这些列表在开头是符号性的我需要用放在Pr_Num <=>数值项目中的实数值替换符号,然后在sub_system_num方法中执行。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34from 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(32)种情况中的一种,因为每一行和每一列都有0,1或2个要实例化的变量。关于if srepr(coeff).__contains__(srepr(Var_row[0])):这检查变量是否确实包含在expr中,如果不是,则手动复制coeff。由于表达式以前是通过符号计算的,因此有时sympy可能会简化它们,从而使该变量不再存在。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73def 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
最后,我用拉姆化的系数填充矩阵,在这种情况下,矩阵的大小较小(Nb_row,Nb_col)。
我将继续自己寻找,不要犹豫,要求更多详细信息!
-编辑3-
我发现如果执行此操作:
1
2
3
4m = 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要少得多)。当我使用数组对象(无论是numpy还是sympy数组类型)调用lambdify创建的lambda函数时,就会出现错误。
您可以显示整个文件Get_System_Num.py吗? 根本无法猜测程序中正在发生的事情。
只要看到它不是不可能的! 我已经编辑了我的帖子,以便您可以看到课程的其余部分。 非常感谢您的宝贵时间。
我可以重现您的错误。 基于以下输入:
1
2
3
4
5m = 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中使用整数:
1
2vm = np.arange(nm, dtype=np.int) +1
result = lambdify(m, expr, modules=libraries)(vm)
但是复数会产生与您相同的错误消息:
1
2vm = np.arange(nm, dtype=np.complex) +1
result = lambdify(m, expr, modules=libraries)(vm)
错误信息:
1
2
3
4
5
6
7
8TypeError Traceback (most recent call last)
in ()
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 (_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的实部。 如果这是您想要的,则应该获得以下结果:
1result = lambdify(m, expr, modules=libraries).(vm.real)
而且有效! 非常感谢你 ! 请问您如何更改arange数据类型?
添加了更新。