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, 2012 at 14:40
  • $\begingroup$ @ruebenko yup. good catch. fixed it. $\endgroup$
    – Eli Lansey
    Sep 6, 2012 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, 2012 at 14:43
  • 1
    $\begingroup$ From the help D[f, {{x1, x2, ...}}] gives the vector derivative $\endgroup$ Sep 6, 2012 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, 2012 at 16:10

3 Answers 3

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, 2012 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$ Sep 7, 2012 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, 2012 at 17:21
  • $\begingroup$ @Verde Whoops, you're right! Let me fix that. $\endgroup$
    – Eli Lansey
    Sep 6, 2012 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, 2012 at 18:54
  • $\begingroup$ @EliLansey Works here ... Mma v 8.0.0. Just in case, start a fresh kernel $\endgroup$ Sep 6, 2012 at 19:03
  • $\begingroup$ yup, that fixed it. $\endgroup$
    – Eli Lansey
    Sep 6, 2012 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, 2012 at 19:25
  • $\begingroup$ It takes 11 seconds to go through 20 points using this method. $\endgroup$
    – Eli Lansey
    Sep 6, 2012 at 19:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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