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

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.

-
Is there a typo in the last component of the curl, should that be D[#[[1]], y]? –  user21 Sep 6 '12 at 14:40
@ruebenko yup. good catch. fixed it. –  Eli Lansey Sep 6 '12 at 14:41
have you looked at Derivative for this. e.g. Derivative[0,0,1][intf][x,y,z]. (I have to go now, so I can not post a solution) –  user21 Sep 6 '12 at 14:43
From the help D[f, {{x1, x2, ...}}] gives the vector derivative –  belisarius Sep 6 '12 at 15:24
@belisarius J.M. pointed that out in chat, but I've been having trouble getting it to work with the interpolated functions. –  Eli Lansey Sep 6 '12 at 16:10

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]]

-
Nice. Using a Module fixes a lot of problems. Not sure why I didn't think of it. It's fast, too! –  Eli Lansey Sep 6 '12 at 21:50
+1 very elegant –  belisarius Sep 7 '12 at 3:43
An alternative that I mentioned to Eli in chat: for the curl, form the skew part jac - Transpose[jac] from the Jacobian returned by intfd[], and then extract the required entries to form the curl. –  Ｊ. Ｍ. Sep 7 '12 at 13:20

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} *)

-
Are you sure about the last component? The function depends only on z, and the z component of the curl contains derivatives in the other two coordinates ... I guess you didn't correct the third component after @ruebenko comment –  belisarius Sep 6 '12 at 17:21
@Verde Whoops, you're right! Let me fix that. –  Eli Lansey Sep 6 '12 at 17:35
(*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 }
*)

-
This looks promising. I tried it, but I got (*Part::partd: "Part specification 0.4[[3]] is longer than depth of object"*). –  Eli Lansey Sep 6 '12 at 18:54
@EliLansey Works here ... Mma v 8.0.0. Just in case, start a fresh kernel –  belisarius Sep 6 '12 at 19:03
yup, that fixed it. –  Eli Lansey Sep 6 '12 at 19:14
From a timing standpoint, this is way slower than my version. I can run through 10000 points in 11 seconds using mine, but I got tired of waiting after 5 minutes. Not sure why yours takes so long. –  Eli Lansey Sep 6 '12 at 19:25
It takes 11 seconds to go through 20 points using this method. –  Eli Lansey Sep 6 '12 at 19:27
show 1 more comment