3
$\begingroup$

The following code is used to compute the gradient of a given 2D data.

computeGrad[data_] := Module[{D1, D2, grad, h = 1, xdim, ydim, n},
   D1 = Append[data[[2 ;; -1, ;;]], data[[-1, ;;]]];
   D2 = Prepend[data[[1 ;; -2, ;;]], data[[1, ;;]]];
   {xdim, ydim} = Dimensions[data];
   grad = ConstantArray[0, {xdim, ydim, 2}];
   (*central differences on X-axis*)
   grad[[;; , ;; , 1]] = (D1 - D2)/(2*h);
   grad[[1, ;; , 
     1]] = (4*data[[2, ;;]] - 3*data[[1, ;;]] - data[[3, ;;]])/(2*h);
   grad[[-1, ;; , 
     1]] = -(4*data[[-2, ;;]] - 3*data[[-1, ;;]] - data[[-3, ;;]])/(2*
       h);
   (*central differences on Y-axis*)
   D1 = Join[data[[;; , 2 ;; -1]], Transpose[{data[[;; , -1]]}], 2];
   D2 = Join[Transpose[{data[[;; , -1]]}], data[[;; , 1 ;; -2]], 2];
   grad[[;; , ;; , 2]] = (D1 - D2)/(2*h);
   grad[[;; , 1, 
      2]] = (4*data[[;; , 2]] - 3*data[[;; , 1]] - data[[;; , 3]])/(2*
       h);
   grad[[;; , -1, 
      2]] = -(4*data[[;; , -2]] - 3*data[[;; , -1]] - 
        data[[;; , -3]])/(2*h);
   n = ArrayDepth[grad];
   Map[Normalize, grad, {n - 1}] // N
   ];

For example:

data = Table[Sin[x]*Cos[y], {x, 0, 2 Pi, 0.2}, {y, 0, 2 Pi, 0.2}];
grad = computeGrad[data];

gives enter image description here How to extend the code to deal with data higher than 2 dimensions? Is there any place to enhance the code to make it more efficient?

$\endgroup$
  • 1
    $\begingroup$ Use ListInterpolation + Grad. $\endgroup$ – Kuba Apr 17 '15 at 11:04
7
$\begingroup$

There are different ways to do this. The easiest, I think, is to use DerivativeFilter:

data = Table[Sin[x]*Cos[y], {x, 0, 2 Pi, 0.2}, {y, 0, 2 Pi, 0.2}];

gradFilter[data_] := Module[{n = ArrayDepth[data]},
  MapThread[List, 
   Table[DerivativeFilter[data, UnitVector[n, i]], {i, n}], n]]

ListVectorPlot[gradFilter[data]]

vectorplot

I made this independent of dimensions by constructing the entries in the gradient from derivatives in the direction of each UnitVector in the range given by the value of ArrayDepth for the data.

DerivativeFilter uses interpolation of the data, and has an option for InterpolationOrder. But if you want to work without interpolation, there is an alternative way of computing gradients that uses finite difference methods, similar to what you programmed in the question.

The following function allows you to set the order of the finite-difference approximation (the number of points for the discrete derivative) explicitly, as a second argument order:

gradFilter[data_, order_] := 
 Module[{n = ArrayDepth[data]}, 
  MapThread[List, 
   Table[NDSolve`FiniteDifferenceDerivative[UnitVector[n, i], 
      Range /@ Dimensions[data], "DifferenceOrder" -> order][
     data], {i, n}], n]]

This uses the built-in functionality explained in the documentation for The Numerical Method of Lines to construct the derivative operator for each direction, and then applies it to the data.

Here is a comparison between the exact gradient and the one calculated with the second function, using derivative order 2:

v[x_, y_, z_] := 
 1/Sqrt[x^2 + y^2 + z^2 + .1] - 1/
  Sqrt[.1 + (x - 1)^2 + (y - 1)^2 + (z - 1)^2]

e[x_, y_, z_] = D[v[x, y, z], {{x, y, z}}];

gradCompare = 
  Table[e[x, y, z], {x, -2, 2, .1}, {y, -2, 2, .1}, {z, -2, 2, .1}];

ListVectorPlot3D[gradCompare, AxesLabel -> {x, y, z}, 
 VectorColorFunction -> (Hue[#7] &)]

original

data = Table[
   v[x, y, z], {x, -2, 2, .1}, {y, -2, 2, .1}, {z, -2, 2, .1}];

gr2 = gradFilter[data, 2];

ListVectorPlot3D[gr2, AxesLabel -> {x, y, z}, 
 VectorColorFunction -> (Hue[#7] &)]

discrete

So the second-order derivative approximation was sufficient to reproduce the gradient in this case. The first version of the gradientFilter without a specified order requires an additional interpolation run for each directional derivative, and therefore is noticeably slower than the implementation using FiniteDifferenceDerivative - I didn't show the plot for that result because it looks visually the same as the other ones.

Just remember that you have to specify the derivative order explicitly to invoke the faster version of the function.

$\endgroup$
  • $\begingroup$ `@Jens,Very general and comprehensive. I think the 2 in the first gradFilter function should be n. The code I provided is a little more efficient for 2 dimensional data though. :-) $\endgroup$ – novice Apr 18 '15 at 11:16
  • $\begingroup$ For the first function I think the bottom two levels need to be transposed - at least, I need Transpose[gradFilter[data], {1, 2, 4, 3}] for it to work with ListVectorPlot3D and 3D data. For the second function, I think instead of Partition[..., Length[data]] you need ArrayReshape[..., Append[Dimensions[data], n]]. $\endgroup$ – Michael E2 Apr 18 '15 at 12:14
  • $\begingroup$ @MichaelE2 Oh sorry - I din't notice that. I'll double check and make the change. Thanks for pointing it out! $\endgroup$ – Jens Apr 18 '15 at 14:49
  • $\begingroup$ @novice Thanks for the correction, I edited the answer to also include the comment by MichaelE2. $\endgroup$ – Jens Apr 18 '15 at 15:27
  • $\begingroup$ +1 on the thorough answer, particularly note on specifying order. $\endgroup$ – ciao Apr 20 '15 at 6:07

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.