Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
Thanks for the Accept. I take it my suggestions helped? Any follow-up questions? –  Mr.Wizard Aug 13 '13 at 14:45

2 Answers 2

up vote 9 down vote accepted

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}]
   ]
  ]
 ]
share|improve this answer
    
Thanks. Very educational. –  rogerl Aug 13 '13 at 14:49
    
@rogerl I just noticed something as I looked at this again; fn should be in a Module not Block, in case the symbol fn might appear in f; alternatively fn should be replaced with a Formal symbol. I need to revise my answer. –  Mr.Wizard Aug 13 '13 at 14:52
    
I would go even as far as setting var_ to var_Symbol to add an extra bit of protection, especially as your passing it into Block. –  rcollyer Aug 13 '13 at 15:22
    
@rcollyer That's not a bad idea, however curvature[foo, {1, 2, 3}] produces Block::lvsym: Local variable specification {1} contains 1, which is not a symbol or an assignment to a symbol. >> and halts, which IMHO is also not a bad behavior; the problem is safely made clear without having to add such messages to curvature itself. Do you disagree? –  Mr.Wizard Aug 13 '13 at 15:25
1  
@rcollyer Archimedes continues to laugh at me frequently. :o) –  Mr.Wizard Aug 13 '13 at 15:48

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}}$$

share|improve this answer
    
Ack. Thanks for catching this. –  rogerl Aug 13 '13 at 14:49

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.