Derivatives (divergence, gradient, curl) of interpolated 3D data

Question

How can I calculate a vector derivative (i.e. divergence, gradient, curl) of interpolated data? For sample data, you can use:

f[x_, y_, z_] := Exp[I z] {1, 0, 0}
testdata=Flatten[Table[N@{x,y,z,f[x,y,z]},{x,0,4 Pi,Pi/10},{y,0,4 Pi,Pi/10},{z,0,4 Pi,Pi/10}],2];
intf = Interpolation[testdata]

I know that for 1D data, like

dim1 = Table[N@{z, f[0, 0, z]}, {z, 0, 2 Pi, Pi/10}];
int1 = Interpolation@dim1;

you can do D[int1[x],x]. However, I can't seem to get convince MMA that the interpolated function intf actually returns a "vector" quantity (i.e. Length@intf[x,y,z]->3), so that I can do something like:

curl = {D[#[[3]],y]-D[#[[2]],z],-(D[#[[3]],x]-D[#[[1]],z]),D[#[[2]],x]-D[#[[1]],y]}&;
curl@intf[x, y, z]
(* Out[] := {0,0,0} *)

Similarly, this should work for gradient:

grad = {D[#[[1]], x], D[#[[2]], y], D[#[[3]], z]} &;

and divergence:

div = (D[#[[1]], x] + D[#[[2]], y] + D[#[[3]], z])&;

I've posted my best attempt in an answer below, but I'm wondering what other solution approaches there are.

Answer 1

Here's another possibility. First compute the derivatives (this should use formal symbols for safety, but I've used x,y,z here to avoid cluttering the code)

intfd[x_, y_, z_] = D[intf[x, y, z], {{x, y, z}}];

then define the curl by

intfcurl[x_, y_, z_] := Module[{q = intfd[x, y, z]},
  {q[[2, 3]] - q[[3, 2]], q[[3, 1]] - q[[1, 3]], q[[1, 2]] - q[[2, 1]]}]

The results seem okay:

intfcurl[0.4, 0.2, 0.3]
(*  {0., -0.297934 + 0.954217 I, 0. + 0. I}  *)

and the timing is reasonable:

intfcurl @@@ testcoords; // Timing
(*  {7.863, Null}  *)

---

You can similarly define the gradient and the divergence using intfd:

intfgrad[x_, y_, z_] := Module[{q = intfd[x, y, z]}, Diagonal[q]]
intfdiv[x_, y_, z_] := Module[{q = intfd[x, y, z]}, Tr[q]]

Answer 2

The best I've been able to come up with is to separately interpolate each component of the vector data:

testdatacomps =Table[{Sequence @@ #[[1 ;; 3]], #[[4, n]]} & /@ testdata, {n, 1, 3}];
intcomps = Interpolation /@ testdatacomps;

Then

interpolatedcurl[xx_, yy_, zz_] := 
 curl@Through[intcomps[x, y, z]] /. Thread[{x, y, z} -> {xx, yy, zz}]

This result compares well to the actual curl:

curl@f[x, y, z] /. Thread[{x, y, z} -> {.4, .2, .3}]
interpolatedcurl[.4, .2, .3]
(* Out[1] := {0, -0.29552 + 0.955336 I, 0}
   Out[2] := {0., -0.297934 + 0.954217 I, 0. + 0. I} *)

Edit It is also reasonably fast. Evaluating at 10000 random coordinates:

testcoords = RandomReal[{0, 4 \[Pi]}, {10000, 3}];

interpolatedcurl @@@ testcoords; // Timing
(* Out[] := {11.451, Null} *)

Answer 3

(*Your definitions*)
f[x_, y_, z_] := Exp[I z] {1, 0, 0}
testdata =Flatten[Table[N@{x,y,z,f[x,y,z]},{x,0,Pi,Pi/10},{y,0,Pi,Pi/10},{z,0,Pi, Pi/10}],2];
intf = Interpolation[testdata]


(* A derivative operator *)

d[i_Integer][args__] :=(Derivative[Sequence@@ RotateRight[{1, 0, 0}, i - 1]] [intf])[x, y, z] 
                       /. Thread[{x, y, z} -> {args}]

(*Curl*)

curl[a__] := Chop/@ {d[2][a][[3]] - d[3][a][[2]], 
                     d[3][a][[1]] - d[1][a][[3]], 
                     d[1][a][[2]] - d[2][a][[1]]}

(*Test*)
curl[.4, .2, .3]
(* 
{0, -0.297934 + 0.954217 I, 0 }
*)