Мне нужно вычислить завиток векторного поля и построить его с помощью matplotlib. Простой пример того, что я ищу, можно было бы сделать так:
Как я могу вычислить и построить завиток векторного поля в quiver3d_demo.py в галерее matplotlib?
Мне нужно вычислить завиток векторного поля и построить его с помощью matplotlib. Простой пример того, что я ищу, можно было бы сделать так:
Как я могу вычислить и построить завиток векторного поля в quiver3d_demo.py в галерее matplotlib?
Вы можете использовать sympy.curl()
для вычисления скручивания векторного поля.
Пример:
Предположим, F(x, y, z) = y 2 z i - xy j + z 2k, тогда:
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 = 0 i + y 2j + (-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()
И окончательный результат таков:
Чтобы вычислить завиток векторной функции, вы также можете использовать 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.])
Вычисление происходит так, как указано выше.