使用Python计算矢量场的旋度并使用matplotlib绘制。

Question

我需要计算矢量场的旋度并使用matplotlib绘制它。一个简单的例子如下：

如何在matplotlib图库中的quiver3d_demo.py中计算和绘制矢量场的旋度？

- celestos

---

似乎在numpy或scipy中没有内置的函数来计算旋度。如果没有，您将不得不编写它；在三维中，结果也将是一个矢量场，因此matplotlib会像示例中一样绘制它。类似的关于散度的问题。 - cphlewis

谢谢您的回复。好的，我知道如何绘制了。您有关于编写curl函数的好建议吗？ - celestos

也许使用Sympy（http://docs.sympy.org/dev/modules/physics/vector/fields.html）会很有用。它具有一些内置的矢量场功能。 - Dietrich

我看到Sympy有一个curl函数，并且还处理绘图（使用matplotlib作为后端），但我从未使用过它。 - cphlewis

Answer 1

你可以使用sympy.curl()来计算矢量场的旋度。

示例：

假设F(x,y,z) = y²zi - xyj + z²k，则：

• y将是R[1]，x是R[0]，z是R[2]
• 三个轴的单位矢量i，j和k分别为R.x，R.y和R.z。

计算矢量场旋度的代码如下：

from sympy.physics.vector import ReferenceFrame
from sympy.physics.vector import curl
R = ReferenceFrame('R')

F = R[1]**2 * R[2] * R.x - R[0]*R[1] * R.y + R[2]**2 * R.z

G = curl(F, R)  

那么，如果G等于R_y**2*R.y + (-2*R_y*R_z - R_y)*R.z，换句话说,
G = 0i + y²j + (-2yz-y)k.

要绘制它，您需要将上述结果转换为3个单独的函数：u、v和w。

(以下示例改编自此 matplotlib 示例):

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = fig.gca(projection='3d')

x, y, z = np.meshgrid(np.arange(-0.8, 1, 0.2),
                      np.arange(-0.8, 1, 0.2),
                      np.arange(-0.8, 1, 0.8))

u = 0
v = y**2
w = -2*y*z - y

ax.quiver(x, y, z, u, v, w, length=0.1)

plt.show()

最终结果如下：

- user

---

请解释一下为什么会有负评。我很乐意修正任何错误。 - user

Answer 2

为计算向量函数的旋度，您也可以使用numdifftools进行自动数值微分，而无需通过符号微分的绕路。Numdifftools不提供curl()函数，但它可以计算一个或多个变量的向量值函数的Jacobian矩阵，并提供相对于所有变量的向量场所有分量的导数；这就是计算旋度所必需的全部内容。

import import scipy as sp
import numdifftools as nd

def h(x):
    return sp.array([3*x[0]**2,4*x[1]*x[2]**3, 2*x[0]])

def curl(f,x):
    jac = nd.Jacobian(f)(x)
    return sp.array([jac[2,1]-jac[1,2],jac[0,2]-jac[2,0],jac[1,0]-jac[0,1]])

x = sp.array([1,2,3)]
curl(h,x)

这将返回 x 处的 curl 值: array([-216., -2., 0.])。按照上述建议进行绘图即可。

- Martin Ligare

---

我更喜欢这个解决方案而不是被接受的答案，因为我很少处理分析向量场。然而，它有点低效，因为卷积没有使用Jacobian的迹（因此在nd.Jacobian中的1/3工作在这里没有使用）。您认为使用nd.directionaldiff()独立计算Jacobian的非对角线元素会更快吗？还是您认为额外的Python代码开销使其不值得？ - NLi10Me

Answer 3

这里有一份基于Octave / Matlab实现的Python代码。（链接）

import numpy as np

def curl(x,y,z,u,v,w):
    dx = x[0,:,0]
    dy = y[:,0,0]
    dz = z[0,0,:]

    dummy, dFx_dy, dFx_dz = np.gradient (u, dx, dy, dz, axis=[1,0,2])
    dFy_dx, dummy, dFy_dz = np.gradient (v, dx, dy, dz, axis=[1,0,2])
    dFz_dx, dFz_dy, dummy = np.gradient (w, dx, dy, dz, axis=[1,0,2])

    rot_x = dFz_dy - dFy_dz
    rot_y = dFx_dz - dFz_dx
    rot_z = dFy_dx - dFx_dy

    l = np.sqrt(np.power(u,2.0) + np.power(v,2.0) + np.power(w,2.0));

    m1 = np.multiply(rot_x,u)
    m2 = np.multiply(rot_y,v)
    m3 = np.multiply(rot_z,w)

    tmp1 = (m1 + m2 + m3)
    tmp2 = np.multiply(l,2.0)

    av = np.divide(tmp1, tmp2)

    return rot_x, rot_y, rot_z, av

- BBSysDyn