I have some data defined in an array (an image) and I need to find the curl of a certain function.

Wikipedia has an integral definition of curl that I like, maybe it can be discrete.

$$ \nabla \times F = \frac{1}{|A|} \oint F \cdot d \mathbf{r}$$

Can we write this as some kind of finite difference. Let $F = (F_1, F_2, F_3)$. I can think of two different ways to write curl:

$$ \frac{1}{(2\epsilon)^2} \big[ F_2(x+\epsilon, y, z) + F_1(x, y+\epsilon, z) - F_2(x-\epsilon, y, z) + F_1(x, y-\epsilon, z) \Big] $$

Hopefully I wrote that average correctly. What if I used the corners of the square instead? Do we get curl still?

$$ \frac{1}{(2\epsilon)^2}\Big[ F(x+\epsilon, y + \epsilon, z) + F(x+\epsilon, y - \epsilon, z) - F(x-\epsilon, y - \epsilon, z) - F(x-\epsilon, y + \epsilon, z) \Big] $$

I don't know what to call these two operators $\nabla_1 $ and $\nabla_2$ who do these operators differ as $\epsilon \ll 1$ ?

My second definition of curl breaks down a little bit. However all the points $(x,y,z) + \epsilon \,(\pm 1, \pm 1, 0)$ all lie on a circle.

$$ \frac{1}{ \pi \epsilon^2} \sum_{\theta = \{ \frac{2\pi k}{n}: 0 \leq k < n \}} F\big[ (x,y,z) + \epsilon( \cos \theta, \sin \theta, 0 )\big] \cdot (-\sin \theta, \cos \theta, 0) $$

This operator should be proportional to curl for small $\epsilon \approx 0$, I think, but this is harder to put into a computer.

  • 2
    What if $f$? Curl is applied to a vector field, not a scalar $f$. The first formula is $\frac{f}{2\epsilon^2}+O(\epsilon^{-1})$, the second $\frac{f_x}{\epsilon}+O(\epsilon)$. Both diverge as $\epsilon\to0$ and neither is related to curl. – Kirill Jan 26 '16 at 16:08
  • 2
    Have a look at the paper of Hyman and Shashkov - you'll find a comprehensive set of discretizations of vector fields (including curl, on page 11) – GoHokies Jan 26 '16 at 16:18
  • @Kirill I made a mistake - or possibly dug my self deeper now. $F$ should be a vector. – john mangual Jan 26 '16 at 16:21
  • If $F$ is a vector, then $\nabla\times F$ also needs to be a vector. But your formulas only show a single element!? – Wolfgang Bangerth Jan 28 '16 at 17:28
  • The lhs of the first equation should be $(\nabla \times F) \cdot n $ where $n$ is a normal vector. – Biswajit Banerjee Jan 28 '16 at 23:49


It looks like you only need 2d curl, so let's start with a simpler continuous definition:

$$ \omega = \nabla \times \mathbf{u} = \frac{\delta v}{\delta x} - \frac{\delta u}{\delta y} $$

where 2d vector field $\mathbf{u}=(u,v)$ (same as your $\mathbf F = (F_1,F_2,F_3)$, dropping $F_3$). Note that curl is a vector and in the 2d version, and is purely in the z direction, and is customarily written as a scalar.


We'll assume $(u,v)$ is sampled on a grid (i.e. at "vertices", where lines cross, and indices have integer values), with spacing $\Delta x$ and $\Delta y$ (equal for a square grid, but helps "typecheck" the discretization).

Forward differences

The simplest Finite Difference discretization method is forward differences: $$ \omega_{i,j} = \frac{ v_{i+1,j} - v_{i,j} }{ \Delta x } - \frac{ u_{i,j+1}-u_{i,j} }{ \Delta y } $$

Central differences

But central differences is more accurate (note the 2 in the denominator, because we're discretizing over 2 grid cells):

$$ \omega_{i,j} = \frac{ v_{i+1,j} - v_{i-1,j} }{ 2\Delta x } - \frac{ u_{i,j+1}-u_{i,j-1} }{ 2\Delta y } $$

Finer central differences

If you look again at the first forward differences method, you'll see it's the same as central differences, just for $\omega_{i+\frac{1}{2},j+\frac{1}{2}}$ instead of $\omega_{ i,j}$. EDIT Actually not, because the central differences are for different locations: the first for $\omega_{i+\frac{1}{2},j}$ and the second for $\omega_{i,j+\frac{1}{2}}$. Instead, for each term, we can average the other co-ordinate values to interpolate the center. That is,

$$ \omega_{i+\frac{1}{2},j+\frac{1}{2}} = \frac{ \frac{ v_{i+1,j} - v_{i,j} }{ \Delta x } + \frac{ v_{i+1,j+1} - v_{i,j+1} }{ \Delta x } }{2} - \frac{ \frac{ u_{i,j+1}-u_{i,j} }{ \Delta y } + \frac{ u_{i+1,j+1}-u_{i+1,j} }{ \Delta y } }{2} $$ At this point, the advantages of the Calculus of Finite Differences (in Philip Roe's answer) become stark. The above is: $$ \omega_{i+\frac{1}{2},j+\frac{1}{2}} = \mu_y \delta_x v_{i,j} - \mu_x \delta_y u_{i,j} $$

Let me unpack that first term: the $\mu_y$ averages over values at $j$ and $j+1$ in y-direction (i.e. the two expressions divided by $2$), the $\delta_x$ is the difference between values at $i$ and $i+1$ in x-direction (i.e. each of the aforementioned expressions, which are divided by $\Delta x$). I'm not sure what happens to the $\Delta x$; I think the idea is to let $\Delta x=1$, so it can be ignored.

Or (if you're better than me at remembering how the offset indices go) just: $$ \omega = \mu_y \delta_x v - \mu_x \delta_y u $$

These seem what you're going for with your $\nabla_1$ and $\nabla_2$.

EDIT Now I think you're going for finite differences along diagonal lines for your $\nabla_2$ (the corners of a square). IDK, but I think that would work, but you'd need to calculate $(u,v)$ perpendicular to each diagonal, and also adjust the lengths that you're dividing by (though I guess the wikipedia integral version may take care of this, since it uses area, not lengths). Of course, like the above, it would be for the center of a cell, $\omega_{ i+\frac{1}{2},j+\frac{1}{2} }$, rather than at the data points.


I won't go into it here, but the accuracy of the above Finite Difference methods is usually analysed with Taylor Series. Forward differences is first order $O(\Delta x)$, central differences is second order $O( \Delta x^2 )$.

This is made a lot easier by introducing the Calculus of Finite Differences.

If $u_{i,j}$ is grid function defined for integer $i,j$ then the y-differencing operator $$\delta_y u=u_{j+1}-u_j$$ is defined at $i,j+1/2$, so no need to write indices. Similarly, $$\mu_xu=(u_{i,j}+u_{i+1,j})/2 $$ is the x-averaging operator and is defined at $i+1/2,j$. We complete the list of operators with $\delta_x$ and $\mu_y$ These operators are binary and the output is defined halfway between the inputs. In the examples so far, they are on the edges of the cells that have primary quantities at vertices.

The operators can be combined according to most laws of algebra, added, multiplied, or raised to powers. For example $\mu_x\delta_x$ is a central difference on the original grid, and $\mu_y\delta_x$ approximates an x-derivative at a cell center. An important rule however, is that two operators can only be equated or combined if their outputs are at the same location.

A nice discrete representation of curl in the 2D case is $$\mu_y\delta_xv-\mu_x\delta_yu$$ at the cell centers if the data were at vertices or vice versa.

You could also use $$\mu_x\delta_xv-\mu_y\delta_yu$$ A drawback of this this formula is that if you imagine the array of data to be colored like a chessboard, it only uses data of one color. In some circumstances this leads to trouble.

  • Some typos: $ around u_{i,j}; an = in the $\mu_y$ definition(?); "You can" (third last para). – hyperpallium Jul 31 '17 at 8:16
  • 1
    Thanks for the typos! Imagine you have a chessboard. You can define things at the centers of the squares $i,j$ or at the vertical interfaces $i+1/2,j$ or the horizontal interfaces $i,J+1/2$ or at the vertices $I+1/2,j+1/2$. The operators are binary, like a+b. The inputs are always the same (cell, interface,vertex) and the output is halfway between them. The operators commute, which is very powerful. Every operator contracts the domain by one cell size in one direction. – Philip Roe Jul 31 '17 at 9:24
  • 1
    The operator $\mu_x\mu_y$ averages four cells meeting at a vertex, or four vertices forming a cell, or four horizontal edges surrounding one vertical edge. The standard five-point Laplacian is just $\delta_x^2+\delta_y^2.$ The nine-point Laplacian is $\mu_y^2\delta_x^2+\mu_x^2\delta_y^2$. there are rules like $$\delta(uv)=\mu u \delta v+\mu v \delta u$$ – Philip Roe Jul 31 '17 at 9:26
  • 1
    Just try writing out few FD expressions that you are familiar with and see how much simpler they are. The thing that takes a bit of getting used to is how the operations change the locations, but they must do so consistently. There is also a very nice tie-in to Fourier analysis. $\mu_x$ multiplies everything by $\cos{(\theta_x/2)}$ and $\delta x$ by $2i\sin{(\theta_x/2})$ – Philip Roe Jul 31 '17 at 13:05
  • 1
    @hyperpallium What you define in your answer is not $\omega_{I,j}$ It IS $\omega{1+1/2,,j+1/2}$ With data in cells, $\delta_x$ puts the output on vertical faces and $\mu_y$ takes it to a vertex. Then when you go to shorthand, I'm afraid you give expressions for the divergence. :( Needs an edit. You can also use what are called divided differences where division by the mesh size is included in the differencing operators. These are usually not written differently, but an author must explain his convention. Works fine if $\Delta x$ and $\Delta y$ are both constant, but not otherwise. – Philip Roe Aug 7 '17 at 18:18

Your Answer


By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.