12
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ Is there a typo in the last component of the curl, should that be D[#[[1]], y]? $\endgroup$
    – user21
    Sep 6 '12 at 14:40
  • $\begingroup$ @ruebenko yup. good catch. fixed it. $\endgroup$
    – Eli Lansey
    Sep 6 '12 at 14:41
  • $\begingroup$ 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) $\endgroup$
    – user21
    Sep 6 '12 at 14:43
  • 1
    $\begingroup$ From the help D[f, {{x1, x2, ...}}] gives the vector derivative $\endgroup$ Sep 6 '12 at 15:24
  • $\begingroup$ @belisarius J.M. pointed that out in chat, but I've been having trouble getting it to work with the interpolated functions. $\endgroup$
    – Eli Lansey
    Sep 6 '12 at 16:10
8
$\begingroup$

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]]
$\endgroup$
2
  • $\begingroup$ Nice. Using a Module fixes a lot of problems. Not sure why I didn't think of it. It's fast, too! $\endgroup$
    – Eli Lansey
    Sep 6 '12 at 21:50
  • 4
    $\begingroup$ 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. $\endgroup$
    – J. M.'s torpor
    Sep 7 '12 at 13:20
6
$\begingroup$

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} *)
$\endgroup$
2
  • $\begingroup$ 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 $\endgroup$ Sep 6 '12 at 17:21
  • $\begingroup$ @Verde Whoops, you're right! Let me fix that. $\endgroup$
    – Eli Lansey
    Sep 6 '12 at 17:35
4
$\begingroup$
(*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 }
*)
$\endgroup$
6
  • $\begingroup$ This looks promising. I tried it, but I got (*Part::partd: "Part specification 0.4[[3]] is longer than depth of object"*). $\endgroup$
    – Eli Lansey
    Sep 6 '12 at 18:54
  • $\begingroup$ @EliLansey Works here ... Mma v 8.0.0. Just in case, start a fresh kernel $\endgroup$ Sep 6 '12 at 19:03
  • $\begingroup$ yup, that fixed it. $\endgroup$
    – Eli Lansey
    Sep 6 '12 at 19:14
  • $\begingroup$ 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. $\endgroup$
    – Eli Lansey
    Sep 6 '12 at 19:25
  • $\begingroup$ It takes 11 seconds to go through 20 points using this method. $\endgroup$
    – Eli Lansey
    Sep 6 '12 at 19:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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