Looking for feedback on this function

Question

I'm not sure this is an appropriate question; if not, I'm sure someone will close it :)

I've got the following function, listed below; its job is to accept a graph (function) presented either as a single-variable function (e.g., x^2), or as a parametrized function (e.g. {t Cos[t], t Sin[t]}) and compute its curvature at some set of values of the variable:

SetAttributes[curvature,HoldAll];
SyntaxInformation[curvature]={"LocalVariables"->{"Table",{2,2}},
    "ArgumentsPattern"->{_,_}};
curvature[f_, pts_] := Block[{extvar,var,fn,i,res},
    extvar = ReleaseHold[Hold[pts]/.{x_,y__}:>HoldPattern[x]];
    fn = ReleaseHold[Hold[f]/.extvar:>var];
    If[!ListQ[fn], fn={var,fn}];
    res = Table[1/Sqrt[Total[D[fn,var]^2]]/.var->i,
                Evaluate[Join[{i},Rest[pts]]]];
    If[Length[res]==1,res=res[[1]]];res]

The second parameter is in "Table" form, and is interpreted as for Table. The first two lines of the function extract the formal variable used in f and substitute occurrences of it in f by a local variable. The third line converts a single-variable function to parametric form. The next line actually does the work, computing the relevant formula over the set of values defined by the second parameter. Finally, if the result is a singleton, it gets converted back to a scalar by the final line of the function.

So for example

curvature[ x^2, {x,3} ]
{1/Sqrt[5], 1/Sqrt[17], 1/Sqrt[37]}

I'm looking for feedback on the way this function is written. For example (but not only), have I properly extracted the formal parameter (x in the example above) and the function? Is there a better way to evaluate the formula at the values specified by the second parameter? Finally, how about the manipulation in the last line? Any other feedback is welcome as well.

Answer 1

Shotgun thoughts:

• You don't need the Hold/ReleaseHold pair; Unevaluated will do: Unevaluated[f] /. rules

• You can use direct destructuring to extract extvar: curvature[f_, range : {var_, __}] :=

• By extracting var as above you can Block it directly, simplifying everything.

• You can leave the Table variable out of the main Block as it is already localized.

• res /. {z_} :> z can replace the line If[Length[res] == 1, res = res[[1]]]; res

• In fact res can be eliminated and the rule applied to the Table

• You can pre-evaluate your 1/Sqrt[Total[D[fn,var]^2]] expression; this should be faster, and also eliminates the need for the replacement inside the Table

Putting it together I would write:

SetAttributes[curvature, HoldAll];

curvature[f_, range : {var_, __}] :=
  Block[{var},
    Module[{fn},
      fn = If[ListQ[f], f, {var, f}];
      fn = 1/Sqrt[Total[D[fn, var]^2]];
      Table[fn, range] /. {z_} :> z
    ]
  ]

You will notice red highlighting of var signifying a possible conflict; this is not a problem and is in fact exactly what I desire: the sharing of var between the function, the Block, and the Table.

---

My answer originally had fn in the Block declaration. This is bad because conceivably fn could appear in the function expression and this would conflict. I have moved it to a Module now.

An alternative is to do without that symbol entirely by writing the lines of code as Functions:

curvature[f_, range : {var_, __}] :=
 Block[{var},
  Table[#, range] /. {z_} :> z &[
   1/Sqrt[Total[D[#, var]^2]] &[
    If[ListQ[f], f, {var, f}]
   ]
  ]
 ]

Answer 2

The main problem, as I see it: your formula for the curvature is incorrect. Here is a derivation of the curvature using Mathematica:

Clear[x, y, t, vec, velocity, tangent, speed, curvature]

tangent[t_] = Simplify[velocity[t]/speed[t]];

curvature[t_] = 
  Sqrt[#.#] &@FullSimplify[D[tangent[t], t]/speed[t], speed[t] > 0];

vec[t_] = Through[{x, y}[t]]

(* ==> {x[t], y[t]} *)

velocity[t_] = D[vec[t], t];

speed[t_] = Sqrt[velocity[t].velocity[t]];

Simplify[curvature[t]]

$$\sqrt{\frac{\left(x''(t) y'(t)-x'(t) y''(t)\right)^2}{\left(x'(t)^2+y'(t)^2\right)^3}}$$