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

Question

我需要计算矢量场的旋度并使用 matplotlib 绘制它。我正在寻找的一个简单示例可以这样说：

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

Answer 1

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

示例：

假设 F(x,y,z) = y²zi - xyj + z²k, 那么:

y 是 R[1]，x 是 R[0]，z 是 R[2]
3 个轴的单位向量 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.

（以下示例改编自 this 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()

最后的结果是这样的：

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)

returns x 处的卷曲值：array([-216., -2., 0.]) 绘图如上所示。

Answer 3

这是一个基于 Octave / Matlab implementation,

的 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

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

Calculate curl of a vector field in Python and plot it with matplotlib

python

matplotlib

sympy