Classify 3D by Jared Antrobus

Modifications by Andy Long, Fall, 2015, to add contour plot and streamlines to the surface.

Classify3D

classify3D[expr_,var1_,var2_]:=Module[
{func,solution,critical,discriminant,fxx,dValues,fxxValues,extrema,plotLimits,classification,f,f1,g1,sP,vP},
(* Define the function *)
func=Function[Evaluate[{var1,var2}],Evaluate[expr]];
(* Solve for a solution set *)
solution={var1,var2,func[var1,var2]}/.Solve[{
D[expr,var1]==0 ,
D[expr,var2]==0
},{var1,var2},Reals];
(* Check to see if Solve was successful *)
If[Apply[Or,Map[Head[#]===Symbol&,Flatten[solution]]],Column[{"Error: Could not solve for a finite set of critical points."}],
(* Define two dimensional critical points *)
critical=Take[solution,All,{1,2}];
(* Construct the discriminant and (Subscript[f, xx]) *)
discriminant=Function[Evaluate[{var1,var2}],Evaluate[D[expr,var1,var1]*D[expr,var2,var2]-((D[expr,var1,var2])^2)]];
fxx=Function[Evaluate[{var1,var2}],Evaluate[D[expr,var1,var1]]];
f=Function[Evaluate[{var1,var2}],Evaluate[expr]];
g1=Function[Evaluate[{var1,var2}],Evaluate[D[expr,{{var1,var2}}]]];
(* Find discriminant and fxx values *)
dValues=Apply[discriminant,critical,1];
fxxValues=Apply[fxx,critical,1];
(* Determine plot limits to use *)
extrema={
{Min[Take[critical,All,{1}]],Max[Take[critical,All,{1}]]},
{Min[Take[critical,All,{2}]],Max[Take[critical,All,{2}]]}
};
plotLimits=extrema+{{-1,1},{-1,1}}+{
{extrema[[1,1]]-extrema[[1,2]],extrema[[1,2]]-extrema[[1,1]]},
{extrema[[2,1]]-extrema[[2,2]],extrema[[2,2]]-extrema[[2,1]]}
}/10;
(* Classify each critical point *)
classification=Table[
If[dValues[[i]]<0,"saddle",
If[dValues[[i]]==0,"inconclusive",
If[dValues[[i]]>0,
If[fxxValues[[i]]>0,"minimum",
If[fxxValues[[i]]<0,"maximum",
"inconclusive"]]]]],
{i,Length[critical]}];
(* Output to screen *)
vP=VectorPlot[g1[var1,var2],{var1,plotLimits[[1,1]],plotLimits[[1,2]]},{var2,plotLimits[[2,1]],plotLimits[[2,2]]},VectorStyle->{Arrowheads[0.005]},VectorColorFunction->Hue,ImageSize->400];
sP=StreamPlot[g1[var1,var2],{var1,plotLimits[[1,1]],plotLimits[[1,2]]},{var2,plotLimits[[2,1]],plotLimits[[2,2]]},StreamStyle->"Line",StreamPoints->80,StreamColorFunction->Function[{x,y},Darker@Yellow]];
Column[{
TableForm[
{critical,
classification,
solution[[All,3]],
dValues,
fxxValues},
TableDepth->2,TableAlignments->Center,TableHeadings->{{"Critical Point","Classification","Function Value","D Value","(Subscript[f, xx]) Value"},None}
],,
Row[{"Discriminant: ",discriminant[var1,var2]}],
Row[{"(Subscript[f, xx]): ",fxx[var1,var2]}],
Show[
Plot3D[expr,{var1,plotLimits[[1,1]],plotLimits[[1,2]]},{var2,plotLimits[[2,1]],plotLimits[[2,2]]},ImageSize->{400,Automatic},MeshFunctions->{#3&},Mesh->10,ColorFunction->ColorData["StarryNightColors"]],
Graphics3D[vP[[1]]/.{{x_Real,y_Real}:>{x,y,f[x,y]}}],
Graphics3D[sP[[1]]/.{{x_Real,y_Real}:>{x,y,f[x,y]}}],
Table[
If[classification[[i]]==="maximum",ListPointPlot3D[Tooltip[{solution[[i]]}],PlotStyle->Directive[PointSize[.04],Red]],
If[classification[[i]]==="minimum",ListPointPlot3D[Tooltip[{solution[[i]]}],PlotStyle->Directive[PointSize[.04],Blue]],
If[classification[[i]]==="saddle",ListPointPlot3D[Tooltip[{solution[[i]]}],PlotStyle->Directive[PointSize[.04],Purple]],
ListPointPlot3D[Tooltip[{solution[[i]]}],PlotStyle->Directive[PointSize[.04],Orange]]
]]]
,{i,Length[classification]}]
],
(* debug *)

(* /debug *)
}]
]
]

ael: Thanks to
http://mathematica.stackexchange.com/questions/18471/adding-contour-lines-to-a-3d-plot
I learned the best way to add the contour lines to a Plot3D; 
from 
https://www.deltacollege.edu/emp/nbykov/VectorCalculus.html
I got the tips for the vectors.

Plot3D[Sin[x]+Sin[y],{x,-5,5},{y,-5,5},MeshFunctions->{#3&},Mesh->7,ColorFunction->ColorData["StarryNightColors"]]

Directions

The classify3D module is of the form classify3D[expression,variable1,variable2].
It is used for classifying the critical points of fairly nice two-variable functions.
The expression should represent a two-variable function.  (i.e. 2(x^2)+(y^4))
The variables should be the ones used in the function.  (i.e. x and y)
The module will fail if the number of critical points is infinite or zero.
NOTE:  classify3D has issues handling functions that have undefined segments.  (i.e. (x^2)+((y^3)/x))

Examples: 
we start with some known examples, and then try some 		stranger cases....

classify3D[(x^2)-(y^2),x,y]

classify3D[(x^2)+(y^2),x,y]

classify3D[(x^3)+(y^4)-6x-2(y^2),x,y]

classify3D[((x^2)-7x)((y^2)-2y),x,y]

classify3D[(x^3)+(y^2),x,y]

classify3D[(x^3)+y,x,y]

Example 6, p. 1084:

classify3D[(x^2)y-(y^3),x,y]

h[x_,y_]=1/Sqrt[x^2+y^2];
f1[x_,y_]=h[x-7,y-10]+h[x-21,y-10];
classify3D[f1[x,y],x,y]