(C) Copyright 1994 D.C.M. Burbulla and  C.T.J. Dodson
       Notebook to support the book
     Self-Tutor for Computer  Calculus

Using Mathematica 2

by
     D. C. M. Burbulla and C. T. J. Dodson, University of Toronto
 Published by:    
    Prentice-Hall Canada  (English)   ISBN 0130152803     and 
     Prentice-Hall Japan  (Japanese) ISBN4-8101-8558-3 C3055 P2880E

Properties of Differentiable Functions

Consider the following polynomial function and investigate critical points:

Clear[f]; f[x_] := 1 - 3 x + x^3

Solve[f'[x] == 0, x]  (* critical points *)

Solve[f''[x] == 0, x]  (* possible inflection points *)

Plot[{f[x], f'[x], f''[x]}, {x, -3, 3}];

ListPlot[{{-1, f[-1]}, {0, f[0]}, {1, f[1]}}];

Show[%, %%, PlotLabel -> {f[x], f'[x], f''[x]}];

Curvature

Now we introduce some new functions. Recall that  curvature is defined in terms of 
the angles that tangents make with respect to the  x-axis.

Clear[f, h]; Limit[(ArcTan[f'[x + h]] - ArcTan[f'[x]])/
        Sqrt[h^2 + (f[x + h] - f[x])^2], h -> 0, Analytic -> True]

Chap5.m

curvature/: curvature[f_][x_] :=
            Derivative[1][Derivative[1][f]][x]/(1 + Derivative[1][f][x]^2)^(3/2)
	    

radiusofcurvature/: radiusofcurvature[f_][x_] := Abs[1/curvature[f][x]]


viewcurvature/: viewcurvature[f_, {a_, b_}] := 
                Block[{x}, Plot[{f[x], curvature[f][x]}, {x, a, b},
		PlotLabel -> {f[x], "its graph and its curvature"}]]
		
viewcurvature::usage = "viewcurvature[f, {a, b}] plots f and its curvature function 
on the interval {a, b}."


centreofcurvature/: centreofcurvature[f_][s_] := 
                    {s - (Derivative[1][f][s]*(1 + Derivative[1][f][s]^2))/
		    Derivative[1][Derivative[1][f]][s],
		    f[s] + (1 + Derivative[1][f][s]^2)/
		    Derivative[1][Derivative[1][f]][s]}
		    
centreofcurvature::usage = "centreofcurvature[f][s] finds the co-ordinates of the centre 
of curvature of f at x = s."


viewevolute/: viewevolute[f_, {a_, b_}] := Block[{x, t},
              ParametricPlot[centreofcurvature[f][t], {t, a, b},
	      PlotLabel -> {"the evolute of", f[x]}]]
	      
viewevolute::usage = "viewevolute[f, {a, b}] plots the centres of curvature of f for all x 
in the interval {a, b}. This curve is called the evolute of f."


viewosculatingcircle/: viewosculatingcircle[f_][s_] := (one = Block[{t, x},
                       ParametricPlot[{radiusofcurvature[f][s]*Cos[t] +
		       centreofcurvature[f][s][[1]], radiusofcurvature[f][s]*
		       Sin[t] + centreofcurvature[f][s][[2]]}, {t, 0, 2*Pi}, 
		       AspectRatio -> Automatic, DisplayFunction -> Identity]];
		       two = ListPlot[{centreofcurvature[f][s], {s, f[s]}},
		       DisplayFunction -> Identity];
		       three = Block[{x}, Plot[f[x], {x, 
		      centreofcurvature[f][s][[1]] - 3*radiusofcurvature[f][s],
		     centreofcurvature[f][s][[2]] + 3*radiusofcurvature[f][s]},
		       DisplayFunction -> Identity]];
		       Block[{x}, Show[one, two, three, PlotLabel -> 
		       {f[x], "its osculating circle at x = ", s},
		       DisplayFunction -> $DisplayFunction]])
		       
viewosculatingcircle::usage = "viewosculatingcircle[f][s] plots the osculating circle of f 
tangent to f at x = s."


survey/: survey[f_, {a_, b_}, s_] := 
    (one = Block[{t}, ParametricPlot[
    {radiusofcurvature[f][s]*Cos[t] + centreofcurvature[f][s][[1]], 
     radiusofcurvature[f][s]*Sin[t] + centreofcurvature[f][s][[2]]},
       {t, 0, 2*Pi},  AspectRatio -> Automatic, 
     DisplayFunction -> Identity]];
    two = Block[{v, u}, If[N[f'[s]] == 0, 
    ParametricPlot[{s, f[s] + v},
       {v, -N[radiusofcurvature[f][s]], 
       N[radiusofcurvature[f][s]]},
                 DisplayFunction -> Identity],
 Plot[(s - u)/
f'[s] + f[s], {u, s - N[Abs[s - centreofcurvature[f][s][[1]]]],
                   s + N[Abs[s - centreofcurvature[f][s][[1]]]]},
                DisplayFunction -> Identity]]];
 three = Block[{m}, Plot[f'[s]*(m - s) + f[s], {m, s - 1, s + 1},
               DisplayFunction -> Identity]];
 four = Block[{x}, Plot[{f[x], curvature[f][x]}, {x, a, b},
                DisplayFunction -> Identity]];
 five = Block[{z}, ParametricPlot[centreofcurvature[f][z], {z, a, b},
                DisplayFunction -> Identity]];
 six = ListPlot[{{s, f[s]}, centreofcurvature[f][s]},
                DisplayFunction -> Identity];
 Block[{x}, Show[one, two, three, four, five, six,
  PlotLabel -> {f[x], "survey of curvature details at x =", s},
                DisplayFunction -> $DisplayFunction ]])

survey::usage = 
"survey[f, {a, b}, s] plots f, its curvature function 
and its evolute for x in the interval {a, b}. 
In addition, it includes the osculating circle tangent to f at 
x = s, plus short normal and tangent lines to f at x = s."
 
d/: d[p_, q_] := 4*p^3 + 27*q^2
 
 
cubicplot/: cubicplot[p_, q_] := Block[{x}, If[p < 0, 
Plot[x^3 + p*x + q, 
{x, N[-(3*Abs[Sqrt[-p]/Sqrt[3]])], N[3*Abs[Sqrt[-p]/Sqrt[3]]]},
     PlotLabel -> Which[d[p, q] > 0,
{x^3 + p x + q, "one real root", If[d[p, q] == 1, "1 = D",
                                                  d[p, q] " = D"]}, 
                       d[p, q] == 0,
{x^3 + p x + q, "two real roots", "0 = D"},
                       d[p, q] < 0,
{x^3 + p x + q, "three real roots", If[d[p, q] == -1, "-1 = D",
                                              d[p, q] " = D"]}]],
Plot[x^3 + p*x + q, {x, N[-3*q^(1/3)] - 1, N[3*q^(1/3)] + 1},
                PlotLabel ->
{x^3 + p x + q, "one real root", " p is non-negative"}]]]

cubicplot::usage = 
"cubicplot[p, q] plots the cubic x^3 + p*x + q 
on a suitably chosen interval. 
The PlotLabel indicates how many real 
solutions there are to the equation x^3 + p*x + q == 0, 
and also indicates the value
of the discriminant, D = 4*p^3 + 27*q^2."

?curvature

curvature[Sin][x]

curvature[Sin][{Pi/2, 0, -3 Pi/2}]    (* listable  in  x *)

?viewcurvature

viewcurvature[Sin, {-4, 4}];

Clear[f]; f[x_] := x^2

curvature[f][x]

curvature[f][0]  (* recall: f(x) = x^2, still *)

?radiusofcurvature

radiusofcurvature[f][0]

?centreofcurvature

centreofcurvature[f][0]

viewcurvature[f, {-3, 3}];

?viewevolute

viewevolute[f, {-3, 3}];  (* f(x)   is  still  x^2  *);

Plot[f[x], {x, -3, 3}];

Show[%%, %];

?viewosculatingcircle

viewosculatingcircle[Sin][-0.5];

viewosculatingcircle[Sin][0.5];

Show[%%, %,PlotLabel->None,PlotRange->{{-10,10},{-9,9}}];

Plot[Sin[x], {x, -4, 4}];

Plot[curvature[Sin][x], {x, -4, 4}];

viewevolute[Sin, {-Pi, Pi}];

Show[%%%, %%, %];

?survey

survey[f, {-3, 3}, 0];  (*  f(x)  is  still   x^2  *)

survey[f, {-3, 3}, 0.5];