Obliczenie zawinięcia pola wektorowego w Pythonie i narysowanie go za pomocą matplotlib

Question

Potrzebuję obliczyć zwijanie pola wektorowego i narysować go za pomocą matplotlib. Prosty przykład tego, czego szukam, może wyglądać następująco:

Jak mogę obliczyć i wykreślić zwijanie pola wektorowego w galerii quiver3d_demo.py w galerii matplotlib?

Answer 1

Możesz użyć sympy.curl(), aby obliczyć zwijanie pola wektorowego.

Przykład:

Załóżmy, że:
F = (Y oo, -xy, z) = Y z x - y xy + z z, następnie y będzie R[1], x jest R[0] i z jest R[2] natomiast wektory 3 osi byłoby R.x, R.y, R.z i kod obliczenie pola wektorowego wywinięcia to:

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) 

W przypadku G byłby równy R_y**2*R.y + (-2*R_y*R_z - R_y)*R.z lub, innymi słowy,
G = (0, y , -2yz-y).

Aby wykreślić to, należy przekonwertować powyższy wynik na 3 osobne funkcje; u, v, w.

(przykład poniżej zaadaptowany z matplotlib example on this link):

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

i ostateczny wynik jest taki:

Answer 2

Aby obliczyć Curl funkcji wektora można także użyć numdifftools dla automatyczne numeryczne różnicowanie bez objazdu poprzez symboliczne różnicowanie. Numdifftools nie zapewnia funkcji curl(), ale oblicza ona jakobianową macierz funkcji wektorowej jednej lub więcej zmiennych, co zapewnia pochodne wszystkich składników pola wektorowego w odniesieniu do wszystkich zmiennych; to wszystko, co jest konieczne do obliczenia zawinięcia.

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) 

ta zwraca wartość zwinięcie w x: array([-216., -2., 0.]) kreślenia jak zasugerowano powyżej.