用Python计算向量场的旋度并用matplotlib绘制

Question

我需要计算向量场的旋度，并用matplotlib绘制出来。我要找的一个简单的例子可以这样说：

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

Answer 1

可以使用^{}计算向量场的旋度。

示例：

假设F（x，y，z）=y²zi-xyj+z²k，那么：

• y是R[1]，x是R[0]，z是R[2]
• 三个轴的单位向量分别为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

（下面的示例改编自thismatplotlib示例）：

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()

最终的结果是：

Answer 2

要计算向量函数的旋度，还可以使用numdifftools进行自动数值微分，而无需绕道符号微分。Numdifftools不提供curl()函数，但它计算一个或多个变量的向量值函数的雅可比矩阵，这提供了向量场的所有分量相对于所有变量的导数；这是计算旋度所必需的。

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：array([-216., -2., 0.])的值 绘制如上文所述。