Подтвердить что ты не робот

Вычислить завиток векторного поля в Python и построить его с помощью matplotlib

Мне нужно вычислить завиток векторного поля и построить его с помощью matplotlib. Простой пример того, что я ищу, можно было бы сделать так:

Как я могу вычислить и построить завиток векторного поля в quiver3d_demo.py в галерее matplotlib?

4b9b3361

Ответ 1

Вы можете использовать sympy.curl() для вычисления скручивания векторного поля.

Пример:

Предположим, F(x, y, z) = y 2 z i - xy j + z 2k, тогда:

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

И окончательный результат таков:

enter image description here

Ответ 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.]) Вычисление происходит так, как указано выше.