1
$\begingroup$

I have the following code to create a list of circles:

size = 7;
range = 30;
r = RandomReal[{.1, 6}, size];
loc = Table[
RandomReal[{-range + r[[i]], range - r[[i]]}, 2], {i, size}];

and then a function based on the distance from the other circles:

surf[x_, y_, j_] := -1/RegionDistance[
  RegionUnion[
    Table[
      Circle[loc[[i]], r[[i]]]
      , {i, Drop[Range[size], {j}]}
    ]
  ], {x, y}
]^2;

Now I want to get the partial derivatives at the location of each of the circles:

grad = Table[{D[surf[loc[[i, 1]], loc[[i, 2]], i], x], 
              D[surf[loc[[i, 1]], loc[[i, 2]], i], y]}, {i, size}]

but this always returns a list of {0,0}'s. Am I doing something wrong? Is this a bug? The function is clearly not constant everywhere:

enter image description here

$\endgroup$
  • 1
    $\begingroup$ The values are plugged into surf before the derivative is taken. Hence you're taking derivatives of constants. $\endgroup$ – Michael E2 Aug 9 '18 at 14:47
  • $\begingroup$ @MichaelE2 I believe I understand the problem now but still don't know the solution. I tried replacing the grad with grad[x_, y_, i_] := {D[surf[x, y, i], x], D[surf[x, y, i], y]} and then evaluating at grad[loc[[1, 1]], loc[[1, 2]], 1], but I get D(10.184)[-0.0025], D(-9.826)[-0.0025] as the output. ({10.1,-9.8} is loc[[1]]) $\endgroup$ – Sam Aug 9 '18 at 15:16
  • $\begingroup$ You probably want ClearAll[grad]; grad[x_, y_, i_] := {Derivative[1, 0, 0][surf][x, y, i], Derivative[0, 1, 0][surf][x, y, i]}. Be sure to clear grad. You also need clear surf and define it with ?NumericQ, surf[x_?NumericQ, y_?NumericQ, 1] :=..., since the expression seems not to be symbolically differentiable, but only numerically. $\endgroup$ – Michael E2 Aug 9 '18 at 15:35
  • $\begingroup$ @MichaelE2 Thanks for this, I was not aware of this very helpful notation for Derivative. I'm using Remove["Global*"] to clear everything at the beginning of the notebook, so I believe it's not necessary to clear grad again. $\endgroup$ – Sam Aug 9 '18 at 15:47
1
$\begingroup$

Here is my way of fixing up the OP's functions:

ClearAll[surf];
surf[x_?NumericQ, y_?NumericQ, j_] := -1/RegionDistance[RegionUnion @@
      Table[Circle[loc[[i]], r[[i]]], {i, Drop[Range[size], {j}]}], {x, y}]^2;

ClearAll[grad];
grad[x_, y_, i_] := {Derivative[1, 0, 0][surf][x, y, i], 
  Derivative[0, 1, 0][surf][x, y, i]}

And here's a speed improvement that factors the Regions because surf is very slow. Also, it uses an order-2 approximation to the derivative instead of the built-in order-8 method that Derivative uses when it cannot do the derivative symbolically.

ClearAll[surf2];
surf2[x_?NumericQ, y_?NumericQ, j_] := -1/Min[RegionDistance[#, {x, y}] & /@ 
      Table[Circle[loc[[i]], r[[i]]], {i, Drop[Range[size], {j}]}]]^2;

ClearAll[grad2];
grad2[x_, y_, i_] := 
  With[{dx = Max[Abs[x] Sqrt@$MachineEpsilon, $MachineEpsilon],
    dy = Max[Abs[y] Sqrt@$MachineEpsilon, $MachineEpsilon]},
   {(surf2[x + dx, y, i] - surf2[x - dx, y, i])/(2 dx), (
    surf2[x, y + dy, i] - surf2[x, y - dy, i])/(2 dy)}
   ];

And here's the difference in timing:

grad[loc[[1, 1]], loc[[1, 2]], 1] // AbsoluteTiming
(*  {75.8591, {-0.00866156, -0.00262315}}  *)

grad2[loc[[1, 1]], loc[[1, 2]], 1] // AbsoluteTiming
(*  {0.002797, {-0.00866156, -0.00262315}}  *)

The order-8 derivative used internally in grad[] with surf2 instead of surf, takes about 25 times longer than the last one (0.07 sec.). It may be sufficiently fast.

Update: I realized just before I had to leave that surf2[] could be symbolically differentiated if I removed the ?NumericQ. This is because RegionDistance evaluates to a symbolic expression on a simple circle but not on a union of them. So here's the fastest way:

ClearAll[surf3];
surf3[x_, y_, j_] := -1/
   Min[ComplexExpand@RegionDistance[#, {x, y}] & /@ 
      Table[Circle[loc[[i]], r[[i]]], {i, Drop[Range[size], {j}]}]]^2;

ClearAll[grad3];
Block[{x, y}, (* protects x,y since we pre-evaluate the derivatives *)
 grad3[x_, y_, 1] = {D[surf3[x, y, 1], x], D[surf3[x, y, 1], y]};
 ]

grad3[loc[[1, 1]], loc[[1, 2]], 1] // AbsoluteTiming
(*  {0.000364, {-0.00866156, -0.00262315}}  *)
$\endgroup$

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.