(* Content-type: application/vnd.wolfram.mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 10.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       158,          7]
NotebookDataLength[     31675,        913]
NotebookOptionsPosition[     29817,        856]
NotebookOutlinePosition[     30170,        872]
CellTagsIndexPosition[     30127,        869]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{

Cell[CellGroupData[{
Cell["Newton\[CloseCurlyQuote]s method", "Title",
 CellChangeTimes->{{3.626602737536145*^9, 3.626602753959387*^9}}],

Cell["\<\
A system of m equations in n unknowns can always be thought of as a single \
equation of the form
\tF(x) = b
where F:R^n -> R^m is a function (with components F_i(x) equal to the left \
side of the ith equation), x is the unknown vector in R^n that we wish to \
solve for, and b is a vector in R^m with components b_i equal to the number \
of the right side of the ith equation.

For instance, if b1, b2 are given real numbers and we want to find the \
intersections between the vertical parabola y=x^2- b1 and the horizontal \
parabola x=y^2 - b2, then we must find the solutions (x,y) of the 2 x 2 \
system
\ty - x^2 = b1
\tx - y^2 = b2
which may be written F(x,y) = b where b=(b1,b2) and\
\>", "Text",
 CellChangeTimes->{{3.6266027636078587`*^9, 3.6266032400841627`*^9}, {
  3.6266033527220984`*^9, 3.6266033689153357`*^9}, {3.626603446625832*^9, 
  3.62660386500418*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{"F", "[", 
   RowBox[{"{", 
    RowBox[{"x_", ",", "y_"}], "}"}], "]"}], " ", "=", " ", 
  RowBox[{"{", 
   RowBox[{
    RowBox[{"y", "-", 
     RowBox[{"x", "^", "2"}]}], ",", 
    RowBox[{"x", "-", 
     RowBox[{"y", "^", "2"}]}]}], "}"}]}]], "Input",
 CellChangeTimes->{{3.626603867607401*^9, 3.6266038893277063`*^9}}],

Cell["\<\

Try to solve this system exactly, and you'll likely end up trying to identify \
the roots of a degree four polynomial.  There's a known way to do this, but \
it's not pretty.   Take b = (-2,-3) for example.  The following picture \
clearly shows four solutions to the system:\
\>", "Text",
 CellChangeTimes->{{3.6266027636078587`*^9, 3.6266032400841627`*^9}, {
  3.6266033527220984`*^9, 3.6266033689153357`*^9}, {3.626603446625832*^9, 
  3.62660386500418*^9}, {3.626604033195619*^9, 3.6266040721238728`*^9}}],

Cell[BoxData[
 RowBox[{"ContourPlot", "[", 
  RowBox[{
   RowBox[{"{", 
    RowBox[{
     RowBox[{
      RowBox[{"y", "-", 
       RowBox[{"x", "^", "2"}]}], "\[Equal]", 
      RowBox[{"-", "2"}]}], ",", 
     RowBox[{
      RowBox[{"x", "-", 
       RowBox[{"y", "^", "2"}]}], "\[Equal]", 
      RowBox[{"-", "3"}]}]}], "}"}], ",", 
   RowBox[{"{", 
    RowBox[{"x", ",", 
     RowBox[{"-", "5"}], ",", "5"}], "}"}], ",", 
   RowBox[{"{", 
    RowBox[{"y", ",", 
     RowBox[{"-", "5"}], ",", "5"}], "}"}]}], "]"}]], "Input",
 CellChangeTimes->{{3.626603972619898*^9, 3.626604017775572*^9}}],

Cell[TextData[{
 "And ",
 StyleBox["Mathematica",
  FontSlant->"Italic"],
 " will find them if you ask it to."
}], "Text",
 CellChangeTimes->{{3.626604090435196*^9, 3.62660411331464*^9}}],

Cell[BoxData[
 RowBox[{"TableForm", "[", 
  RowBox[{
   RowBox[{"{", 
    RowBox[{"x", ",", "y"}], "}"}], "/.", 
   RowBox[{"Simplify", "[", 
    RowBox[{"Solve", "[", 
     RowBox[{
      RowBox[{
       RowBox[{"F", "[", 
        RowBox[{"{", 
         RowBox[{"x", ",", "y"}], "}"}], "]"}], "\[Equal]", 
       RowBox[{"{", 
        RowBox[{
         RowBox[{"-", "2"}], ",", 
         RowBox[{"-", "3"}]}], "}"}]}], ",", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y"}], "}"}]}], "]"}], "]"}]}], "]"}]], "Input",
 CellChangeTimes->{{3.626603833442603*^9, 3.626603838110785*^9}, {
  3.626603896862872*^9, 3.626603957767727*^9}, {3.626604146144355*^9, 
  3.6266041807581778`*^9}}],

Cell["\<\

And this is about as easy as non-linear systems of equations get.  Typically, \
it's impossible to find exact solutions to non-linear systems. 

Nevertheless, when the function F is differentiable and there is a reasonably \
good initial guess x=x0 at a solution of F(x) = b, one can often use linear \
approximation to improve the guess repeatedly, thereby obtaining arbitrarily \
good approximations of an exact solution.  This is known as Newton's method.  \
 Specificly, one replaces the function F in F(x) = b with its linear \
approximation at the initial guess x0 to get a new equation

F(x0) + Df(x0) (x-x0) = b

The great thing is that this is a linear system (as you can see by moving the \
F(x0) term to the other side of the equation), and we have a fairly \
straightforward way to solve linear systems.  To illustrate, let\
\[CloseCurlyQuote]s take the equation F(x,y) = (-2,-3) above and observe that \
if x0 = (2,2), then F(x0) = (-3,-3) is fairly close to our target vector.  We \
then employ Newton\[CloseCurlyQuote]s method as follows:\
\>", "Text",
 CellChangeTimes->{{3.6266027636078587`*^9, 3.6266032400841627`*^9}, {
  3.6266033527220984`*^9, 3.6266033689153357`*^9}, {3.626603446625832*^9, 
  3.626603827148531*^9}, {3.626604188962425*^9, 3.626604260906702*^9}, {
  3.626604368817635*^9, 3.626604467129157*^9}}],

Cell[BoxData[
 RowBox[{"MatrixForm", "[", 
  RowBox[{
   RowBox[{"DF", "[", 
    RowBox[{"{", 
     RowBox[{"x_", ",", "y_"}], "}"}], "]"}], " ", "=", " ", 
   RowBox[{"D", "[", 
    RowBox[{
     RowBox[{"F", "[", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y"}], "}"}], "]"}], ",", 
     RowBox[{"{", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y"}], "}"}], "}"}]}], "]"}]}], "]"}]], "Input",
 CellChangeTimes->{{3.6266042783173523`*^9, 3.626604305548889*^9}, {
  3.626604557963168*^9, 3.626604565402893*^9}}],

Cell[BoxData[{
 RowBox[{"x0", " ", "=", " ", 
  RowBox[{"{", 
   RowBox[{"2.", ",", "2"}], "}"}]}], "\[IndentingNewLine]", 
 RowBox[{"b", "=", 
  RowBox[{"{", 
   RowBox[{
    RowBox[{"-", "2"}], ",", 
    RowBox[{"-", "3"}]}], "}"}]}]}], "Input",
 CellChangeTimes->{{3.626604472598579*^9, 3.626604477018469*^9}, {
   3.626604522225604*^9, 3.6266045427450857`*^9}, {3.626604883727392*^9, 
   3.6266048931663723`*^9}, 3.6266311381595497`*^9}],

Cell[BoxData[
 RowBox[{
  RowBox[{"Linapprox", "[", 
   RowBox[{"{", 
    RowBox[{"x_", ",", "y_"}], "}"}], "]"}], "=", " ", 
  RowBox[{"Simplify", "[", 
   RowBox[{
    RowBox[{"F", "[", "x0", "]"}], "+", 
    RowBox[{
     RowBox[{"DF", "[", "x0", "]"}], ".", 
     RowBox[{"(", 
      RowBox[{
       RowBox[{"{", 
        RowBox[{"x", ",", "y"}], "}"}], "-", "x0"}], ")"}]}]}], 
   "]"}]}]], "Input",
 CellChangeTimes->{{3.626604578426538*^9, 3.626604644242029*^9}, {
  3.626604823857235*^9, 3.6266048264806833`*^9}}],

Cell["Our new guess is", "Text",
 CellChangeTimes->{{3.626604707861247*^9, 3.626604711573153*^9}}],

Cell[BoxData[
 RowBox[{"x1", "=", 
  RowBox[{
   RowBox[{"{", 
    RowBox[{"x", ",", "y"}], "}"}], "/.", 
   RowBox[{"First", "[", 
    RowBox[{"Solve", "[", 
     RowBox[{
      RowBox[{
       RowBox[{"Linapprox", "[", 
        RowBox[{"{", 
         RowBox[{"x", ",", "y"}], "}"}], "]"}], "\[Equal]", "b"}], ",", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y"}], "}"}]}], "]"}], "]"}]}]}]], "Input",
 CellChangeTimes->{{3.6266044812866573`*^9, 3.62660451809032*^9}, {
  3.626604654641302*^9, 3.626604689665265*^9}, {3.626604727313478*^9, 
  3.626604763304555*^9}, {3.626604878264426*^9, 3.626604907007791*^9}}],

Cell["Note that F(x1) is much closer to (-2,-3):", "Text",
 CellChangeTimes->{{3.626604776044773*^9, 3.6266047943322573`*^9}}],

Cell[BoxData[
 RowBox[{"N", "[", 
  RowBox[{
   RowBox[{"F", "[", "x1", "]"}], ",", "25"}], "]"}]], "Input",
 CellChangeTimes->{{3.626604797782711*^9, 3.62660480267865*^9}, {
  3.6266050055338917`*^9, 3.626605014653029*^9}}],

Cell["\<\
So already we have a lot of improvement.  However, the really great thing \
about Newton\[CloseCurlyQuote]s method is that you can repeat this step to \
improve your guess further:\
\>", "Text",
 CellChangeTimes->{{3.626604838731378*^9, 3.626604874851967*^9}}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"Linapprox", "[", 
   RowBox[{"{", 
    RowBox[{"x_", ",", "y_"}], "}"}], "]"}], "=", " ", 
  RowBox[{"Simplify", "[", 
   RowBox[{
    RowBox[{"F", "[", "x1", "]"}], "+", 
    RowBox[{
     RowBox[{"DF", "[", "x1", "]"}], ".", 
     RowBox[{"(", 
      RowBox[{
       RowBox[{"{", 
        RowBox[{"x", ",", "y"}], "}"}], "-", "x1"}], ")"}]}]}], 
   "]"}]}], "\[IndentingNewLine]", 
 RowBox[{"x1", "=", 
  RowBox[{
   RowBox[{"{", 
    RowBox[{"x", ",", "y"}], "}"}], "/.", 
   RowBox[{"First", "[", 
    RowBox[{"Solve", "[", 
     RowBox[{
      RowBox[{
       RowBox[{"Linapprox", "[", 
        RowBox[{"{", 
         RowBox[{"x", ",", "y"}], "}"}], "]"}], "\[Equal]", "b"}], ",", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y"}], "}"}]}], "]"}], 
    "]"}]}]}], "\[IndentingNewLine]", 
 RowBox[{"N", "[", 
  RowBox[{
   RowBox[{"F", "[", "x1", "]"}], ",", "25"}], "]"}]}], "Input",
 CellChangeTimes->{{3.626604469132805*^9, 3.62660447057718*^9}, {
  3.626604931445557*^9, 3.626604966759165*^9}, {3.62660502266361*^9, 
  3.6266050245508137`*^9}, {3.6266051189262114`*^9, 3.626605119389862*^9}}],

Cell["...and so on as often as you like", "Text",
 CellChangeTimes->{{3.626604995330578*^9, 3.6266050012346163`*^9}}],

Cell[CellGroupData[{

Cell["Nonlinear systems with more equations than unknowns", "Subsection",
 CellChangeTimes->{{3.6267756109421997`*^9, 3.626775625487554*^9}}],

Cell["\<\
Both in theory and in practice, one often needs to solve (or at least \
understand solutions of) non-linear systems with more equations than \
unknowns.  Consider for instance the single equation f[x,y] = 0 where \
\>", "Text",
 CellChangeTimes->{{3.626775643145893*^9, 3.626775780994412*^9}, {
  3.6267764537049303`*^9, 3.626776453856785*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{"f", "[", 
   RowBox[{"{", 
    RowBox[{"x_", ",", "y_"}], "}"}], "]"}], "=", 
  RowBox[{
   RowBox[{"x", "^", "3"}], "+", 
   RowBox[{"y", "^", "3"}], "-", 
   RowBox[{"2", "x", " ", "y"}]}]}]], "Input",
 CellChangeTimes->{{3.6267764637093067`*^9, 3.626776482507841*^9}, {
  3.626776520083549*^9, 3.626776526706479*^9}}],

Cell["\<\
There are many points (x,y) that solve f(x,y)=0.  In fact the set of all such \
points, which includes (1,1) is a classical curve known as a `folium of \
Descartes\[CloseCurlyQuote], which looks like this\
\>", "Text",
 CellChangeTimes->{{3.626775883539014*^9, 3.626775924951612*^9}, {
  3.626776625634391*^9, 3.626776695246066*^9}}],

Cell[BoxData[
 RowBox[{"Show", "[", 
  RowBox[{
   RowBox[{"ContourPlot", "[", 
    RowBox[{
     RowBox[{
      RowBox[{"f", "[", 
       RowBox[{"{", 
        RowBox[{"x", ",", "y"}], "}"}], "]"}], "\[Equal]", "0"}], ",", 
     RowBox[{"{", 
      RowBox[{"x", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"y", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}]}], "]"}], ",", 
   RowBox[{"Graphics", "[", 
    RowBox[{"{", 
     RowBox[{
      RowBox[{"PointSize", "\[Rule]", "Large"}], ",", "Red", ",", 
      RowBox[{"Point", "[", 
       RowBox[{"{", 
        RowBox[{"1", ",", "1"}], "}"}], "]"}]}], "}"}], "]"}]}], 
  "]"}]], "Input",
 CellChangeTimes->{{3.626776559670504*^9, 3.626776579005047*^9}, {
  3.626776699859497*^9, 3.62677674331491*^9}, {3.626776774873802*^9, 
  3.6267767911763067`*^9}}],

Cell[TextData[{
 "The red point is (1,1).  From the picture it\[CloseCurlyQuote]s evident \
that if we vary either the x or the y coordinate of (1,1) a ",
 StyleBox["little bit ",
  FontWeight->"Bold",
  FontSlant->"Italic"],
 "then we can change the other coordinate a little bit to compensate and \
obtain another point on the curve.  For instance, we can hope that there is a \
point x near 1 such that (x,1..05) is still on the curve.  We can again use \
Newton\[CloseCurlyQuote]s method to help find this point.  But since we \
already have a value for y in mind, we might as well think of f as a function \
h(x) = f(x,1.05) as a function of x only"
}], "Text",
 CellChangeTimes->{{3.6267768325498047`*^9, 3.626776891148275*^9}, {
  3.626776931192301*^9, 3.6267770479261217`*^9}, {3.626777092414194*^9, 
  3.626777164475293*^9}}],

Cell[BoxData[{
 RowBox[{"y0", "=", ".8"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"h", "[", "x_", "]"}], " ", "=", 
  RowBox[{"f", "[", 
   RowBox[{"{", 
    RowBox[{"x", ",", "y0"}], "}"}], "]"}]}]}], "Input",
 CellChangeTimes->{{3.626777169562212*^9, 3.626777182009955*^9}, {
  3.62677741197337*^9, 3.62677744768827*^9}, {3.626777479539241*^9, 
  3.6267774892451982`*^9}, {3.627038714703485*^9, 3.627038740910533*^9}, {
  3.627038793493779*^9, 3.627038822861706*^9}}],

Cell[BoxData["`"], "Input",
 CellChangeTimes->{{3.627038788994595*^9, 3.627038789508079*^9}}],

Cell["\<\
We want to solve h[x]=0 and again resort to linear approximation:\
\>", "Text",
 CellChangeTimes->{{3.626777206298498*^9, 3.6267772407870407`*^9}}],

Cell[BoxData[{
 RowBox[{"step", " ", "=", " ", "0"}], "\[IndentingNewLine]", 
 RowBox[{"x0", " ", "=", "1"}], "\[IndentingNewLine]", 
 RowBox[{"Show", "[", 
  RowBox[{
   RowBox[{"ContourPlot", "[", 
    RowBox[{
     RowBox[{
      RowBox[{"f", "[", 
       RowBox[{"{", 
        RowBox[{"x", ",", "y"}], "}"}], "]"}], "\[Equal]", "0"}], ",", 
     RowBox[{"{", 
      RowBox[{"x", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"y", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}]}], "]"}], ",", 
   RowBox[{"Graphics", "[", 
    RowBox[{"{", 
     RowBox[{
      RowBox[{"PointSize", "\[Rule]", "Large"}], ",", "Red", ",", 
      RowBox[{"Point", "[", 
       RowBox[{"{", 
        RowBox[{"1", ",", "1"}], "}"}], "]"}]}], "}"}], "]"}], ",", 
   "\[IndentingNewLine]", 
   RowBox[{"Graphics", "[", 
    RowBox[{"{", 
     RowBox[{
      RowBox[{"PointSize", "\[Rule]", "Large"}], ",", "Blue", ",", 
      RowBox[{"Point", "[", 
       RowBox[{"{", 
        RowBox[{"x0", ",", "y0"}], "}"}], "]"}]}], "}"}], "]"}]}], 
  "]"}]}], "Input",
 CellChangeTimes->{{3.626777242439152*^9, 3.6267773312589273`*^9}, {
  3.6267775275710583`*^9, 3.6267775306753674`*^9}}],

Cell[BoxData[{
 RowBox[{"step", "+=", "1"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"lapprox", "[", "x_", "]"}], " ", "=", " ", 
  RowBox[{
   RowBox[{"h", "[", "x0", "]"}], "+", 
   RowBox[{
    RowBox[{
     RowBox[{"h", "'"}], "[", "x0", "]"}], 
    RowBox[{"(", 
     RowBox[{"x", "-", "x0"}], ")"}]}]}]}], "\[IndentingNewLine]", 
 RowBox[{"x0", "=", 
  RowBox[{"x", "/.", 
   RowBox[{"First", "[", 
    RowBox[{"Solve", "[", 
     RowBox[{
      RowBox[{
       RowBox[{"lapprox", "[", "x", "]"}], "\[Equal]", "0"}], ",", "x"}], 
     "]"}], "]"}]}]}], "\[IndentingNewLine]", 
 RowBox[{"h", "[", "x0", "]"}], "\[IndentingNewLine]", 
 RowBox[{"Show", "[", 
  RowBox[{
   RowBox[{"ContourPlot", "[", 
    RowBox[{
     RowBox[{
      RowBox[{"f", "[", 
       RowBox[{"{", 
        RowBox[{"x", ",", "y"}], "}"}], "]"}], "\[Equal]", "0"}], ",", 
     RowBox[{"{", 
      RowBox[{"x", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"y", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}]}], "]"}], ",", 
   RowBox[{"Graphics", "[", 
    RowBox[{"{", 
     RowBox[{
      RowBox[{"PointSize", "\[Rule]", "Large"}], ",", "Red", ",", 
      RowBox[{"Point", "[", 
       RowBox[{"{", 
        RowBox[{"1", ",", "1"}], "}"}], "]"}]}], "}"}], "]"}], ",", 
   "\[IndentingNewLine]", 
   RowBox[{"Graphics", "[", 
    RowBox[{"{", 
     RowBox[{
      RowBox[{"PointSize", "\[Rule]", "Large"}], ",", "Blue", ",", 
      RowBox[{"Point", "[", 
       RowBox[{"{", 
        RowBox[{"x0", ",", "y0"}], "}"}], "]"}]}], "}"}], "]"}]}], 
  "]"}]}], "Input",
 CellChangeTimes->{{3.626777242439152*^9, 3.626777336069889*^9}, {
  3.626777533277302*^9, 3.626777536289687*^9}}],

Cell["\<\
The implicit function theorem says (in the present circumstance) that we can \
solve the equation f[x,y]=0 for x in terms of y near any point at with the \
partial derivative of f with respect to x is non-zero.  So let\
\[CloseCurlyQuote]s find the bad points:\
\>", "Text",
 CellChangeTimes->{{3.626969955192374*^9, 3.626970043995327*^9}}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"fx", "[", 
   RowBox[{"{", 
    RowBox[{"x_", ",", "Y_"}], "}"}], "]"}], "=", 
  RowBox[{"D", "[", 
   RowBox[{
    RowBox[{"f", "[", 
     RowBox[{"{", 
      RowBox[{"x", ",", "y"}], "}"}], "]"}], ",", "x"}], 
   "]"}]}], "\[IndentingNewLine]", 
 RowBox[{"badpts", "=", 
  RowBox[{
   RowBox[{"{", 
    RowBox[{"x", ",", "y"}], "}"}], "/.", 
   RowBox[{"Solve", "[", 
    RowBox[{
     RowBox[{"{", 
      RowBox[{
       RowBox[{
        RowBox[{"f", "[", 
         RowBox[{"{", 
          RowBox[{"x", ",", "y"}], "}"}], "]"}], "\[Equal]", "0"}], ",", 
       RowBox[{
        RowBox[{"fx", "[", 
         RowBox[{"{", 
          RowBox[{"x", ",", "y"}], "}"}], "]"}], "\[Equal]", "0"}]}], "}"}], 
     ",", 
     RowBox[{"{", 
      RowBox[{"x", ",", "y"}], "}"}]}], "]"}]}]}], "\[IndentingNewLine]", 
 RowBox[{"Show", "[", 
  RowBox[{
   RowBox[{"ContourPlot", "[", 
    RowBox[{
     RowBox[{
      RowBox[{"f", "[", 
       RowBox[{"{", 
        RowBox[{"x", ",", "y"}], "}"}], "]"}], "\[Equal]", "0"}], ",", 
     RowBox[{"{", 
      RowBox[{"x", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"y", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}]}], "]"}], ",", 
   RowBox[{"ListPlot", "[", 
    RowBox[{"badpts", ",", 
     RowBox[{"PlotStyle", "\[Rule]", 
      RowBox[{"{", 
       RowBox[{"Red", ",", 
        RowBox[{"PointSize", "[", "Large", "]"}]}], "}"}]}]}], "]"}]}], 
  "]"}]}], "Input",
 CellChangeTimes->{{3.6269700453420153`*^9, 3.6269702047162647`*^9}, {
  3.626970278047146*^9, 3.626970290248948*^9}, {3.626977945372415*^9, 
  3.626977961443203*^9}}],

Cell["\<\
The general idea is that if f[x1,...,xn] is a decent (i.e. C^1) scalar \
function of n variables and f[a] = b, then one can solve the equation f[x] = \
b for xj in terms of the other components of x near a provided the jth \
partial derivative of f at a is non-zero.  The 2 dimensional sphere in R^3 \
provides a nice illustration\
\>", "Text",
 CellChangeTimes->{{3.626777827176696*^9, 3.626778007297243*^9}, {
  3.626969839688383*^9, 3.6269698511826878`*^9}}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"f", "[", 
   RowBox[{"{", 
    RowBox[{"x_", ",", "y_", ",", "z_"}], "}"}], "]"}], "=", 
  RowBox[{
   RowBox[{"x", "^", "2"}], "+", 
   RowBox[{"y", "^", "2"}], "+", 
   RowBox[{"z", "^", "2"}]}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"grad", "[", 
   RowBox[{"{", 
    RowBox[{"x_", ",", "y_", ",", "z_"}], "}"}], "]"}], "=", 
  RowBox[{"D", "[", 
   RowBox[{
    RowBox[{"f", "[", 
     RowBox[{"{", 
      RowBox[{"x", ",", "y", ",", "z"}], "}"}], "]"}], ",", 
    RowBox[{"{", 
     RowBox[{"{", 
      RowBox[{"x", ",", "y", ",", "z"}], "}"}], "}"}]}], 
   "]"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"plot1", "=", 
   RowBox[{"ContourPlot3D", "[", 
    RowBox[{
     RowBox[{
      RowBox[{"f", "[", 
       RowBox[{"{", 
        RowBox[{"x", ",", "y", ",", "z"}], "}"}], "]"}], "\[Equal]", "1"}], 
     ",", 
     RowBox[{"{", 
      RowBox[{"x", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"y", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"z", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", 
     RowBox[{"AxesLabel", "\[Rule]", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y", ",", "z"}], "}"}]}]}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"plot2", "=", 
   RowBox[{"ContourPlot3D", "[", 
    RowBox[{
     RowBox[{
      RowBox[{
       RowBox[{"grad", "[", 
        RowBox[{"{", 
         RowBox[{"x", ",", "y", ",", "z"}], "}"}], "]"}], "[", 
       RowBox[{"[", "2", "]"}], "]"}], "\[Equal]", "0"}], ",", 
     RowBox[{"{", 
      RowBox[{"x", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"y", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"z", ",", 
       RowBox[{"-", "2"}], ",", "2"}], "}"}], ",", 
     RowBox[{"AxesLabel", "\[Rule]", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y", ",", "z"}], "}"}]}]}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{"Show", "[", "plot1", "]"}]}], "Input",
 CellChangeTimes->{{3.6269698566435413`*^9, 3.6269699002190123`*^9}, {
  3.626977978231236*^9, 3.626978086785389*^9}, {3.626978448221264*^9, 
  3.6269785503119698`*^9}, {3.6269787169040728`*^9, 3.626978720759378*^9}, {
  3.6269788298462954`*^9, 3.626978843495346*^9}, {3.6269788979233027`*^9, 
  3.626978907931682*^9}, {3.626980052281767*^9, 3.626980107960671*^9}}],

Cell["The torus in R^3 gives a more elaborate example.", "Text",
 CellChangeTimes->{{3.626777827176696*^9, 3.626778007297243*^9}, 
   3.626969839688383*^9, {3.626978924952669*^9, 3.626978941404779*^9}}],

Cell[BoxData[{
 RowBox[{"R", "=", "2"}], "\[IndentingNewLine]", 
 RowBox[{"r", "=", "1"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"f", "[", 
    RowBox[{"{", 
     RowBox[{"x_", ",", "y_", ",", "z_"}], "}"}], "]"}], "=", 
   RowBox[{
    RowBox[{
     RowBox[{"(", 
      RowBox[{
       RowBox[{"x", "^", "2"}], "+", 
       RowBox[{"y", "^", "2"}], "+", 
       RowBox[{"z", "^", "2"}], "+", 
       RowBox[{"R", "^", "2"}], "-", 
       RowBox[{"r", "^", "2"}]}], ")"}], "^", "2"}], "-", 
    RowBox[{"4", 
     RowBox[{"R", "^", "2"}], " ", 
     RowBox[{"(", 
      RowBox[{
       RowBox[{"x", "^", "2"}], "+", 
       RowBox[{"y", "^", "2"}]}], ")"}]}]}]}], " "}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"grad", "[", 
   RowBox[{"{", 
    RowBox[{"x_", ",", "y_", ",", "z_"}], "}"}], "]"}], "=", 
  RowBox[{"Simplify", "[", 
   RowBox[{"D", "[", 
    RowBox[{
     RowBox[{"f", "[", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y", ",", "z"}], "}"}], "]"}], ",", 
     RowBox[{"{", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y", ",", "z"}], "}"}], "}"}]}], "]"}], 
   "]"}]}]}], "Input",
 CellChangeTimes->{{3.626775785854348*^9, 3.626775854221464*^9}, {
   3.626775943713763*^9, 3.626775945890586*^9}, {3.6267760966413717`*^9, 
   3.626776116752543*^9}, {3.626776153955818*^9, 3.626776162157522*^9}, {
   3.626776195275976*^9, 3.626776203732314*^9}, {3.626776266146792*^9, 
   3.6267762663918333`*^9}, {3.626776319041646*^9, 3.626776328172225*^9}, {
   3.6267763795979757`*^9, 3.6267764243794107`*^9}, {3.626778012455068*^9, 
   3.62677804785802*^9}, {3.6267781736087847`*^9, 3.626778176800837*^9}, {
   3.6267786698303328`*^9, 3.6267787185358686`*^9}, {3.626778839767812*^9, 
   3.626778850546812*^9}, {3.626778885279518*^9, 3.6267789004883947`*^9}, {
   3.626778978221078*^9, 3.6267789787074633`*^9}, {3.626978945218978*^9, 
   3.626978970475305*^9}, {3.626979891893981*^9, 3.626979896847517*^9}, 
   3.626979934252185*^9}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"plot1", "=", 
   RowBox[{"ContourPlot3D", "[", 
    RowBox[{
     RowBox[{
      RowBox[{"f", "[", 
       RowBox[{"{", 
        RowBox[{"x", ",", "y", ",", "z"}], "}"}], "]"}], "\[Equal]", "0"}], 
     ",", 
     RowBox[{"{", 
      RowBox[{"x", ",", 
       RowBox[{"-", "3"}], ",", "3"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"y", ",", 
       RowBox[{"-", "3"}], ",", "3"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"z", ",", 
       RowBox[{"-", "3"}], ",", "3"}], "}"}], ",", 
     RowBox[{"AxesLabel", "\[Rule]", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y", ",", "z"}], "}"}]}]}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"plot2", " ", "=", 
   RowBox[{"ContourPlot3D", "[", 
    RowBox[{
     RowBox[{
      RowBox[{
       RowBox[{"grad", "[", 
        RowBox[{"{", 
         RowBox[{"x", ",", "y", ",", "z"}], "}"}], "]"}], "[", 
       RowBox[{"[", "1", "]"}], "]"}], "\[Equal]", "0"}], ",", " ", 
     RowBox[{"{", 
      RowBox[{"x", ",", 
       RowBox[{"-", "3"}], ",", "3"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"y", ",", 
       RowBox[{"-", "3"}], ",", "3"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"z", ",", 
       RowBox[{"-", "3"}], ",", "3"}], "}"}], ",", 
     RowBox[{"AxesLabel", "\[Rule]", 
      RowBox[{"{", 
       RowBox[{"x", ",", "y", ",", "z"}], "}"}]}]}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{"Show", "[", "plot1", "]"}]}], "Input",
 CellChangeTimes->{{3.6267759279297113`*^9, 3.626775969355301*^9}, {
  3.62677603017314*^9, 3.62677604257673*^9}, {3.626776302083514*^9, 
  3.626776345593766*^9}, {3.626778026599371*^9, 3.626778161538435*^9}, {
  3.62677820025244*^9, 3.626778352034292*^9}, {3.626778394091424*^9, 
  3.62677840454002*^9}, {3.626778444674733*^9, 3.62677845586093*^9}, {
  3.626778603602969*^9, 3.626778612203734*^9}, {3.62677890324829*^9, 
  3.626778907896566*^9}, {3.626978988533622*^9, 3.62697903338741*^9}, {
  3.626979961526348*^9, 3.6269799918435287`*^9}, {3.6269800278884373`*^9, 
  3.626980029334762*^9}, {3.626980080025342*^9, 3.626980080888034*^9}}],

Cell["\<\
Here is an example with two equations in three unknowns, an intersection of \
two cylinders.\
\>", "Text",
 CellChangeTimes->{{3.6269801929201727`*^9, 3.6269802428296633`*^9}}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"f", "[", 
   RowBox[{"{", 
    RowBox[{"x_", ",", "y_", ",", "z_"}], "}"}], "]"}], "=", 
  RowBox[{"{", 
   RowBox[{
    RowBox[{
     RowBox[{"x", "^", "2"}], "+", 
     RowBox[{"y", "^", "2"}]}], ",", 
    RowBox[{
     RowBox[{
      RowBox[{"(", 
       RowBox[{"x", "-", "1"}], ")"}], "^", "2"}], "+", 
     RowBox[{"z", "^", "2"}]}]}], "}"}]}], "\[IndentingNewLine]", 
 RowBox[{"a", "=", 
  RowBox[{"{", 
   RowBox[{"1", ",", 
    RowBox[{"Sqrt", "[", "3", "]"}], ",", "1"}], 
   "}"}]}], "\[IndentingNewLine]", 
 RowBox[{"b", "=", 
  RowBox[{"f", "[", "a", "]"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"plot1", "=", 
   RowBox[{"ContourPlot3D", "[", 
    RowBox[{
     RowBox[{"{", 
      RowBox[{
       RowBox[{
        RowBox[{
         RowBox[{"f", "[", 
          RowBox[{"{", 
           RowBox[{"x", ",", "y", ",", "z"}], "}"}], "]"}], "[", 
         RowBox[{"[", "1", "]"}], "]"}], "\[Equal]", 
        RowBox[{"b", "[", 
         RowBox[{"[", "1", "]"}], "]"}]}], ",", 
       RowBox[{
        RowBox[{
         RowBox[{"f", "[", 
          RowBox[{"{", 
           RowBox[{"x", ",", "y", ",", "z"}], "}"}], "]"}], "[", 
         RowBox[{"[", "2", "]"}], "]"}], "\[Equal]", "1"}]}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"x", ",", 
       RowBox[{"-", "3"}], ",", "3"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"y", ",", 
       RowBox[{"-", "3"}], ",", "3"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"z", ",", 
       RowBox[{"-", "3"}], ",", "3"}], "}"}]}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{"Show", "[", 
  RowBox[{"plot1", ",", 
   RowBox[{"Graphics3D", "[", 
    RowBox[{"{", 
     RowBox[{
      RowBox[{"PointSize", "\[Rule]", "Large"}], ",", "Red", ",", 
      RowBox[{"Point", "[", "a", "]"}]}], "}"}], "]"}]}], 
  "]"}], "\[IndentingNewLine]"}], "Input",
 CellChangeTimes->{{3.626980244793541*^9, 3.6269803660483103`*^9}, {
   3.626980414294622*^9, 3.626980457199945*^9}, {3.626980522502513*^9, 
   3.626980652213614*^9}, {3.626980693843371*^9, 3.626980696901113*^9}, {
   3.6269809717323513`*^9, 3.626980974009583*^9}, {3.626981102843829*^9, 
   3.626981128582198*^9}, 3.626981474743848*^9}],

Cell[BoxData[{
 RowBox[{
  RowBox[{"Df", "[", 
   RowBox[{"{", 
    RowBox[{"x_", ",", "y_", ",", "z_"}], "}"}], "]"}], " ", "=", 
  RowBox[{"D", "[", 
   RowBox[{
    RowBox[{"f", "[", 
     RowBox[{"{", 
      RowBox[{"x", ",", "y", ",", "z"}], "}"}], "]"}], ",", 
    RowBox[{"{", 
     RowBox[{"{", 
      RowBox[{"x", ",", "y", ",", "z"}], "}"}], "}"}]}], 
   "]"}]}], "\[IndentingNewLine]", 
 RowBox[{"vec", "=", 
  RowBox[{"First", "[", 
   RowBox[{"NullSpace", "[", 
    RowBox[{"Df", "[", "a", "]"}], "]"}], "]"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"plot2", "=", 
   RowBox[{"ParametricPlot3D", "[", 
    RowBox[{
     RowBox[{"a", "+", 
      RowBox[{"t", " ", "vec"}]}], ",", 
     RowBox[{"{", 
      RowBox[{"t", ",", 
       RowBox[{"-", "1"}], ",", "1"}], "}"}], ",", 
     RowBox[{"PlotRange", "\[Rule]", 
      RowBox[{"{", 
       RowBox[{
        RowBox[{"{", 
         RowBox[{
          RowBox[{"-", "3"}], ",", "3"}], "}"}], ",", 
        RowBox[{"{", 
         RowBox[{
          RowBox[{"-", "3"}], ",", "3"}], "}"}], ",", 
        RowBox[{"{", 
         RowBox[{
          RowBox[{"-", "3"}], ",", "3"}], "}"}]}], "}"}]}], ",", 
     RowBox[{"PlotStyle", "\[Rule]", 
      RowBox[{"{", "Thick", "}"}]}]}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{"Show", "[", 
  RowBox[{"plot1", ",", "plot2", ",", 
   RowBox[{"Graphics3D", "[", 
    RowBox[{"{", 
     RowBox[{
      RowBox[{"PointSize", "[", ".03", "]"}], ",", "Red", ",", 
      RowBox[{"Point", "[", "a", "]"}]}], "}"}], "]"}], ",", 
   RowBox[{"AxesLabel", "\[Rule]", 
    RowBox[{"{", 
     RowBox[{"x", ",", "y", ",", "z"}], "}"}]}]}], 
  "]"}], "\[IndentingNewLine]"}], "Input",
 CellChangeTimes->{{3.626980736471838*^9, 3.626980935589201*^9}, {
   3.6269809791247673`*^9, 3.6269810805169*^9}, {3.626981167715695*^9, 
   3.626981192032774*^9}, {3.626981228105604*^9, 3.62698128580958*^9}, 
   3.626981471371785*^9}],

Cell[BoxData[
 RowBox[{"MatrixForm", "[", 
  RowBox[{"Df", "[", 
   RowBox[{"{", 
    RowBox[{"2", ",", "0", ",", "0"}], "}"}], "]"}], "]"}]], "Input",
 CellChangeTimes->{{3.6269812704868717`*^9, 3.6269812716198606`*^9}, {
  3.6269813062854013`*^9, 3.626981323128317*^9}}]
}, Open  ]]
}, Open  ]]
},
WindowSize->{1920, 1145},
WindowMargins->{{0, Automatic}, {Automatic, 0}},
Magnification->2.,
FrontEndVersion->"10.0 for Linux x86 (64-bit) (June 27, 2014)",
StyleDefinitions->"Default.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[CellGroupData[{
Cell[580, 22, 115, 1, 191, "Title"],
Cell[698, 25, 885, 19, 442, "Text"],
Cell[1586, 46, 355, 11, 61, "Input"],
Cell[1944, 59, 519, 9, 148, "Text"],
Cell[2466, 70, 592, 19, 61, "Input"],
Cell[3061, 91, 187, 6, 65, "Text"],
Cell[3251, 99, 685, 20, 61, "Input"],
Cell[3939, 121, 1345, 24, 610, "Text"],
Cell[5287, 147, 517, 15, 61, "Input"],
Cell[5807, 164, 441, 11, 101, "Input"],
Cell[6251, 177, 521, 16, 61, "Input"],
Cell[6775, 195, 98, 1, 64, "Text"],
Cell[6876, 198, 616, 16, 61, "Input"],
Cell[7495, 216, 126, 1, 64, "Text"],
Cell[7624, 219, 224, 5, 61, "Input"],
Cell[7851, 226, 270, 5, 106, "Text"],
Cell[8124, 233, 1138, 34, 140, "Input"],
Cell[9265, 269, 117, 1, 64, "Text"],
Cell[CellGroupData[{
Cell[9407, 274, 141, 1, 91, "Subsection"],
Cell[9551, 277, 354, 6, 106, "Text"],
Cell[9908, 285, 355, 10, 61, "Input"],
Cell[10266, 297, 343, 6, 106, "Text"],
Cell[10612, 305, 853, 25, 61, "Input"],
Cell[11468, 332, 834, 15, 191, "Text"],
Cell[12305, 349, 474, 10, 129, "Input"],
Cell[12782, 361, 93, 1, 61, InheritFromParent],
Cell[12878, 364, 157, 3, 64, "Text"],
Cell[13038, 369, 1207, 34, 207, "Input"],
Cell[14248, 405, 1710, 51, 285, "Input"],
Cell[15961, 458, 350, 6, 106, "Text"],
Cell[16314, 466, 1658, 53, 179, "Input"],
Cell[17975, 521, 471, 8, 148, "Text"],
Cell[18449, 531, 2458, 72, 218, "Input"],
Cell[20910, 605, 202, 2, 64, "Text"],
Cell[21115, 609, 1961, 48, 179, "Input"],
Cell[23079, 659, 2118, 55, 140, "Input"],
Cell[25200, 716, 186, 4, 64, "Text"],
Cell[25389, 722, 2200, 64, 257, "Input"],
Cell[27592, 788, 1922, 56, 218, "Input"],
Cell[29517, 846, 272, 6, 61, "Input"]
}, Open  ]]
}, Open  ]]
}
]
*)

(* End of internal cache information *)