まず、SymPy と NumPy の違いをご理解ください。SymPy はシンボリック計算のためのライブラリで、NumPy は数値計算のためのライブラリです。
その上で、SymPy にはベクトル場を表すための関数が用意されています。このドキュメントをご覧ください:Scalar and Vector Field Functionality
サンプルプログラムです:
>>> from sympy.vector import CoordSys3D, Del, divergence
>>> # 座標系を決めます
... C = CoordSys3D('C')
>>> # ナブラは Del() として定義されています
... nabla = Del()
>>> # C.i, C.j, C.k が標準的な基底、C.x, C.y, C.z がベクトル場の仮引数です
... b = C.x * C.i + C.y * C.j + C.z * C.k
>>> # 内積は dot です
... nabla.dot(b)
Derivative(C.x, C.x) + Derivative(C.y, C.y) + Derivative(C.z, C.z)
>>> # 微分計算をしてもらうためには doit() を使います
... nabla.dot(b).doit()
3
>>> # ナブラとの内積は divergence としても定義されています
... divergence(b).doit()
3