Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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
1  
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

3 Answers 3

up vote 7 down vote accepted

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]]
share|improve this answer
    
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
4  
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. –  J. M. Sep 7 '12 at 13:20
(*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 }
*)
share|improve this answer
    
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

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} *)
share|improve this answer
    
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 Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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