Solución:
Puedes usar sympy.curl()
para calcular la curvatura de un campo vectorial.
Ejemplo:
Suponer F(x, y, z) = y2zI – xyj + z2k, luego:
-
y
seríaR[1]
,x
esR[0]
yz
esR[2]
- los vectores unitarios I, j, k de los 3 ejes, serían respectivamente
R.x
,R.y
,R.z
.
El código para calcular el rizo del campo vectorial es:
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)
En ese caso, G sería igual a R_y**2*R.y + (-2*R_y*R_z - R_y)*R.z
o, en otras palabras,
G = 0I + y2j + (-2yz-y)k.
Para trazarlo necesita convertir el resultado anterior en 3 funciones separadas; u, v, w.
(ejemplo a continuación adaptado de este ejemplo de 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()
Y el resultado final es este:
Para calcular la curvatura de una función vectorial, también puede utilizar numdifftools para la diferenciación numérica automática sin un desvío a través de la diferenciación simbólica. Numdifftools no proporciona un curl()
función, pero calcula la matriz jacobiana de una función de valor vectorial de una o más variables, y esto proporciona las derivadas de todos los componentes de un campo vectorial con respecto a todas las variables; esto es todo lo que se necesita para calcular el rizo.
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)
Esto devuelve el valor del rizo en x
: array([-216., -2., 0.])
El trazado es como se sugirió anteriormente.
Aquí hay un código de Python que se basa en una implementación de Octave / Matlab,
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