Speed up derivative evaluation

Question

I'm trying to calculate the normal vectors and the tangent vectors at the discrete points of a suface. For example:

f[x_, y_] := (Cos[x] Sin[y])/2. (* The surface.*)
fx[a_, b_] := 
  Module[{x, y}, D[f[x, y], x] /. {x -> a, y -> b}]
fy[a_, b_] := 
  Module[{x, y}, D[f[x, y], y] /. {x -> a, y -> b}]
normalVector[a_, b_] := 
  Module[{x, y}, 
   {-D[f[x, y], x], -D[f[x, y], y], 1} /. {x -> a, y -> b}]

tangentVector[x_, y_, θ_] := 
 {Cos[θ], Sin[θ], 
   fx[x, y] Cos[θ] + fy[x, y] Sin[θ]};

xr = yr = 1;
n = 20;
hx = 2 xr/n; 
hy = 2 yr/n; 
hθ = 2 Pi/n;
AbsoluteTiming[
 Table[
    nV = normalVector[x, y];
    tV = tangentVector[x, y, θ];
    {nV, tV}, 
    {x, -xr, xr - hx, hx}, {y, -yr, yr - hy, hy},
    {θ, -Pi, Pi - hθ, hθ}
    ];
 ]

gives (vectors shown only at one point)

The problem is that the calculation is slow. The desired discrete point should be 200 per axis. I think the time is mostly consumed by derivative evaluation. So how could I speed up the calculation ? One choice is to transform the code into C++. However the surface function might come from interpolation. I want to stay in Mathematica.

Answer 1

As simple as it sound, use Set instead of delayed assignment. i.e., replace your := in the function definitions with =.

First, make sure everything is safe by clearing x,y,a,b as:

ClearAll[x,y,a,b]

Then, define your functions as:

f[x_, y_] = (Cos[x] Sin[y])/2.; (*The surface.*)

fx[a_, b_] = Module[{x, y}, D[f[x, y], x] /. {x -> a, y -> b}];
fy[a_, b_] = Module[{x, y}, D[f[x, y], y] /. {x -> a, y -> b}];
normalVector[a_, b_] = 
  Module[{x, 
    y}, {-D[f[x, y], x], -D[f[x, y], y], 1} /. {x -> a, y -> b}];

I got 5x improvement in speed just by doing that.

This works because the function f[x,y] is fully defined and all the derivatives can be obtained symbolically beforehand. What is happening with the delayed assignment, is basically having D[f[x,y],x] being calculated each time a call is made for fx[a,b] is made. Repetitive evaluation get cashed, but apparently still not good enough in this case.

Also, I am not sure about your C++ comment, but a great idea would be to compile your code from Mathematica. Define f[x,y] as above and then use:

f[x_, y_] = (Cos[x] Sin[y])/2.; (*The surface.*)   
fx = Compile[{{a, _Real}, {b, _Real}}, 
   Module[{x, y}, D[f[x, y], x] /. {x -> a, y -> b}]];
fy = Compile[{{a, _Real}, {b, _Real}}, 
   Module[{x, y}, D[f[x, y], y] /. {x -> a, y -> b}]];

Update: After experimenting with Compile, there is no performance gain obtained. Only the use of Set (=) instead of the delayed assignment (:=).

Answer 2

Here's one speed-up: Use approximate machine real inputs. Outer is a bit faster than Table here.

Clear[x, y];
f[x_, y_] := (Cos[x] Sin[y])/2 (*The surface.*)
df[x_, y_] = D[f[x, y], {{x, y}}];

normalVector[a_, b_] := Join[-df[a, b], {1}];
tangentVector[x_, y_, θ_] := Join[#, {df[x, y].#}] &@{Cos[θ], Sin[θ]};

xr = yr = 1;
n = 20.;        (* < this makes the Range's below machine reals *)
hx = 2 xr/n;
hy = 2 yr/n;
hθ = 2 Pi/n;
RepeatedTiming[
 foo = Outer[{normalVector[#1, #2], tangentVector[##]} &, 
    Range[-xr, xr - hx, hx],
    Range[-yr, yr - hy, hy], 
    Range[-Pi, Pi - hθ, hθ]];]
(*   {0.11, Null}   *)

If you value speed over readability, here's this, which takes some advantage of the vectorization of Mathematica functions. More savings may be possible.

RepeatedTiming[
 glurg = Outer[
    Transpose@{
       ConstantArray[normalVector[#1, #2], Round[2 Pi/hθ]], 
       Transpose @ tangentVector[#1, #2, Range[-Pi, Pi - hθ, hθ]]} &, 
    Range[-xr, xr - hx, hx], Range[-yr, yr - hy, hy]];]
(*   {0.019, Null}   *)

foo == glurg
(*  True  *)

Answer 3

Just for fun:

f[x_, y_] := Cos[x] Sin[y]/2
grd[u_, v_] := Grad[z - f[x, y], {x, y, z}] /. {x -> u, y -> v};
tg[x0_, y0_, {a_, b_}] := 
 With[{ru = {1, 0, D[f[x, y], x] /. {x -> x0, y -> y0}}, 
   rv = {0, 1, D[f[x, y], y] /. {x -> x0, y -> y0}}}, a ru + b rv]
tn[x_, y_, {s_, t_}] := 
 With[{pt = {x, y, f[x, y]}}, 
  Graphics3D[{Point[pt], PointSize[0.04], Arrowheads[0.02], Red, 
    Arrow[{pt, pt + Normalize[grd[x, y]]}], Blue, 
    Arrow[{pt, pt + Normalize[tg[x, y, {s, t}]]}]}]]
tab = Flatten[
   Table[Show[
     Plot3D[f[x, y], {x, 0, 2 Pi}, {y, 0, 2 Pi}, 
      PlotRange -> {-1, 1.5}, PlotStyle -> Opacity[0.5], 
      ColorFunction -> "DarkRainbow"], 
     tn[#, #, {Cos[t], Sin[t]}] &@pnt], {pnt, 0, 2 Pi, 0.3}, {t, 0, 
     2 Pi, 0.3}], 1];

An individual plot does not take long. The following gif has been cut down from tab to display: