如何在sympy中实现div和grad?
最近我在玩一个叫做sympy的工具,遇到了一个问题,就是如何计算标量场和向量场的散度和梯度。我现在想做的是求解热方程,也就是:
d/dt u(x,t) - a * lap(u(x,t)) = 0,其中lap(x) = (div(grad(x))
这是在一个标量场上进行的。因为我在sympy.physics.mechanics里找不到lap、div和grad这几个函数,所以我尝试自己实现它们。
from sympy import *
from sympy.physics.mechanics import *
o = ReferenceFrame('o')
x,y,z = symbols('x y z')
class div(Function):
@classmethod
def eval(cls, F):
return F.diff(x, o).dot(o.x)+F.diff(y, o).dot(o.y)+F.diff(z, o).dot(o.z)
class grad(Function):
@classmethod
def eval(cls, F):
return o.x * F.diff(x) + o.y * F.diff(y) + o.z * F.diff(z)
f = x**2 + y**3 + z*0.5
print grad(f)
print type(grad(f))
print div(grad(f))
不幸的是,这样做的结果是:
2*x*o.x + 3*y**2*o.y + 0.500000000000000*o.z
Traceback (most recent call last):
File "/home/fortmeier/Desktop/autokernel/autokernel/tools/GenerateCode.py", line 24, in <module>
print div(grad(f))
File "/usr/local/lib/python2.7/dist-packages/sympy/core/cache.py", line 93, in wrapper
r = func(*args, **kw_args)
File "/usr/local/lib/python2.7/dist-packages/sympy/core/function.py", line 368, in __new__
result = super(Function, cls).__new__(cls, *args, **options)
File "/usr/local/lib/python2.7/dist-packages/sympy/core/cache.py", line 93, in wrapper
r = func(*args, **kw_args)
File "/usr/local/lib/python2.7/dist-packages/sympy/core/function.py", line 188, in __new__
args = list(map(sympify, args))
File "/usr/local/lib/python2.7/dist-packages/sympy/core/sympify.py", line 313, in sympify
expr = parse_expr(a, local_dict=locals, transformations=transformations, evaluate=evaluate)
File "/usr/local/lib/python2.7/dist-packages/sympy/parsing/sympy_parser.py", line 757, in parse_expr
return eval_expr(code, local_dict, global_dict)
File "/usr/local/lib/python2.7/dist-packages/sympy/parsing/sympy_parser.py", line 691, in eval_expr
code, global_dict, local_dict) # take local objects in preference
File "<string>", line 1, in <module>
AttributeError: 'Symbol' object has no attribute 'x'
[Finished in 0.3s with exit code 1]
我知道我可以使用galgebra模块来做一些事情,但我首先想更清楚地理解这里发生了什么。所以我的问题是,我缺少了什么呢?
1 个回答
2
这看起来像是SymPy里的一个错误。sympify(ReferenceFrame('o').x)
这个代码不管用。