DE.08  Linearizations
TUTORIALS

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]
]

T.1) Using formulas resulting from the linearization of a pendulum oscillator to estimate the amplitude and frequency of a pendulum oscillator

T.1.a.i) Approximate formulas: The good news

The truth is that nothing, including Mathematica, has ever succeeded in coming up with an exact formula for the actual pendulum oscillator. Try it in the sample case L = 7.5 and g = 9.8 with starter data y[0] = 0.4 and y'[0]=-0.3:

L = 7.5; g = 9.8; starterangle = 0.4; startervelocity = -0.3; Clear[t, y];
pendulumoscillator=y''[t]+((g Sin[y[t]])/L)==0

actual=DSolve[{pendulumoscillator,y[0] == starterangle,y'[0]==startervelocity},y[t],t]

Out of luck. Formulas for the solution do not exist.
Lots of formula-hungry folks turn to the linearized version for a quick formula:

linpendulumoscillator=y''[t]+((g y[t])/L)==0Clear[approxformulay];
linearized=Chop[DSolve[{linpendulumoscillator,y[0] == starterangle,y'[0]==startervelocity},y[t],t]];approxformulay[t_] = 
      N[ComplexExpand[
          y[t] /. linearized[[1]]], 6]

How do you know whether or not to trust this approximate formula?

Look at the starter data:

{starterangle, startervelocity}

The starter point is not far from the point of linearization at   {0, 0}; so you can expect this formula to give good results. Give it the acid test:

endtime = 12; realthing=NDSolve[{pendulumoscillator,y[0] == starterangle,y'[0]==startervelocity},y[t],{t, 0, endtime}];yreal[t_] = y[t] /. realthing[[1]]; 
    Plot[{yreal[t], approxformulay[t]}, {t, 0, endtime}, 
      PlotStyle -> {{Thickness[0.015], Red}, {Thickness[0.008], Blue}}, 
      AxesLabel -> {"t", ""}, 
      AspectRatio -> (1/GoldenRatio)];

Nailed it fairly well. You could probably use this formula for t between 0 and 12.
See more:

endtime = 50; realthing=NDSolve[{pendulumoscillator,y[0] == starterangle,y'[0]==startervelocity},y[t],{t, 0, endtime}];yreal[t_] = y[t] /. realthing[[1]]; 
    Plot[{yreal[t], approxformulay[t]}, {t, 0, endtime}, 
      PlotStyle -> {{Thickness[0.015], Red}, {Thickness[0.008], Blue}}, 
      AxesLabel -> {"t", ""}, 
      AspectRatio -> (1/(2 GoldenRatio))];

Whether to use the formula for a long time interval is a question of your own judgment.

T.1.a.ii)

Look at the approximate formula for the pendulum oscillator you got in part i) above:

approxformulay[t]

Use it to estimate the amplitude, period and frequency of this pendulum oscillator.

Look at the approximate formula again:

approxformulay[t]

Write down the amplitude:

approxamplitude = (Sqrt[(0.4^2) + (((-0.262445))^2)])

Confirm with a plot:

amplitudeplot = 
      Plot[{yreal[t], approxamplitude, -approxamplitude}, {t, 0, 
          endtime}, 
        PlotStyle 
          -> {{Thickness[0.01], Red}, {Thickness[0.008], Blue}, {
              Thickness[0.008], Blue}}, AxesLabel -> {"t", ""}, 
        AspectRatio -> (1/(2 GoldenRatio))];

Estimate the period of the pendulum oscillator:

approxformulay[t]

approxperiod = Solve[1.1431 t == N[2 π], t]

This pendulum oscillator repeats itself approximately every 5.5 time units.

approxperiod = 5.5; \(periodlines = 
      Table[Graphics[{UltramarineViolet, Thickness[0.01], 
            Line[{{t, \(-approxamplitude\)}, {t, approxamplitude}}]}], {t, 0, 
          8\ approxperiod, approxperiod}]; 

    \(Show[amplitudeplot, periodlines]; \)\)

Lots of folks like to talk about frequency. The approximate frequency of this pendulum oscillator is:

appoxfrequency = (1/approxperiod)

This oscillator goes through about  0.18 complete oscillations per unit of time.

T.1.a.iii) Approximate formulas: The bad news

Can you always use the approximate formula you get from the linearization to estimate amplitude, period, and frequency of the pendulum oscillator?

For crying out loud, this sounds like a question that wuss Calculus Cal would ask (if he had the brains to come up with it).  
The quality of estimates you get from the linearization depends on the quality of the approximation resulting from the linearization.  When you blindly trust the stuff you get from linearizations, you are running the risk of walking in deep do-do.
Look at this one:

L = 7.5; g = 9.8; Clear[t, y];
pendulumoscillator=y''[t]+((g Sin[y[t]])/L)==0starterangle = 0.5; startervelocity = -2.0; simplependulumoscillator=y''[t]+((g y[t])/L)==0Clear[approxformulay]; linearized=Chop[DSolve[{simplependulumoscillator,y[0] == starterangle,y'[0]==startervelocity},y[t],t]];approxformulay[t_] = 
      N[y[t] /. linearized[[1]], 7]

endtime = 12; realthing=NDSolve[{pendulumoscillator,y[0] == starterangle,y'[0]==startervelocity},y[t],{t, 0, endtime}];yreal[t_] = y[t] /. realthing[[1]]; 
    Plot[{yreal[t], approxformulay[t]}, {t, 0, endtime}, 
      PlotStyle -> {{Thickness[0.015], Red}, {Thickness[0.008], Blue}}, 
      AxesLabel -> {"t", ""}, 
      AspectRatio -> (1/GoldenRatio)];

The "approximate formula" resulting from the linearization here is:

approxformulay[t]

Any one who would try to use this approximate formula to estimate anything about the actual pendulum oscillator should look for the clue phone or run the risk of getting the reputation of a total dork.  

In fact the only DiffEq&Mathematica student who might be tempted  to use this formula to estimate anything is Calculus Cal himself.

Cal, get yourself and your brain-rot out of here right now.

T.2) Why most folks don't care about  linearizing at points other than equilibrium points

T.2.a)

Given a nonlinear system, folks calculate the equilibrium points and then check out the eigenvalues of the linearization matrix to determine the nature of the equilibrium points.
Why don't folks do this at points other than equilibrium points?

It doesn't give you useful information.
To see why, go with 
        x'[t]=m[x[t], y[t]]
        y'[t]=n[x[t], y[t]].
If{xnoteq, ynoteq} is not an equilibrium point, then the flow at {xnoteq, ynoteq} is given by
       {m[xnoteq, ynoteq], n[xnoteq, ynoteq]},
And because {xnoteq, ynoteq} is not an equilibrium point, this flow is not {0, 0}.
So:
The flow at {xnoteq, ynoteq} goes through {xnoteq, ynoteq} in the direction of the vector:
       {m[xnoteq, y noteq], n[xnoteq, y noteq]}.
As a result, 
→{xnoteq, ynoteq} cannot be an attractor. 
→{xnoteq, ynoteq} cannot be a repellor. 

Check it out in this example:
The system:

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

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

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

    ColumnForm[Thread[diffeqsystem]]\)

The equilbrium points:

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 flow with the equilibrium points:

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

    \({xlow, xhigh} = {\(-3\), 3}; \)\){ylow, yhigh} = {-3, 3}; jump = ((xhigh - xlow)/12); flowplot = 
      Table[Arrow[Field[x, y], Tail -> {x, y}, VectorColor -> Blue, 
          ScaleFactor -> 0.3, HeadSize -> 0.3], {x, xlow, xhigh, 
          jump}, {y, ylow, yhigh, jump}]; equilibriumplots = 
      Table[Graphics[{Red, PointSize[0.04], Point[{xeq[k], yeq[k]}]}], {k, 1, 
          2}]; setup = 
      Show[flowplot, equilibriumplots, Axes -> True, 
        AxesLabel -> {"x", "y"}];

Away from the equilibrium points, the flow is just like "Old Man River."  The flow just keeps on rolling along. 
Throw in a random point that is not an equilibrium point:

{xnoteq, ynoteq} = {Random[Real, {-2, 2}], 
        Random[Real, {-1, 1}]}notequilibriumplot = 
      Graphics[{Magenta, PointSize[0.04], Point[{xnoteq, ynoteq}]}]; Show[setup, notequilibriumplot];

The flow goes right through the new point.

T.3) Using Lyapunov's rules to investigate the predator- prey model

T.3.a.i)

The predator-prey model is
           x'[t]=  a x[t] -  b x[t] y[t], 
           y'[t]=-c y[t] +  d y[t] x[t].
 with a, b, c and d all positive.

When you feed extra prey in at a rate of r units per time unit the model becomes
          x'[t]=   a x[t] -  b x[t] y[t] +  r
          y'[t]= -c y[t] +  d y[t] x[t].

Go with the sample case of
         a  =  0.7,  b  =  0.3,  c  =  0.44  and  d  =  0.08
and calculate the equilibrium points as functions of r. 
Discuss the results.

Set it up with cleared a, b, c, d and r:

\(Clear[m, n, x, y, a, b, c, d, r, t]; 

    m[x_, y_]\  = \ a\ x\  - \ b\ x\ y\  + \ r; 

    n[x_, y_]\  = \ \(-c\)\ y\  + \ d\ y\ x; 

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

    ColumnForm[Thread[DiffEqsystem]]\)

The equilibrium points are:

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

The only equilbrium point of any physical interest is:

\(Clear[xeq, yeq, r]; 
    
{xeq[r_], yeq[r_]} = {c\/d, \(a\ c + d\ r\)\/\(b\ c\)}\)

Now incorporate the given values of a, b, c and d to get the equilibrium points as functions of r:

a = 0.7; b = 0.3; c = 0.44; d = 0.08; {xeq[r], yeq[r]}

This indicates that the equilibrim prey population is 5.5  REGARDLESS of what r is. 
The equilibrium predator population varies with r this way:

Plot[yeq[r], {r, -5, 12}, 
      PlotStyle -> {{Red, Thickness[0.01]}}, 
      AxesLabel -> {"r", "predator equilibrium"}];

If you feed in prey at a rate of r prey units per time unit and you make r more negative than:

Solve[yeq[r] == 0, r]

Then the equilibrium predator population becomes negative.
This is not good.

T.3.a.ii)

Keep everything the same as in part i) and keep everything in part i) live in your computer's memory.
The equilibrium points as functions of r are:

{xeq[r], yeq[r]}

Use Lyapunov's rules to determine the nature of the equilibrium points.

The linearization (Jacobian) 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]]

Calculate the eigenvalues of the linearization matrix (Jacobian) at {xeq[r], yeq[r]}:

{eigenvalue1[r_], eigenvalue2[r_]} = 
      Chop[Eigenvalues[A[xeq[r], yeq[r]]]]

Check a few out:

ColumnForm[
      Table[{r, eigenvalue1[r], eigenvalue2[r]}, {r, -2, 16, 2}]]

All but the large r's.
Look at the oh-so important real number that determines whether   
         {xeq[r], yeq[r]} 
is an attractor or a repellor:

ColumnForm[
      Chop[Table[{r, (eigenvalue1[r] + Conjugate[eigenvalue1[r]])/2, 
            (eigenvalue2[r] + Conjugate[eigenvalue2[r]])/2}, {r, -2, 
            16, 2}]]]

Plot as a function of r:

eigenvalplot = 
      Plot[{(1/2) (eigenvalue1[r] + Conjugate[eigenvalue1[r]]), 
          (1/2) (eigenvalue2[r] + Conjugate[eigenvalue2[r]])}, {r, 
          -2, 20}, 
        PlotStyle -> {{Red, Thickness[0.015]}, {Blue, Thickness[0.01]}}, 
        AspectRatio -> (1/GoldenRatio), 
        AxesLabel -> {"r", ""}];

The plot signals that:
-> The equillibrium points {xeq[r], yeq[r]} are attractors for r > 0.
-> The equillibrium points{xeq[r], yeq[r]} are repellors for  r < 0.
Here's a sample for positive r:

\(\(r\  = \ 2.5; 

h = 2; 

    
{xstarter, ystarter}\  = 
      \ {Random[Real, {xeq[r]\  - \ h, xeq[r]\  + \ h}], 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        Random[Real, {yeq[r]\  - \ h, yeq[r]\  + \ 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[r], yeq[r]}]}]; 

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


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

    \ \ \ \ \ \ \ \ \ \ \ \ 

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

    \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\)

A genuine attractor.
See what happens for a sample negative r:

\(\(r\  = \ \(-0.3\); 

h\  = \ 2; 
    
{xstarter, ystarter}\  = 
      \ {Random[Real, {xeq[r]\  - \ h, xeq[r]\  + \ h}], 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        Random[Real, {yeq[r]\  - \ h, yeq[r]\  + \ 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[r], yeq[r]}]}]; 

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


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

    \ \ \ \ \ \ \ \ \ \ \ \ 

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

    \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\)

Looks like a genuine repellor.

T.3.a.iii)

Look at the plot of the real parts of the eigenvalues of the linearization matrix:

Show[eigenvalplot];

How do you interpret that fork?

Take another look:

Show[eigenvalplot];

This bifurcation happens for r somewhere between 12 and 13:
to see what the fork means, look at:

r = 12.62; Eigenvalues[A[xeq[r], yeq[r]]]

r = 12.63; Eigenvalues[A[xeq[r], yeq[r]]]

r = 12.64; Eigenvalues[A[xeq[r], yeq[r]]]

The bifurcation point r between swirl and no swirl is between 12.63 and 12.64.
Take another look at the fork:

Show[eigenvalplot];

The handle of the fork on the left corresponds to (theoretical) swirl. 
The right hand tines of the fork correspond to no swirl at all.