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;
    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)

enter image description here

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.

  • $\begingroup$ You know that the gradient can be computed in one go, right? D[f[x, y], {{x, y}}] $\endgroup$ – J. M.'s ennui May 26 '15 at 4:04
  • $\begingroup$ Nice tip for speedup. $\endgroup$ – novice May 26 '15 at 4:09

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:


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_] = 
    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 (:=).

  • $\begingroup$ `@Bichoy, great improvement ! another 15x speedup. Does compile support interpolation functions? $\endgroup$ – novice May 26 '15 at 3:07
  • $\begingroup$ I am not sure :) You have to try ... You should also read the documentation for Compile, you can set the compilation target to "C", set the RuntimeAttributes to {Listable} and enable parallelization for even much more improvement... Just check Compile documentation. $\endgroup$ – Bichoy May 26 '15 at 3:11
  • $\begingroup$ @novice Sorry, I had a syntax error in the compiled version (when copying) please recheck that. $\endgroup$ – Bichoy May 26 '15 at 3:14
  • $\begingroup$ @@Bichoy, I think Compile fails for functions containing D. The compiled version is slower. $\endgroup$ – novice May 26 '15 at 4:59
  • 2
    $\begingroup$ @novice, @bichoy, always worth checking this list of compilable functions before looking into Compile $\endgroup$ – dr.blochwave May 26 '15 at 6:45

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;
 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.

 glurg = Outer[
       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  *)
  • $\begingroup$ @@Michael E2, the idea of vectorization is great. I get 20x speedup. I am curious about where does the time-saving come from? What if normalVector is also a function of θ? $\endgroup$ – novice May 27 '15 at 1:14
  • $\begingroup$ @novice, for many built-in math functions, f[array] is handled in the math library (in low-level code) and is several times faster than, say, Map[f, array, {-1}]. See mathematica.stackexchange.com/a/21863 for more tips. As for normalVector, it would depend on how it is a function of θ. $\endgroup$ – Michael E2 May 27 '15 at 1:25

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[
     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: enter image description here


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.