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.

Improve this question
$\endgroup$
7
  • $\begingroup$ Is there a typo in the last component of the curl, should that be D[#[[1]], y]? $\endgroup$ Commented Sep 6, 2012 at 14:40
  • $\begingroup$ @ruebenko yup. good catch. fixed it. $\endgroup$ Commented 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$ Commented Sep 6, 2012 at 14:43
  • 1
    $\begingroup$ From the help D[f, {{x1, x2, ...}}] gives the vector derivative $\endgroup$ Commented 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$ Commented Sep 6, 2012 at 16:10

3 Answers 3

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.