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

Question

I need to calculate the curl of a vector field and plot it with matplotlib. A simple example of what I am looking for could be put like that:

How can I calculate and plot the curl of the vector field in the quiver3d_demo.py in the matplotlib gallery?

There doesn't seem to be a built-in function in numpy or scipy to calculate curl. If so you'll have to write that; in 3D the result will also be a vector field so matplotlib plots it just as in the example. [Similar question about divergence](http://stackoverflow.com/questions/11435809/compute-divergence-of-vector-field-using-python) — cphlewis, May 21 '15 at 17:36
Thanks for your reply. Ok, I got it how to plot. Do you have any suggestions about a good way of writing the curl function? — celestos, May 22 '15 at 20:28
Maybe using Sympy (http://docs.sympy.org/dev/modules/physics/vector/fields.html) could be useful. It has some built-in functionality for vector fields. — Dietrich, May 24 '15 at 12:10
I see that Sympy has a `curl` function, and also handles plotting (using matplotlib as the back end), but I've never used it. — cphlewis, Jun 04 '15 at 19:03

user · Accepted Answer · 2019-08-08T07:48:12.557

You can use sympy.curl() to calculate the curl of a vector field.

Example:

Suppose F(x,y,z) = y²zi - xyj + z²k, then:

y would be R[1], x is R[0] and z is R[2]
the unit vectors i, j, k of the 3 axes, would be respectively R.x, R.y, R.z.

The code to calculate the vector field curl is:

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)

In that case G would be equal to R_y**2*R.y + (-2*R_y*R_z - R_y)*R.z or, in other words,
G = 0i + y²j + (-2yz-y)k.

To plot it you need to convert the above result into 3 separate functions; u,v,w.

(example below adapted from this matplotlib example):

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

And the final result is this:

enter image description here

Please explain the downvotes. I'd love to fix any mistakes. – user Aug 07 '19 at 19:09 — user, Aug 07 '19 at 19:09

score 5 · Answer 2 · edited Jan 01 '17 at 16:02

To calculate the curl of a vector function you can also use numdifftools for automatic numerical differentiation without a detour through symbolic differentiation. Numdifftools doesn't provide a curl() function, but it does compute the Jacobian matrix of a vector valued function of one or more variables, and this provides the derivatives of all components of a vector field with respect to all of the variables; this is all that's necessary for the calculation of the 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)

This returns the value of the curl at x: array([-216., -2., 0.]) Plotting is as suggested above.

I prefer this solution to the accepted answer, as it's infrequent that I deal with analytical vector fields. However, it is a bit inefficient as the curl does not use the trace of the Jacobian (thus 1/3 of the work in nd.Jacobian is not used here). Do you think it would be faster to compute the off-diagonal elements of the Jacobian independently using nd.directionaldiff()? Or do you think the additional python code overhead would make it not worth it? — NLi10Me, Jun 10 '19 at 17:33

score 5 · Answer 3 · answered Sep 22 '20 at 18:52

Here is a Python code that is based on an Octave / Matlab implementation,

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

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

3 Answers3

Linked