如何在sympy中实现div和grad?

1 投票
1 回答
1026 浏览
提问于 2025-04-18 01:49

最近我在玩一个叫做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) 这个代码不管用。

撰写回答