DE.08  Linearizations
 BASICS

Needs["Graphics`Colors`"];
Needs["Graphics`FilledPlot`"];

If[MemberQ[{"Global`","Graphics`Common`GraphicsCommon`","Vector3D`"}, 
    Context[ScaleFactor]],Remove["*`ScaleFactor"]];
<<"Graphics`Legend`";
<<"Graphics`Common`GraphicsCommon`";
Remove["*`ScaleFactor"];
Graphics`Colors`GosiaGreen=RGBColor[0, 0.392187, 0];
Off[General::spell];
Off[General::spell1];
Off[ParametricPlot3D::ppcom];
Off[ParametricPlot::ppcom];
Off[Plot::plnr];
Needs["Calculus`Master`"];
Off[Unique::"usym"];
SetOptions[Limit,Analytic->True];

CMView = {2.7, 1.6, 1.2};

Off[Graphics::gpt]

New Vector Drawers

(* :Title:  NewVectorDraw  *)
(* :Author: Bruce Carpenter *)
(* :Summary:  Construct 2 and 3 dimensional vectors,
     arrows, and arrowheads as graphics 
     objects.
*)
(* :Date:  Copyright 1997, Math Everywhere, Inc. 
*)
(*  :Drawbacks:  Since the vectors are drawn out of 
     context, the shapes of the arrowheads are only 
     correct with AspectRatio set to Automatic.
*)
BeginPackage["Vector3D`","Graphics`Colors`"];

Axes3D::usage = "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, y, and z running from -a/3 to a, and with axes labels b units beyond the tips of the axes.  Axes3D[a] is Axes3D[a, a/8].";

Perpend::usage = "Perpend[a] returns a unit vector which is perpendicular to the vector a.";

Arrow::usage = "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. Arrow[a, Tail -> tail] gives a vector from tail to a.";
	
Vector::usage = "Vector[a, b] gives a 2D or 3D vector from a to b.";

ArrowHead::usage = "ArrowHead[a, vec] produces an arrowhead with its tip placed at the point a, pointing in the direction of the vector vec.";

Tail::usage = "Tail -> point puts the tail of the vector at point.";

Norm::usage = "Norm[vec] returns the Euclidean length of vec.";

Normalize::usage = "Normalize[vec] returns a unit vector in the direction of vec, provided vec is nonzero.";

Aperture::usage = "Aperture is the ratio of the radius of the base to the length of the head of a vector.";

TipSize::usage = "TipSize specifies an absolute size for the inner tip length of the head of a vector.";

TipRatio::usage = "TipRatio is the ratio of the inner length of the arrowhead to the length of the shaft of a vector.";

HeadRatio::usage = "HeadRatio is the ratio of the outer length of the arrowhead to the length of the shaft of a vector.";

HeadSize::usage = "HeadSize specifies an absolute size for the head of a vector.";

ScaleFactor::usage = "ScaleFactor specifies the amount to scale a vector in length and may be either a number or a pure function.";

ZeroVectorPointSize::usage = "ZeroVectorPointSize is the size of the point used to represent a zero vector.";

Options[Vector] = Options[Arrow] = Options[ArrowHead] =
{HeadRatio -> 0.18, TipRatio -> 0.14, Aperture -> 0.3, ScaleFactor -> 1, HeadType -> Polygon, ShaftQ -> True, EdgesQ -> True, VectorColor -> RGBColor[0,0,1], ShaftWidth -> 0.005, ZeroVectorPointSize -> 0.01};
		
Begin["`Private`"];

SetOptions[ParametricPlot,AspectRatio -> Automatic];
SetOptions[Plot,AspectRatio -> Automatic];
SetOptions[Graphics,AspectRatio -> Automatic];

Axes3D[u_,v_] := 
      Graphics3D[{
            {Blue,Line[{{-u/3,0,0},{u,0,0}}]},
                    Text["x",{u + v,0,0}],
            {Blue,Line[{{0,-u/3,0},{0,u,0}}]},
                    Text["y",{0, u + v,0}],
            {Blue,Line[{{0,0,-u/3},{0,0,u}}]},
                    Text["z",{0,0,u + v}]}]; 
Axes3D[u_] := Axes3D[u,u/8];
				
Norm[a_?VectorQ] := N[Sqrt[a.a]];

Normalize[a_?VectorQ] := If[a.a == 0, a, a/Norm[a]];

Perpend[a_?VectorQ] := Normalize[{a[[2]],-a[[1]]}]/;Length[a]==2
Perpend[a_?VectorQ] := Normalize[If[a[[1]]==0,{1,0,0},{0,a[[3]],-a[[2]]}]]/;Length[a]==3

Base3D = Table[{Cos[2 k Pi/8.], Sin[2 k Pi/8.]},{k,0,8}]; 

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, tip},
   {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, TipRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]; tr = 0.75 hr];
   If[NumberQ[ts], tr = ts/Norm[diff]];
   Graphics[{If[sq,{vc,Thickness[sw],Line[{from,tip - tr diff}]},{}],
      With[{trans = ap hr Norm[diff] Perpend[diff]},
      {vc, Thickness[sw], ht[{tip,
                   (tip - hr diff + trans),
                   (tip - tr diff),
                   (tip - hr diff - trans),tip}]}]}]
] /; Length[from]==Length[to]==2

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, eq, sw, hs, hr, tr,ap, sf, vc, zvps, diff = to - from, tip},
   {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics3D[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]];
   tr = hr/2;
   Graphics3D[{If[sq,{vc,Thickness[sw],Line[{from,tip-tr diff}]},{}], If[eq, {}, EdgeForm[]], {SurfaceColor[vc],
      With[{temp = Perpend[diff]},
      Polygon[Append[#,tip]]& /@ Partition[(tip - hr diff + #)&/@((ap hr Norm[diff] Base3D) . {temp,Normalize[Cross[temp,diff]]}),2,1]]}}]
] /; Length[from]==Length[to]==3


Arrow[a_, Tail->b_, opts___] := Vector[b,a+b,opts]

Arrow[a_, opts___] := Vector[Table[0,{Length[a]}],a,opts]

ArrowHead[a_, b_, opts___] := Vector[a-b,a,ShaftQ->False,opts];

End[];

EndPackage[];
Null

VectorField and Phase Lines

(*  :Title: VectorAndPhase.m  *)

(*  :Author:  Bill Davis  *)

(*  :Dates:  about March, 1996 *)

(*  :Drawbacks:  I wish someone serious would come and turn these into
                 Really carefully done packages.  *)

BeginPackage["PhaseStuff`","Vector3D`",
		"Utilities`FilterOptions`"];

PhaseLine::usage = "PhaseLine[fcn,{var,a,b},tint,pos]
	makes a vertical phase line for the function, fcn, of the variable, var,
	for the variable going from a to b, with 
	a color given by tint,  at the x-coordinate pos.";
	
PhaseLine3D::usage = "PhaseLIne3D[fcn,{var,a,b},tint,loc]
	makes a vertical phase line for the function, fcn, of the variable, var,
	for the variable going from a to b, with 
	a color given by tint, and base point at loc = {x,y}";
	
VectorField::usage = "VectorField[{function,{var1,a,b,space1},{var2,c,d,space2}]
	will make a vector field of normalized vectors, colored by
	magnitude, for the function of the variables on the
	range {r1,r2} with spacing between points given by 
	{space1, space2}.  If a single function of a single variable is given,
	the vector field {1,f} will be plotted.";
	
StikPikXZ::usage = "StikPikXZ[picture,a:_0] takes a 2 dimensional
	Mathematica picture and embeds it on the XZ plane at y = a.  
	Some optional characteristics of the graphs are lost in the 
	process.  StikPikYZ and StikPikXZ are similar.";
	
StikPikYZ::usage = "StikPikYZ[picture,a:_0] takes a 2 dimensional
	Mathematica picture and embeds it on the XZ plane at y = a.  
	Some optional characteristics of the graphs are lost in the 
	process.  StikPikYZ and StikPikXZ are similar.";
	
StikPikXY::usage = "StikPikXY[picture,a:_0] takes a 2 dimensional
	Mathematica picture and embeds it on the XZ plane at y = a.  
	Some optional characteristics of the graphs are lost in the 
	process.  StikPikYZ and StikPikXZ are similar.";
	
Begin["`Private`"];	

PhaseLine[fcn_,{var_,a_,b_},tint_,pos_:0] := 
	Module[{len,len1,k,zset,zero,zero1,signs,stick,dots,pointers},
	zset = Solve[fcn==0,var];
	zero = If[Length[zset]==0,{},var/.zset];
	zero = Sort[zero];
	zero = Select[zero,((# >= a)&&(# <= b))&];
	zero1 = Sort[Union[zero,{a,b}]];
	len = Length[zero];
	len1 = Length[zero1];
	signs = Table[Sign[fcn/.var->((zero1[[k+1]] + zero1[[k]])/2)],
		{k,1,len1-1}];
	stick = {Thickness[0.007],Line[{{pos,a},{pos,b}}]};
	dots = Point[{pos,#}]&/@zero;
	dots = Join[{PointSize[0.05],tint},dots];
	heads = Table[MakeHead[zero1[[k]],zero1[[k+1]],pos,
		signs[[k]],b-a,tint],{k,1,len1-1}];
	Return[{Graphics[Join[{tint},stick,dots]],heads}]
		]
		
MakeHead[a_,b_,c_,s_,slen_,tint_] := Module[{u,v,w},
	u = Min[2(b-a),3/4 slen];
	Return[ArrowHead[{c,a(.5-s/5)+b(.5+s/5)},
		{0,s u},VectorColor->tint]]]
		
PhaseLine3D[fcn_,{var_,a_,b_},tint_,pos_:{0,0}] := 
	Module[{len,len1,k,zset,zero,zero1,signs,stick,dots,pointers},
	zset = Solve[fcn==0,var];
	zero = If[Length[zset]==0,{},var/.zset];
	zero = Sort[zero];
	zero = Select[zero,((# >= a)&&(# <= b))&];
	zero1 = Sort[Union[zero,{a,b}]];
	len = Length[zero];
	len1 = Length[zero1];
	signs = Table[Sign[fcn/.var->((zero1[[k+1]] + zero1[[k]])/2)],
		{k,1,len1-1}];
	stick = {Thickness[0.007],Line[{{pos[[1]],pos[[2]],a},
		{pos[[1]],pos[[2]],b}}]};
	dots = Point[{pos[[1]],pos[[2]],#}]&/@zero;
	dots = Join[{PointSize[0.03],tint},dots];
	heads = Table[MakeHead3D[zero1[[k]],zero1[[k+1]],pos,
		signs[[k]],b-a,tint],{k,1,len1-1}];
	Return[{Graphics3D[Join[{tint},stick,dots]],heads}]
		]
		
MakeHead3D[a_,b_,c_,s_,slen_,tint_] := Module[{u,v,w},
	u = Min[2(b-a),1/2 slen];
	Return[ArrowHead[{c[[1]],c[[2]],a(.5-s/5)+b(.5+s/5)},
		{0,0,s u},VectorColor->tint]]]
		
		
VectorField[function_List,{var1_,a_,b_,sp1_:1},{var2_,c_,d_,sp2_:1}] := 
	Module[{i,j,U,V,UU,VV,ulow,uhigh,vlow,vhigh,maxval,minval,
		lengths,directions,locations,mlen},
	mlen = Min[sp1,sp2];
	normf[uu_,vv_] = Sqrt[function[[1]]^2+function[[2]]^2]/.
		{var1->uu,var2->vv};
	lengths = Table[N[normf[U,V]],{U,a,b,sp1},{V,c,d,sp2}];
	maxval = Max[lengths]; minval = Min[lengths];
	If[minval == maxval,maxval += 1];
	locations = Table[{uu,vv},{uu,a,b,sp1},
		{vv,c,d,sp2}];
	directions = Table[If[N[normf[U,V]]==0,{0,0},
		mlen (function/.{var1->U,var2->V})/(1.3 N[normf[U,V]])],
			{U,a,b,sp1},{V,c,d,sp2}];
	{UU,VV} = Dimensions[directions][[{1,2}]];
	Return[Table[Arrow[directions[[i,j]],
			Tail->locations[[i,j]],Aperture->1,
			HeadRatio->.4,ShaftWidth->0.008, VectorColor->
			RGBColor[(lengths[[i,j]] - minval)/(maxval-minval),0,
				(maxval -lengths[[i,j]])/(maxval-minval)]],{i,UU},{j,VV}]]
		]
		
VectorField[function_,{var1_,a_,b_,sp1_:1},{var2_,c_,d_,sp2_:1}] :=
	VectorField[{1,function},{var1,a,b,sp1},{var2,c,d,sp2}]
	
	
Needs["Utilities`FilterOptions`"]

FixOptions[a_] := Module[{picture = a, opt3},
If[Options[picture] != {}, 
	opt3 = FilterOptions[Graphics3D,#]&/@Options[picture];
	(* Fix things that might carry 2D characteristics. *)
	opt3 = opt3//. {(AxesLabel -> _) ->(AxesLabel -> None),
			(PlotRange->_) -> (PlotRange->Automatic),
			(PlotRegion -> _) -> (PlotRegion -> Automatic),
			(Ticks -> _) -> (Ticks -> Automatic)};
	picture = Apply[List,picture];
	picture[[-1]] = opt3;
	picture = Apply[Graphics3D,picture]];
	Return[picture]
	]

StikPikXZ[picture_,a_:0] := Module[{PIK3 = picture},
	PIK3 = PIK3//.{Graphics->Graphics3D,
		{u_?NumberQ,v_?NumberQ}:>{u,a,v}};
	PIK3 = FixOptions[PIK3];
	Return[PIK3]
		]


StikPikYZ[picture_,a_:0] := Module[{PIK3 = picture},
	PIK3 = PIK3//.{Graphics->Graphics3D,
		{u_?NumberQ,v_?NumberQ}:>{a,u,v}};
	PIK3 = FixOptions[PIK3];
	Return[PIK3]
		]


StikPikXY[picture_,a_:0] := Module[{PIK3 = picture},
	PIK3 = PIK3//.{Graphics->Graphics3D,
		{u_?NumberQ,v_?NumberQ}:>{u,v,a}};
	PIK3 = FixOptions[PIK3];
	Return[PIK3]
		]
		
End[];

EndPackage[];
Null

ThreeAxes[u,v]

ThreeAxes[u_,v_] := 
      Graphics3D[{
            {Blue,Line[{{-u,0,0},{u,0,0}}]},
                    Text["x",{u + v,0,0}],
            {Blue,Line[{{0,-u,0},{0,u,0}}]},
                    Text["y",{0, u + v,0}],
            {Blue,Line[{{0,0,0},{0,0,u}}]},
                    Text["z",{0,0,u + v}]}]; 
ThreeAxes[u_] := ThreeAxes[u,u/8];
ThreeAxes::usage = "ThreeAxes[a,b] makes a standard cartesian axis graphics object
	with x, y, and z running from -a to a, and with axis
	labels b units beyond the tips of the axes.  ThreeAxes[a]
	is ThreeAxes[a,a/8].";

Crisscross

Clear[x];
crisscross = Plot[{Cos[x],Sin[x]},{x,Pi/8,3 Pi/8},
    PlotStyle->{{RGBColor[1,0,0],Thickness[0.02]}},
    Axes->None,AspectRatio->1/3,
    DisplayFunction->Identity];

Pendulum Stuff

modticksX[gr_Graphics,m_] :=
   Module[{f,x},
     f[{a_,"",st___}] := {a,""};
     f[{a_,b_,st___}] := {a,b m};
     Show[gr,
        AxesOrigin->Reverse[FullOptions[gr,AxesOrigin]],
        Ticks->{Map[f,FullOptions[gr,Ticks][[1]]],
                Automatic}]
]

modticksY[gr_Graphics,m_] :=
   Module[{f,x},
     f[{a_,"",st___}] := {a,""};
     f[{a_,b_,st___}] := {a,b m};
     Show[gr,
        AxesOrigin->Reverse[FullOptions[gr,AxesOrigin]],
        Ticks->{Automatic,
                Map[f,FullOptions[gr,Ticks][[2]]]
               }]
]

Clear[center,angle,length,a,l,pendulum,calibratedpendulum];

pendulum[length_,angle_]:=
Module[
{center:={length Sin[Pi-angle],length Cos[Pi-angle]}},

Show[

Graphics[{Blue,Thickness[.01],
Line[{{0,0},{0,-(length+3)}}]}],

Graphics[{Blue,
Polygon[{{-5,-(length+3)},{-5,-(length+3.5)},{-5.3,-(length+3.5)},
{-5.3,-(length+4)},{5.3,-(length+4)},{5.3,-(length+3.5)},
{5,-(length+3.5)},{5,-(length+3)},{-5,-(length+3)}}]}],

Graphics[{Black,
Disk[{0,0},length,{(-Pi/2+angle-Pi/84),(-Pi/2+angle+Pi/84)}]}],

Graphics[{Black,
Disk[center,length,{(Pi/2+angle-Pi/500),(Pi/2+angle+Pi/500)}]}],

Graphics[{Red,PointSize[.05],Point[{0,0}]}],

Graphics[{Red,
Disk[center,length/5,{0,360 Degree}]}],

Graphics[{GoldenrodLight,
Disk[center,length/6,{330 Degree,350 Degree}]}],

Graphics[{Red,
Disk[center,length/7,{0,360 Degree}]}],

Graphics[{Black,PointSize[.02],Point[center]}],


AspectRatio->Automatic,
PlotRange->{{length+5,-(length+3)},{(length+4),-(length+4)}},
AxesOrigin->{0,0},
DisplayFunction->Identity]
]

calibratedpendulum[length_,angle_]:=
Module[
{center:={length Sin[Pi-angle],length Cos[Pi-angle]}},

Show[

Graphics[{Blue,Thickness[.01],
Line[{{0,0},{0,-(length+3)}}]}],

Graphics[{Blue,
Polygon[{{-5,-(length+3)},{-5,-(length+3.5)},{-5.3,-(length+3.5)},
{-5.3,-(length+4)},{5.3,-(length+4)},{5.3,-(length+3.5)},
{5,-(length+3.5)},{5,-(length+3)},{-5,-(length+3)}}]}],

Graphics[{Black,
Disk[{0,0},length,{(-Pi/2+angle-Pi/84),(-Pi/2+angle+Pi/84)}]}],

Graphics[{Black,
Disk[center,length,{(Pi/2+angle-Pi/500),(Pi/2+angle+Pi/500)}]}],

Graphics[{Red,PointSize[.05],Point[{0,0}]}],

Graphics[{Red,
Disk[center,length/5,{0,360 Degree}]}],

Graphics[{GoldenrodLight,
Disk[center,length/6,{330 Degree,350 Degree}]}],

Graphics[{Red,
Disk[center,length/7,{0,360 Degree}]}],

Graphics[{Black,PointSize[.02],Point[center]}],

Graphics[Text[StringForm["y[t] =  ``",angle],
{8 length/8 Sin[15 Pi/12], 8 length/8 Cos[5Pi/6]}]],

AspectRatio->Automatic,
PlotRange->{{length+5,-(length+3)},{(length+4),-(length+4)}},
AxesOrigin->{0,0},DisplayFunction->Identity]
]

B.1) Approximating a nonlinear diffeq system by linearizing at equilibrium points: 
        How to do it and how to use it

B.1.a.i)

Here's a nonlinear diffeq system:

\(Clear[m, n, x, y, t]; 
m[x_, y_]\  = 0.4\ x - 0.9\ y\^3 + 1.2; 

    n[x_, y_]\  = \ \ x\  - \ 0.4\ Sin[y]; 

    diffeqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[diffeqsystem]]\)

And a flow plot:

\(Clear[Field]; 
Field[x_, y_] = {m[x, y], n[x, y]}; 

    \({xlow, xhigh} = {\(-2\), 2}; \)\){ylow, yhigh} = {-2, 2}; jump = ((xhigh - xlow)/12); flowplot = 
      Table[Arrow[Field[x, y], Tail -> {x, y}, VectorColor -> Blue, 
          ScaleFactor -> 0.25, HeadSize -> 0.2], {x, xlow, xhigh, 
          jump}, {y, ylow, yhigh, jump}]; nonlinearflow = 
      Show[flowplot, Axes -> True, 
        AxesLabel -> {"x", "y"}];

You can spot an equilibrium point in this plot.
Plot it.

Saying that a point {xeq, yeq} is an equilibrium point is the same as saying that there is no flow out of {xeq, yeq}, and this is the same as saying:
            {m[xeq, yeq], n[xeq, yeq]} = {0, 0}.

This tells you that {xeq, yeq} lies at the intersection of the curve defined by
           m[x, y] = 0 
 and the curve defined by 
           n[x, y] = 0.
Here are these two curves shown with the flow plot:

mcontourplot = 
      ContourPlot[m[x, y], {x, xlow, xhigh}, {y, ylow, yhigh}, 
        Contours -> {0}, ContourStyle -> Thickness[0.01], 
        PlotPoints -> 50, ContourSmoothing -> Automatic, 
        ContourShading -> False, DisplayFunction -> Identity]; 
    ncontourplot = 
      ContourPlot[n[x, y], {x, xlow, xhigh}, {y, ylow, yhigh}, 
        Contours -> {0}, ContourStyle -> Thickness[0.01], 
        PlotPoints -> 50, ContourSmoothing -> Automatic, 
        ContourShading -> False, DisplayFunction -> Identity]; 
    Show[nonlinearflow, mcontourplot, ncontourplot, 
      DisplayFunction -> $DisplayFunction];

You can see  {xeq, yeq}  near  {0.2, 1.3}.
To nail {xeq, yeq} down, do this:

{xeq, yeq}={x, y}/.FindRoot[{m[x, y] == 0, n[x, y] == 0},{x,0.2},{y,1.3}]

See the equilibrium point with the flow:

\(equilibriumplot = Graphics[{Red, PointSize[0.04], Point[{xeq, yeq}]}]; 
    

Show[nonlinearflow, equilibriumplot]; \)

B.1.a.ii) Using the linearization matrix to linearize a nonlinear diffeq system

Stay with the same nonlinear diffeq system as in part i);

\(Clear[m, n, x, y, t]; 
m[x_, y_]\  = \ 0.4\ x - 0.9\ y\^3 + 1.2; 

    n[x_, y_]\  = \ \ x\  - \ 0.4\ Sin[y]; 

    diffeqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[diffeqsystem]]\)

The linearization matrix A[x, y] at {x, y} for this system is:

\(Clear[A, gradm, gradn]; 

    gradm[x_, y_]\  = \ {D[m[x, y], x], D[m[x, y], y]}; 

    gradn[x_, y_]\  = \ {D[n[x, y], x], D[n[x, y], y]}; 

    A[x_, y_]\  = \ {gradm[x, y], gradn[x, y]}; 
MatrixForm[A[x, y]]\)

How do you use the linearization matrix A[xeq, yeq] to linearize this system at {xeq, yeq}?

Very easily.
You go with the linearization matrix A[xeq, yeq] and use it for the coefficient matrix of this system:
       {x'[t],y'[t]}=A[xeq, yeq] . {x[t] - xeq, y[t] - yeq}:

\(linearized\  = \ 

      \(({\(x'\)[t], \(y'\)[t]}\  == \ 
          Chop[Expand[A[xeq, yeq] . {x[t]\  - \ xeq, y[t]\  - \ yeq}]])\); 

    ColumnForm[Thread[linearized]]\)

Done.
You can see A[xeq, yeq] sitting in the system:

MatrixForm[A[xeq, yeq]]

B.1.b)

Here is the flow for the original system:

h = 1.0; {xlow, xhigh} = {xeq - h, xeq + h}; {ylow, yhigh} = {yeq - h, yeq + h}; jump = ((xhigh - xlow)/12); scalefactor = 0.2; sizer = 0.1; nonlinearflow = 
      Table[Arrow[Field[x, y], Tail -> {x, y}, VectorColor -> Blue, 
          ScaleFactor -> scalefactor, HeadSize -> sizer], {x, xlow, 
          xhigh, jump}, {y, ylow, yhigh, jump}]; nonlinearflowplot = 
      Show[nonlinearflow, equilibriumplot, 
        PlotRange 
          -> {{xeq - 1.5 h, xeq + 1.5 h}, {yeq - 1.5 h, 
              yeq + 1.5 h}}, Axes -> True, 
        AxesLabel -> {"x", "y"}, 
        PlotLabel -> "Nonlinear"];

Here is the corresponding flow of the linearized system:

\(Clear[linearizedField]; 

    linearizedField[x_, y_] = A[xeq, yeq] . {x - xeq, y - yeq}; 

    \(linearizedflow = 
      Table[Arrow[linearizedField[x, y], Tail -> {x, y}, 
          VectorColor -> Red, ScaleFactor -> scalefactor, 
          HeadSize -> sizer], {x, xlow, xhigh, jump}, {y, ylow, yhigh, 
          jump}]; \)\)linearizedflowplot = 
      Show[linearizedflow, equilibriumplot, 
        PlotRange 
          -> {{xeq - 1.5 h, xeq + 1.5 h}, {yeq - 1.5 h, 
              yeq + 1.5 h}}, Axes -> True, 
        AxesLabel -> {"x", "y"}, 
        PlotLabel -> "Linearized at equilibrium"];

Here are both flows together:

both = 
      Show[nonlinearflowplot, linearizedflowplot, 
        PlotLabel -> "Both"];

What fundamental relation between the orginal nonlinear system and its linearization at the equilibrium point is illustrated by these plots?

Take another look:

Show[both];

Near the equilibrium point, the flow of the linearization is a deadringer for the original nonlinear flow. As you go further from the equilbrium point, the quality of the approximation begins to deteriorate.
See what is happening in the immediate vicinty of the equilibrium point:

h = 0.1; {xlow, xhigh} = {xeq - h, xeq + h}; {ylow, yhigh} = {yeq - h, yeq + h}; jump = ((xhigh - xlow)/8); scalefactor = 0.2; sizer = 0.01; nonlinearflowplot = 
      Table[Arrow[Field[x, y], Tail -> {x, y}, VectorColor -> Blue, 
          ScaleFactor -> scalefactor, HeadSize -> sizer], {x, xlow, 
          xhigh, jump}, {y, ylow, yhigh, jump}]; linearizedflowplot = 
      Table[Arrow[linearizedField[x, y], Tail -> {x, y}, 
          VectorColor -> Red, ScaleFactor -> scalefactor, 
          HeadSize -> sizer], {x, xlow, xhigh, jump}, {y, ylow, yhigh, 
          jump}]; Show[nonlinearflowplot, linearizedflowplot, equilibriumplot, 
      PlotRange 
        -> {{xeq - 1.5 h, xeq + 1.5 h}, {yeq - 1.5 h, yeq + 1.5 h}}, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      PlotLabel -> "Both"];

Near the equilibrium point, you get a great approximation of the nonlinear flow by the linearized flow.
See what happens when you go far away:

h = 5; {xlow, xhigh} = {xeq - h, xeq + h}; {ylow, yhigh} = {yeq - h, yeq + h}; jump = ((xhigh - xlow)/8); scalefactor = 0.1; sizer = 0.5; nonlinearflowplot = 
      Table[Arrow[Field[x, y], Tail -> {x, y}, VectorColor -> Blue, 
          ScaleFactor -> scalefactor, HeadSize -> sizer], {x, xlow, 
          xhigh, jump}, {y, ylow, yhigh, jump}]; linearizedflowplot = 
      Table[Arrow[linearizedField[x, y], Tail -> {x, y}, 
          VectorColor -> Red, ScaleFactor -> scalefactor, 
          HeadSize -> sizer], {x, xlow, xhigh, jump}, {y, ylow, yhigh, 
          jump}]; Show[nonlinearflowplot, linearizedflowplot, equilibriumplot, 
      PlotRange 
        -> {{xeq - 1.5 h, xeq + 1.5 h}, {yeq - 1.5 h, yeq + 1.5 h}}, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      PlotLabel -> "Both"];

The further you go from the equilibrium point, the lower the quality of the approximation of the original nonlinear system by the linearized system.

B.1.c) Using the linearization matrix to predict the behavior of a nonlinear system

What advantage do you get by approximating this nonlinear system with its linearized version?

You can look at the eigenvalues of the linearization matrix:

Eigenvalues[A[xeq, yeq]]

Propelling swirlers.
Because the linerized flow approximates the original nonlinear flow so well, this information guarantees that trajectories in the original nonlinear system that start near (but not on) the equilibrium point spiral away from the equilibrium point never to return to the vicinity of equilibrium point.
Watch it happen for some random starting points near the equilibrium point:

h = 0.3; {xstarter, ystarter} = {Random[Real, {xeq - h, xeq + h}], 
        Random[Real, {yeq - h, yeq + h}]}; endtime = 14; starterpoint = {xstarter, ystarter}; Clear[x, y];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; approxsolutions = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; \(Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. approxsolutions[[1]], 
        y[t] /. approxsolutions[[1]]}; 

    
trajectoryplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 


    starterplot = Graphics[{Red, PointSize[0.03], Point[starterpoint]}]; 


    Show[starterplot, trajectoryplot, equilibriumplot, PlotRange -> All, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      PlotLabel -> "Nonlinear trajectory", 
      DisplayFunction -> $DisplayFunction]; \)

Just as predicted by the eigenvalues of the linearization matrix (Jacobian matrix) at the equillibrium point:

Eigenvalues[A[xeq, yeq]]

B.2) Lyapunov's rules

B.2.a) Lyapunov's

Given a nonlinear system, when can you trust the prediction you get from the eigenvalues of the linearization matrix A[xeq, yeq] at an equilibrium point {xeq, yeq}?

A juicy piece of advanced mathematics covers this.  A brillant Russian fellow named Lyapunov proved conclusively the following rules:

Attractor rules:

If the eigenvalues of the linearization matrix A[xeq, yeq] are both negative real numbers, then trajectories of the nonlinear system that start near the equilibrium point are sucked more or less directly to the equilibrium point. 

If the eigenvalues of the linearization matrix A[xeq, yeq] are 
        p + I q and p - I q with p < 0 and q ≠ 0, 
then trajectories of the nonlinear system that start near the equilibrium point spiral in and stall at the equilibrium point.

If either of these happens, most folks say the equilbrium point is an attractor.

Repellor rules:

If the eigenvalues of the linearization matrix A[xeq, yeq] are both positive real numbers, then trajectories of the nonlinear system that start near the equilibrium point are propelled more or less directly from the equilibrium point. 

If the eigenvalues of the linearization matrix A[xeq, yeq] are 
        p + I q and p - I q with p > 0 and q ≠ 0, 
then trajectories of the nonlinear system that start near the equilibrium point spiral away from the equilibrium point.

If either of these happens, most folks say the equilbrium point is a repellor.

Saddle rule:

If the eigenvalues of the linearization matrix A[xeq, yeq] consist of one positive real number and one negative real number, then trajectories of the nonlinear system that start near the equilibrium are either propelled directly away from the equilibrium point or move toward the equilibrium point and then move away.
If either of these happens, most folks say the equilbrium point is a saddle.

No information rule:

If the equilibrium point is neither a repellor nor an attractor nor a saddle, then you are walking on thin ice if you trust the information you get from the linearization at the equilibrium point.

B.2.b)

Give some illustrations of Lyapunov's rules.

Here's one:

\(Clear[m, n, x, y, t]; 

    m[x_, y_]\  = \(-1.3\)\ x + 0.17\ y - 0.2\ x\^2\ y; 

    n[x_, y_]\  = \ 0.50\ x\  - \ 0.29\ y\  + \ 0.7; 

    DiffEqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[DiffEqsystem]]\)

The linearization matrix (Jacobian) at a point {x, y} is:

\(Clear[A, gradm, gradn]; 

    gradm[x_, y_] = {∂\_x m[x, y], ∂\_y m[x, y]}; 

    \(gradn[x_, y_] = {∂\_x n[x, y], ∂\_y n[x, y]}; 
    \)\)A[x_, y_] = {gradm[x, y], gradn[x, y]}; MatrixForm[A[x, y]]

The equilibrium points are:

equisols  =  
      ColumnForm[Solve[{m[x, y]  ==  0, n[x, y]  ==  0}, {x, y}]]

The only equilbrium point is:

{xeq, yeq} = {0.339065, 2.99839}

The eigenvalues of the linearization matrix at {xeq, yeq} are:

Eigenvalues[A[xeq, yeq]]

Both negative. The equilibrium point is an attractor.
See a random trajectory that starts near the equilibrium point head right for the equilibrium point:

h = 2; {xstarter, ystarter} = {Random[Real, {xeq - h, xeq + h}], 
        Random[Real, {yeq - h, yeq + h}]}; endtime = 15; starterpoint = {xstarter, ystarter};

Clear[x, y];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; \(approxsolutions = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. approxsolutions[[1]], 
        y[t] /. approxsolutions[[1]]}; 

    
trajectoryplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 

    \(starterplot = Graphics[{Red, PointSize[0.05], Point[starterpoint]}]; 
    \)\)\(equilibriumplot = 
      Graphics[{Red, PointSize[0.05], Point[{xeq, yeq}]}]; 

    starterpointlabel = 
      Graphics[{Red, Text["start", starterpoint, {0, 1.5}]}]; 

    equiliblabel = 
      Graphics[{Blue, Text["equilib point", {xeq, yeq}, {0, 1.5}]}]; 


    Show[starterplot, trajectoryplot, equilibriumplot, equiliblabel, 
      starterpointlabel, 
      PlotRange -> {{xeq - h, xeq + h}, {yeq - h, yeq + h}}, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      AspectRatio -> 1, PlotLabel -> "Nonlinear trajectory", 
      DisplayFunction -> $DisplayFunction]; \)

Down the drain every time.

Here's another:

\(Clear[m, n, x, y, t]; 
m[x_, y_]\  = \ 1.3\ y - 8.3\ x\^2; 

    n[x_, y_]\  = \ \ \(-1.6\)\ x\ ; 

    DiffEqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[DiffEqsystem]]\)

The linearization matrix at a point {x, y} is:

\(Clear[A, gradm, gradn]; 

    gradm[x_, y_] = {∂\_x m[x, y], ∂\_y m[x, y]}; 

    \(gradn[x_, y_] = {∂\_x n[x, y], ∂\_y n[x, y]}; 
    \)\)A[x_, y_] = {gradm[x, y], gradn[x, y]}; MatrixForm[A[x, y]]

The equilibrium points are:

equisols  =  
      ColumnForm[Solve[{m[x, y]  ==  0, n[x, y]  ==  0}, {x, y}]]

The equilbrium point is {0, 0} :

{xeq, yeq} = {0, 0}

The eigenvalues of the linearization matrix at {xeq, yeq} are:

Eigenvalues[A[xeq, yeq]]

The linearization predicts trajectories that oscillate around {xeq, yeq}. 
But Lyapunov's "No Information" rule says not to trust this.

See a trajectory that starts near {xeq, yeq}:

{xstarter, ystarter} = {xeq, yeq}  +  {1, 1};
endtime = 6.2;
starterpoint = {xstarter, ystarter};

Clear[x, y];
xdiffeq=x'[t]==m[x[t], y[t]];
ydiffeq=y'[t]==n[x[t], y[t]];xstartereqn = x[0] == xstarter; ystartereqn = y[0] == ystarter; \(approxsolutions = 
      NDSolve[{xdiffeq, ydiffeq, xstartereqn, ystartereqn}, {x[t], y[t]}, {t, 
          0, endtime}]; 

Clear[x, y, t]; 
    
{x[t_], y[t_]} = {
        x[t] /. approxsolutions[[1]], 
        y[t] /. approxsolutions[[1]]}; 

    
trajectoryplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, endtime}, 
        PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
        DisplayFunction -> Identity]; 

    \(starterplot = Graphics[{Red, PointSize[0.05], Point[starterpoint]}]; 
    \)\)\(equilibriumplot = 
      Graphics[{Red, PointSize[0.05], Point[{xeq, yeq}]}]; 

    starterpointlabel = 
      Graphics[{Red, Text["start", starterpoint, {0, 1.5}]}]; 

    equiliblabel = 
      Graphics[{Blue, Text["equilib point", {xeq, yeq}, {0, 1.5}]}]; 


    Show[starterplot, trajectoryplot, equilibriumplot, equiliblabel, 
      starterpointlabel, PlotRange -> {{\(-1\), 1}, {\(-1\), 1}}, 
      Axes -> True, AxesLabel -> {"x", "y"}, 
      AspectRatio -> 1, PlotLabel -> "Nonlinear trajectory", 
      DisplayFunction -> $DisplayFunction]; \)

The linearization predicts trajectories that oscillate around {xeq, yeq}.  
But this time, the information from the linearization matrix was useless.
This is in harmony with Lyapunov's rules.

B.3) The pendulum oscillator: Damped and undamped

Put your eyes on this snapshot of a pendulum:

L = 10; setup = 
      Show[pendulum[L, (π/6)], DisplayFunction -> $DisplayFunction];

It's just like the pendulum that hangs out the bottom of a cuckoo clock, or the bigger pendulum that hangs in that old clock in your grandmother's hallway. 
To come up with a tractable model for the motion of the pendulum, you have to make a few assumptions that may be hard to swallow at first.  
They are these: 
→ The pendulum hangs from a frictionless pivot. 
→ The pendulum moves in a vacuum so that there is no air resistance. 
→ The arm of the pendulum which holds the bob at the end is both massless and perfectly rigid. 
With those assumptions under your belt, look at this:

\(Clear[y, t, s]; 
L = 10; 
\(y[t] = π\/6; \)\)\(Clear[s]; 

    angle = ParametricPlot[0.5\ L\ {Sin[s], \(-Cos[s]\)}, {s, 0, y[t]}, 
        DisplayFunction -> Identity]; 

    \(labels = 
      Graphics[{
          Text["y[t]", 1\/10\ 5.5\ L\ {Sin[y[t]\/2], \(-Cos[y[t]\/2]\)}], 
          Text["L", 
            0.4\ L\ {Sin[y[t]], \(-Cos[y[t]]\)} + {Cos[y[t]], Sin[y[t]]}]}]; 
    \)\)Show[calibratedpendulum[L, y[t]], angle, labels, 
      DisplayFunction -> $DisplayFunction];

Here's what happens when you go with
        y[t] = (π/3):

\(Clear[y, t, s]; 
y[t] = π\/3; 
Clear[s]; 

    angle = ParametricPlot[0.5\ L\ {Sin[s], \(-Cos[s]\)}, {s, 0, y[t]}, 
        DisplayFunction -> Identity]; 

    \(labels = 
      Graphics[{
          Text["y[t]", 1\/10\ 5.5\ L\ {Sin[y[t]\/2], \(-Cos[y[t]\/2]\)}], 
          Text["L", 
            0.4\ L\ {Sin[y[t]], \(-Cos[y[t]]\)} + {Cos[y[t]], Sin[y[t]]}]}]; 
    \)\)Show[calibratedpendulum[L, y[t]], angle, labels, 
      DisplayFunction -> $DisplayFunction];

Here's what happens when you go with
         y[t] = 0:

\(Clear[y, t, s]; 
y[t] = 0; 
Clear[s]; 

    labels = Graphics[{
          Text["L", 
            0.4\ L\ {Sin[y[t]], \(-Cos[y[t]]\)} + {Cos[y[t]], Sin[y[t]]}]}]; 
    
\(Show[calibratedpendulum[L, y[t]], labels, 
      DisplayFunction -> $DisplayFunction]; \)\)

Here's what happens when you go with
        y[t] = -(π/3):

\(Clear[y, t, s]; 
y[t] = \(-\(π\/3\)\); 
Clear[s]; 

    angle = ParametricPlot[0.5\ L\ {Sin[s], \(-Cos[s]\)}, {s, 0, y[t]}, 
        DisplayFunction -> Identity]; 

    \(labels = 
      Graphics[{
          Text["y[t]", 1\/10\ 5.5\ L\ {Sin[y[t]\/2], \(-Cos[y[t]\/2]\)}], 
          Text["L", 
            0.4\ L\ {Sin[y[t]], \(-Cos[y[t]]\)} + {Cos[y[t]], Sin[y[t]]}]}]; 
    \)\)Show[calibratedpendulum[L, y[t]], angle, labels, 
      DisplayFunction -> $DisplayFunction];

Here's what happens when you go with
      y[t] = ((3 π)/2):

\(Clear[y, t, s]; 
y[t] = \(3\ π\)\/2; 
Clear[s]; 

    angle = ParametricPlot[0.5\ L\ {Sin[s], \(-Cos[s]\)}, {s, 0, y[t]}, 
        DisplayFunction -> Identity]; 

    \(labels = 
      Graphics[{
          Text["y[t]", 1\/10\ 5.5\ L\ {Sin[y[t]\/2], \(-Cos[y[t]\/2]\)}], 
          Text["L", 
            0.4\ L\ {Sin[y[t]], \(-Cos[y[t]]\)} + {Cos[y[t]], Sin[y[t]]}]}]; 
    \)\)Show[calibratedpendulum[L, y[t]], angle, labels, 
      DisplayFunction -> $DisplayFunction];

Get it?  Got it?  Good!
To get ready to come up with the differential equation y[t] solves, look at this set up:

bobpath = 
      ParametricPlot[10 {Sin[s], -Cos[s]}, {s, -π, π}, 
        DisplayFunction -> Identity]; \(Clear[y, t, s]; 
L = 10; 
\(y[t] = 0.3\ π; \)\)\(Clear[s]; 

    angle = ParametricPlot[0.5\ L\ {Sin[s], \(-Cos[s]\)}, {s, 0, y[t]}, 
        DisplayFunction -> Identity]; 

    \(labels = 
      Graphics[{
          Text["y[t]", 1\/10\ 5.5\ L\ {Sin[y[t]\/2], \(-Cos[y[t]\/2]\)}], 
          Text["L", 
            0.4\ L\ {Sin[y[t]], \(-Cos[y[t]]\)} + {Cos[y[t]], Sin[y[t]]}]}]; 
    \)\)setup = 
      Show[pendulum[L, y[t]], bobpath, angle, labels, 
        DisplayFunction -> $DisplayFunction];

That circle is the path the center of the bob moves on.

B.3.a.i) The pendulum oscillator diffeq 
     y''[t]+((g Sin[y[t]])/L)=0.

Use the picture and Newton's law 
         Force = Mass × Acceleration
to come up with the differential equation that models the motion of this pendulum.

Take another look at the set up together with the circular path on which the center of the lump swings:

bobpath = 
      ParametricPlot[10 {Sin[s], -Cos[s]}, {s, -π, π}, 
        DisplayFunction -> Identity]; \(Clear[y, t, s]; 
L = 10; 
\(y[t] = 0.3\ π; \)\)\(Clear[s]; 

    angle = ParametricPlot[0.5\ L\ {Sin[s], \(-Cos[s]\)}, {s, 0, y[t]}, 
        DisplayFunction -> Identity]; 

    \(labels = 
      Graphics[{
          Text["y[t]", 1\/10\ 5.5\ L\ {Sin[y[t]\/2], \(-Cos[y[t]\/2]\)}], 
          Text["L", 
            0.4\ L\ {Sin[y[t]], \(-Cos[y[t]]\)} + {Cos[y[t]], Sin[y[t]]}]}]; 
    \)\)setup = 
      Show[pendulum[L, y[t]], bobpath, angle, labels, 
        DisplayFunction -> $DisplayFunction];

Agree that the length of the rod is L.
You need a circular parameterization of the circle in terms of y[t].  Here's one that works very well:

Clear[position, L, y, t]; position[t_] = L {Sin[y[t]], -Cos[y[t]]}

This is the parametrization of the circle you want to make the pendulum hang down in the vertical position when the angle y[t] = 0:

position[t] /. y[t] -> 0

Take a new look with all the goodies thrown in:

bobpath = 
      ParametricPlot[10 {Sin[s], -Cos[s]}, {s, -π, π}, 
        DisplayFunction -> Identity]; \(Clear[y, t, s]; 
L = 10; 
\(y[t] = 0.3\ π; \)\)\(Clear[s]; 

    angle = ParametricPlot[0.5\ L\ {Sin[s], \(-Cos[s]\)}, {s, 0, y[t]}, 
        DisplayFunction -> Identity]; 

    \(labels = 
      Graphics[{
          Text["y[t]", 1\/10\ 5.5\ L\ {Sin[y[t]\/2], \(-Cos[y[t]\/2]\)}], 
          Text["L", 
            0.4\ L\ {Sin[y[t]], \(-Cos[y[t]]\)} + {Cos[y[t]], Sin[y[t]]}], 
          Text["M g", {L\ Sin[y[t]], \(-L\)\ Cos[y[t]]} + {0, 
                \(-0.3\)\ L} + {0.15\ L, 0}], 
          Text["F[t]", 
            0.85\ {L\ Sin[y[t]], \(-L\)\ Cos[y[t]]} + 
              0.5\ {0, \(-0.6\)\ L} . {Cos[y[t]], Sin[y[t]]}\ {Cos[y[t]], 
                  Sin[y[t]]}]}]; \)\)overallforce = 
      Arrow[{0, -0.6 L}, 
        Tail -> {L Sin[y[t]], -L Cos[y[t]]}]; tangentialforce = 
      Arrow[{0, -0.6 L} . {Cos[y[t]], Sin[y[t]]} {Cos[y[t]], 
            Sin[y[t]]}, Tail -> {L Sin[y[t]], -L Cos[y[t]]}]; 
    factor = 
      (Sqrt[(((0.6 L))^2) - 
          ((({0, -0.6 L} . {Cos[y[t]], Sin[y[t]]}))^2)]); tension = 
      Arrow[factor {Sin[y[t]], -Cos[y[t]]}, 
        Tail -> {L Sin[y[t]], -L Cos[y[t]]}]; all = 
      Show[pendulum[L, y[t]], angle, labels, bobpath, overallforce, 
        tangentialforce, tension, DisplayFunction -> $DisplayFunction];

Agree that:
→ The mass of the bob on the end is M. 
→ The vertical force on the bob is M(ass) g(ravity).  
→ The component of force pointing in the direction of the pendulum rod can only try to stretch the unstretchable rod,
    so the only part of the force that actually moves the pendulum is the tangential component F[t].
Take another look:

Show[all];

In this snapshot, assume the bob is moving from left to right so that 
           y'[t]>0. 
And look at:

Clear[L, t, y];
position'[t]

Because
           y'[t]>0, 
the unit tangent vector pointing in the direction of the path of the center of the bob is:

Clear[unittan]; unittan[t_] = {Cos[y[t]], Sin[y[t]]}

The unit tangent points opposite to the tangential acceleration. Consequently the tangential component of the force on the bob, F[t], is given by:

Clear[F, M]; F[t_]=Simplify[-M position''[t].unittan[t]]

Take one final peek at the picture:

Show[all];

The picture tells you that F[t] is also given by:

\(Clear[sameF, g]; 
sameF[t_] = M\ g\ Sin[y[t]]\)

The result is:

pendulummodel = F[t] == sameF[t]

Solve this for y''[t]:

Solve[pendulummodel,y''[t]]

This gives:

pendulumoscillator= (y''[t]+((g Sin[y[t]])/L)==0)

There it is! 
A very famous differential equation.

Look again at:

position'[t]

When you assume the bob is moving from right to left, then you have to say that y'[t]<0.
In this case, the unit tangent vector pointing in the direction of the path of the center of the bob is:

Clear[unittan]; unittan[t_] = -{Cos[y[t]], Sin[y[t]]}

Take another look at the snapshot:

Show[all];

Because the bob is moving from right to left, the unit tangent above points in the same direction as the tangential component of the acceleration. This gives

Clear[F, M]; F[t_] = Simplify[M position''[t] . unittan[t]]

On the otherhand, the picture still tells you that F[t] is also given by:

\(Clear[sameF, g]; 
sameF[t_] = M\ g\ Sin[y[t]]\)

The result is:

pendulummodel = F[t] == sameF[t]

Solve[pendulummodel,y''[t]]

This gives you:

pendulumoscillator=y''[t]+((g Sin[y[t]])/L)==0

Just as you got in the case that the pendulum is moving from left to right.
Going either way, the pendulum oscillator is a winner.

B.3.a.ii) The damped pendulum oscillator
        y''[t]+((b y'[t])/(L M))+((g Sin[y[t]])/L)=0.

What differential equation do you get when you include a friction term?

You do exactly as you did above at the start.

Show[all];

Clear[position, L, x, t]; position[t_] = L {Sin[y[t]], -Cos[y[t]]}

position'[t]

Assume the bob is moving from left to right so that   y'[t]>0.
Accordingly the unit tangent vector pointing in the direction of the path of the center of the bob is:

Clear[unittan]; unittan[t_] = {Cos[y[t]], Sin[y[t]]}

Without friction, you put:

Clear[F, M]; F[t_]=Simplify[-M position''[t].unittan[t]]

With friction, you put:

Clear[F, M, b]; F[t_]=(F[t_]=Simplify[-M position''[t].unittan[t]])-b y'[t]

Now you go on as in part i):

\(Clear[sameF, g]; 
sameF[t_] = M\ g\ Sin[y[t]]\)

The result is:

dampedpendulummodel = F[t] == sameF[t]

Solve[dampedpendulummodel,  y''[t]]

A little rearranging and cancelling gives:

dampedpendulumoscillator=y''[t]+((b y'[t])/(L M))+((g Sin[y[t]])/L)==0

There it is!
A very famous nonlinear differential equation.

B.4) The ordinary linear oscillator 
                 y''[t]+((g y[t])/L)=0
        is a linearization of the pendulum oscillator    
                 y''[t]+((g Sin[y[t]])/L)=0

B.4.a)

The undamped pendulum oscillator diffeq is:

Clear[y, g, L, t];
pendulumoscillator=y''[t]+((g Sin[y[t]])/L)==0

The ordinary undamped linear oscillator diffeq is:

Clear[y, g, L, t];
ordinaryoscillator=y''[t]+((g y[t])/L)==0

Explain why the ordinary linear oscillator arrives as a linearization of the pendulum oscillator.

Look at:

pendulumoscillator

Make this second order differential equation into a system of two first order differential equations, by putting 
             x[t]=y'[t]
and then replacing y'[t]with x[t] and replacing y''[t]with x'[t]:

\(Clear[x]; 

    ColumnForm[
      Thread[{new\  = \ \((\(y'\)[t]\  == \ x[t])\), 

          pendulumoscillator /. {\(y'\)[t] -> x[t], 
              \(\(y'\)'\)[t] -> \(x'\)[t]}}]\ ]\)

Clean this up:

\(Clear[m, n, x, y, t]; 
m[x_, y_]\  = \ \(-\(\(g\ Sin[y]\)\/L\)\); 

    n[x_, y_]\  = \ \ x; 

    pendulumsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[pendulumsystem]]\)

The linearization matrix at {0, 0} is:

{xeq, yeq} = {0, 0}; \(Clear[A, gradm, gradn]; 

    gradm[x_, y_] = {∂\_x m[x, y], ∂\_y m[x, y]}; 

    \(gradn[x_, y_] = {∂\_x n[x, y], ∂\_y n[x, y]}; 
    \)\)A[x_, y_] = {gradm[x, y], gradn[x, y]}; MatrixForm[A[xeq, yeq]]

The linearization at {0, 0} is:

\(linearized\  = \ 

      \(({\(x'\)[t], \(y'\)[t]}\  == \ 
          Expand[A[xeq, yeq] . {x[t]\  - \ xeq, y[t]\  - \ yeq}])\); 

    ColumnForm[Thread[linearized]]\)

Remembering that x[t]=y'[t], you see that this is the same as:

y''[t]+((g y[t])/L)==0

The upshot:
The ordinary linear oscillator
        y''[t]+((g y[t])/L)=0
is a linearization of the pendulum oscillator
        y''[t]+((g Sin[y[t]])/L)=0.

B.4.b.i) Using the linearization to see why the pendulum oscillator sometimes behaves
            like an ordinary linear oscillator and sometimes it doesn't

Here are three solutions of the pendulum oscillator corresponding to
         g = 9.8 and L = 5:
They all start with y[0] = 0.8, but one starts with a reasonably small y'[0], the next one starts with a bigger y'[0] and the last starts with an even bigger  y'[0]:

Clear[y, yslow, yfaster, yreallyfast];
pendulumoscillator=y''[t]+((g Sin[y[t]])/L)==0;
g = 9.8; 
      L = 5; endtime = 13; slowstart=NDSolve[{pendulumoscillator,y[0] == 0.8,y'[0]==0.5},y[t],{t, 0, endtime}];yslow[t_] = y[t] /. slowstart[[1]]; 
    fasterstart=NDSolve[{pendulumoscillator,y[0] == 0.8,y'[0]==2.57},y[t],{t, 0, endtime}];yfaster[t_] = 
      y[t] /. fasterstart[[1]]; reallyfaststart=NDSolve[{pendulumoscillator,y[0] == 0.8,y'[0]==2.9},y[t],{t, 0, endtime}];yfaster[t_] = 
      y[t] /. fasterstart[[1]]; yreallyfast[t_] = 
      y[t] /. reallyfaststart[[1]]; Plot[{yslow[t], yfaster[t], yreallyfast[t]}, {t, 0, endtime}, 
      PlotStyle 
        -> {{Red, Thickness[0.01]}, {CadmiumOrange, Thickness[0.01]}, {
            HotPink, Thickness[0.01]}}, 
      AxesLabel -> {"t", "y[t]"}, 
      PlotRange -> {-3, 10}, AspectRatio -> (1/2)];

One of these looks just like an ordinary oscillator.
Another looks something like an ordinary oscillator.
The third doesn't look like an ordinary oscillator at all.
Use linearizations to help explain why this happened.

The flow field for the nonlinear pendulum oscillator is:

\(Clear[m, n, x, y, t, Field]; 
m[x_, y_] = \(-\(\(g\ Sin[y]\)\/L\)\); 

    \(n[x_, y_] = x; \)\)Field[x_, y_] = {m[x, y], n[x, y]}

The linearization matrix at {0, 0} is:

{xeq, yeq} = {0, 0}; \(Clear[A, gradm, gradn]; 

    gradm[x_, y_] = {∂\_x m[x, y], ∂\_y m[x, y]}; 

    \(gradn[x_, y_] = {∂\_x n[x, y], ∂\_y n[x, y]}; 
    \)\)A[x_, y_] = {gradm[x, y], gradn[x, y]}; MatrixForm[A[xeq, yeq]]

The linearization at {0, 0} is:

\(linearized\  = \ 

      \(({\(x'\)[t], \(y'\)[t]}\  == \ 
          Expand[A[xeq, yeq] . {x[t]\  - \ xeq, y[t]\  - \ yeq}])\); 

    ColumnForm[Thread[linearized]]\)

The flow field for the linearization is:

\(Clear[linearizedField]; 

    linearizedField[x_, y_] = A[xeq, yeq] . {x - xeq, y - yeq}\)

See the linearized flow and the nonlinear pendulum flow:

h = 3.5; {xlow, xhigh} = {xeq - h, xeq + h}; {ylow, yhigh} = {yeq - h, yeq + h}; jump = ((xhigh - xlow)/12); scalefactor = 0.35; sizer = 0.3; nonlinearflowplot = 
      Table[Arrow[Field[x, y], Tail -> {x, y}, VectorColor -> Blue, 
          ScaleFactor -> scalefactor, HeadSize -> sizer], {x, xlow, 
          xhigh, jump}, {y, ylow, yhigh, jump}]; linearizedflowplot = 
      Table[Arrow[linearizedField[x, y], Tail -> {x, y}, 
          VectorColor -> Red, ScaleFactor -> scalefactor, 
          HeadSize -> sizer], {x, xlow, xhigh, jump}, {y, ylow, yhigh, 
          jump}]; equilibriumplot = 
      Graphics[{CadmiumOrange, PointSize[0.05], Point[{xeq, yeq}]}]; both = 
      Show[nonlinearflowplot, equilibriumplot, linearizedflowplot, 
        PlotRange 
          -> {{xeq - 1.5 h, xeq + 1.5 h}, {yeq - 1.5 h, 
              yeq + 1.5 h}}, Axes -> True, 
        AxesLabel -> {"y'", "y"}, PlotLabel -> "Both"];

Near {0, 0}, the two flows are similar - almost copies of each other.
The upshot:
When the starter data on y[0] and y'[0]are close to {0, 0}, the nonlinear pendulum oscillator behaves very much like the ordinary linear oscillator that comes from its linearization.
Watch it happen for random starter data close to {0, 0}:

h = 0.5; \(Clear[y, ypend, yordinary, t]; 
    
{startery, starteryprime} = {Random[Real, {\(-h\), h}], 
        Random[Real, {\(-h\), h}]}\)pendulumoscillator=y''[t]+((g Sin[y[t]])/L)==0;linearizedoscillator=y''[t]+((g y[t])/L)==0;g = 9.8; L = 5; endtime = 15; pendulum=NDSolve[{pendulumoscillator,y[0] == startery,y'[0]==starteryprime},y[t],{t, 0, endtime}];ypend[t_] = y[t] /. pendulum[[1]]; 
    ordinary=NDSolve[{linearizedoscillator,y[0] == startery,y'[0]==starteryprime},y[t],{t, 0, endtime}];yordinary[t_] = 
      y[t] /. ordinary[[1]]; Plot[{ypend[t], yordinary[t]}, {t, 0, endtime}, 
      PlotStyle -> {{Red, Thickness[0.012]}, {Blue, Thickness[0.008]}}, 
      AxesLabel -> {"t", "y[t]"}, 
      PlotLabel -> "Starter data close to {0,0}", 
      AspectRatio -> (1/3), PlotLegend -> {"pend", "lin"}, 
      LegendSize -> 0.4];

Take another look at both flows near {0, 0}:

Show[both];

The plot indicates that when you go with starter data on {y[0],y'[0]}not so close to {0, 0}, the pendulum oscillator behaves much differently from the ordinary oscillator coming from its linearization. See what happens when you go with some random starter data not so close to {0, 0}:

\(Clear[y, ypend, yordinary, t]; 
    
{startery, starteryprime} = {Random[Real, {\(-1.0\), 1.0}], 
        Random[Real, {1.5, 4.0}]}\)pendulumoscillator=y''[t]+((g Sin[y[t]])/L)==0;linearizedoscillator=y''[t]+((g y[t])/L)==0;g = 9.8; L = 5; endtime = 15; pendulum=NDSolve[{pendulumoscillator,y[0] == startery,y'[0]==starteryprime},y[t],{t, 0, endtime}];ypend[t_] = y[t] /. pendulum[[1]]; 
    ordinary=NDSolve[{linearizedoscillator,y[0] == startery,y'[0]==starteryprime},y[t],{t, 0, endtime}];yordinary[t_] = 
      y[t] /. ordinary[[1]]; Plot[{ypend[t], yordinary[t]}, {t, 0, endtime}, 
      PlotStyle -> {{Red, Thickness[0.012]}, {Blue, Thickness[0.008]}}, 
      AxesLabel -> {"t", "y[t]"}, 
      PlotLabel -> "Starter velocity not close to 0", 
      AspectRatio -> (1/3), PlotLegend -> {"pend", "lin"}, 
      LegendSize -> 0.4];

B.5) Linearizations and gradients:
        Getting an idea about why linearizations give good approximations

B.5.a)

The linearization of a function f[x] at a point (Subscript[x, 0]) on the x - axis is:

Clear[f, x, linearf, x0];
linearf[x_]=f[x0]+f'[x0] (x - x0)

See the plot of a function f[x] together with the plot of its linearization at a random point (Subscript[x, 0]) on the x - axis:

f[x_] = 0.2 x + Sin[2 x]; {xlow, xhigh} = {-1, 4}; \(x0 = Random[Real, {xlow, xhigh}]; 


    Plot[{f[x], linearf[x]}, {x, xlow, xhigh}, 
      PlotStyle -> {{Blue, Thickness[0.01]}, {Red, Thickness[0.01]}}, 
      AxesLabel -> {"x", ""}, PlotRange -> {\(-2\), 2.5}, 
      PlotLabel -> 
        "Function and linearization 
 at indicated point", 
      PlotLegend -> {"f[x]", "linf[x]"}, LegendSize -> 0.6, 
      Epilog -> {Red, PointSize[0.04], Point[{x0, f[x0]}]}]; \)

That's right! 
The linearization of f[x] at  (Subscript[x, 0]) is nothing but the tangent      
       y=f[(Subscript[x, 0])]+f'[(Subscript[x, 0])] (x - (Subscript[x, 0]))
at  (Subscript[x, 0]).
It gives a really good approximation of f[x] near (Subscript[x, 0]), the point of linearization.
How do you use gradients to carry this off for functions of two variables?

One of the great drivers of mathematics is the power of analogy.
To linearize a function f[x] of one variable at (Subscript[x, 0]), you go with the derivative 
       f'[x]
and write 
      linf[x]=f[(Subscript[x, 0])]+f'[(Subscript[x, 0])] (x - (Subscript[x, 0])).
         
To linearize a function f[x, y] of two variables at {(Subscript[x, 0]), (Subscript[y, 0])}, you go with the gradient 
       ∇f[x, y] = {∂f[x, y]/∂x, 
          ∂f[x, y]/∂y}

which is the same as
       ∇f[x, y] = {(∂_x) f[x, y], (∂_y) f[x, y]}
and write
      linf[x, y] = 
        f[(Subscript[x, 0]), (Subscript[y, 0])] + ∇f[(Subscript[x, 0]), (Subscript[y, 0])] . {x - (Subscript[y, 0]), y - (Subscript[y, 0])}
which is the same as
      linf[x, y] = 
        f[(Subscript[x, 0]), (Subscript[y, 0])] + gradf[(Subscript[x, 0]), (Subscript[y, 0])] . {x - (Subscript[y, 0]), y - (Subscript[y, 0])}

To help you stick your mental tongs into the idea behind using the gradient for linearization, look at the following plot of a sample function f[x, y]:

\(Clear[f, x, y]; 
f[x_, y_] = 3 - 0.4\ x\^2 - 0.5\ y\^2; 

    \(actualsurfaceplot = 
      ParametricPlot3D[{x, y, f[x, y]}, {x, \(-1\), 1}, {y, \(-1\), 1}, 
        Axes -> Automatic, 
        AxesLabel -> {"x", "y", "f[x,y]"}, 
        PlotRange -> All, ViewPoint -> CMView, 
        DisplayFunction -> $DisplayFunction]; \)\)

Now look at plots of linearizations at the indicated points:

{x0, y0} = {Random[Real, {0, 1}], Random[Real, {-1, 1}]}; Clear[gradf, linearf, point]; gradf[x_, y_] = {(∂_x) f[x, y], (∂_y) f[x, y]}; 
    linearf[x_, y_] = f[x0, y0] + gradf[x0, y0] . {x - x0, y - y0}; 
    h = 0.4; linearizationplot = 
      ParametricPlot3D[{x, y, linearf[x, y]}, {x, x0 - h, x0 + h}, {y, 
          y0 - h, y0 + h}, PlotPoints -> {2, 2}, 
        DisplayFunction -> Identity]; linpoint = 
      Graphics3D[{Red, PointSize[0.02], Point[{x0, y0, f[x0, y0]}]}]; Show[actualsurfaceplot, linearizationplot, linpoint];

That flat sheet hugging onto the actual surface plot is a plot of the linearization at the indicated point. 
In fact, that sheet is tangent to the surface and approximates the actual surface very well close to the indicated point.

Some calculous courses call the flat sheet by the name "tangent plane."

B.5.b)

Here's a nonlinear diffeq system:

\(Clear[m, n, x, y]; 
m[x_, y_]\  = \ 0.4\ x\  - \ 0.1\ x\ y; 

    n[x_, y_]\  = \ \(-0.3\)\ y\  + \ 0.5\ y\ x\  - \ 0.5; 

    DiffEqsystem\  = \ 
      \(({\(x'\)[t], \(y'\)[t]}\  == \ {m[x[t], y[t]], n[x[t], y[t]]})\); 

    ColumnForm[Thread[DiffEqsystem]]\)

The equilibrium points are:

equisols = Solve[{m[x, y] == 0, n[x, y] == 0}, {x, y}]; \(Clear[xeq, yeq]; 
    
{xeq[1], yeq[1]} = {x, y} /. 
        equisols[[1]]\){xeq[2], yeq[2]} = {x, y} /. 
        equisols[[2]]

The linearized system at {xeq[2], yeq[2]} is:

\(Clear[A, gradm, gradn]; 

    gradm[x_, y_]\  = \ {D[m[x, y], x], D[m[x, y], y]}; 

    gradn[x_, y_]\  = \ {D[n[x, y], x], D[n[x, y], y]}; 

    A[x_, y_]\  = \ {gradm[x, y], gradn[x, y]}; 


    linearized2\  = \ 

      \(({\(x'\)[t], \(y'\)[t]}\  == \ 
          Chop[Expand[
              A[xeq[2], yeq[2]] . {x[t]\  - \ xeq[2], y[t]\  - \ yeq[2]}]])
        \); 
ColumnForm[Thread[linearized2]]\)

How is this linearized system related to the linearizations of m[x, y] and n[x, y] at {xeq[2], yeq[2]}?
How does the answer explain why the linearized system approximates the original system very well near {xeq[2], yeq[2]}?

Very intimately.
To get the idea, linearize m[x, y] and linearize n[x, y] at {xeq[2], yeq[2]}:

Clear[linearm, linearn]; linearm[x_, y_] = 
      Chop[m[xeq[2], yeq[2]] + 
          gradm[xeq[2], yeq[2]] . {x - xeq[2], y - yeq[2]}]\(
\)linearn[x_, y_] = 
      Chop[n[xeq[2], yeq[2]] + 
          gradn[xeq[2], yeq[2]] . {x - xeq[2], y - yeq[2]}]

Use linearm[x, y] and linearn[x, y] to make a new diffeq system:

\(newDiffEqsystem\  = \ 

      \(({\(x'\)[t], \(y'\)[t]}\  == \ 
          Expand[{linearm[x[t], y[t]], linearn[x[t], y[t]]}])\); 

    ColumnForm[Thread[newDiffEqsystem]]\)

This is the same as the linearized system you got via the linearization matrix:

ColumnForm[Thread[linearized2]]

It works this way everytime. 
And that's why linearized systems approximate the original system very, very well at the point of linearization.

B.5.c)

In what other mathematical issues is the idea of linearization important?

The idea of linearization sits at the heart of the chain rule for functions of more than one variable.
And the idea of linearization sits at the heart of the way you change variables in double and triple integrals.
It's around all over the place.