DE.02 Transition from Calculus to DiffEq:
The Forced Oscillator DiffEq
       y''[t]+ b y'[t]+ c y[t]=f[t]
 LITERACY

Needs["Graphics`Colors`"];
Graphics`Colors`GosiaGreen=RGBColor[0, 0.392187, 0];
Needs["Graphics`FilledPlot`"];
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]
]

L.1)

Here are plots of two unforced oscillators:

b = 0.4; c = 5.1; Clear[y, t];
unforcedoscillator=y''[t]+b y'[t]+c y[t]==0;
initialdisplacement = 1.5; initialvel = 3.0; endtime = 20; Clear[s, Derivative, y, notforcedy, t];
notforcedsolution=NDSolve[{unforcedoscillator,y[0] == initialdisplacement,y'[0]==initialvel},y[t],{t, 0, endtime}];
notforcedy[t_] = 
        y[t] /. notforcedsolution[[1]]; 
      notforcedplot1 = 
      Plot[notforcedy[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Blue}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}, AspectRatio -> (1/3)]; 
    b = 0.2; c = 5.1; Clear[y, t];
unforcedoscillator=y''[t]+b y'[t]+c y[t]==0;
initialdisplacement = 1.5; initialvel = 3.0; Clear[s, Derivative, y, notforcedy, t];
notforcedsolution=NDSolve[{unforcedoscillator,y[0] == initialdisplacement,y'[0]==initialvel},y[t],{t, 0, endtime}];
notforcedy[t_] = 
        y[t] /. notforcedsolution[[1]]; 
      notforcedplot2 = 
      Plot[notforcedy[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Blue}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}, AspectRatio -> (1/3)];

One of these plots is the solution of
       y''[t] + 0.2 y'[t] + 5.1 y[t]=0
with 
        y[0] = 1.5 and y'[0]=3.0.

The other is a plot of the solution of
        y''[t] + 0.4 y'[t] + 5.1 y[t]=0;
with 
        y[0] = 1.5 and y'[0]=3.0.

Which is which?

L.2)

Here are plots of the solutions of
         diffeq1):   y''[t] + 4.7 y[t]=0
          diffeq2):  y''[t] + 0.3 y'[t] + 4.7 y[t]=0,
          diffeq3):  y''[t] + 1.3 y'[t] + 4.7 y[t]=0,
          diffeq4):  y''[t] + 7.3 y'[t] + 4.7 y[t]=0,
all with the same starter data
          y[0] = 2 and y'[0]=3.
The plots are not in order.  Your job is to match the differential equation with the plot of its solution.

Clear[y, t, y1, y2, y3, y4];
sol1=DSolve[{y''[t]+0.3 y'[t]+4.7 y[t]==0,y[0] == 2,y'[0]==3},y[t],t];
y1[t_] = 
        y[t] /. sol1[[1]]; sol2=DSolve[{y''[t]+1.3 y'[t]+4.7 y[t]==0,y[0] == 2,y'[0]==3},y[t],t];y2[t_] = y[t] /. sol2[[1]]; sol3=DSolve[{y''[t]+7.3 y'[t]+4.7 y[t]==0,y[0] == 2,y'[0]==3},y[t],t];y3[t_] = y[t] /. sol3[[1]]; sol4=DSolve[{y''[t]+4.7 y[t]==0,y[0] == 2,y'[0]==3},y[t],t];y4[t_] = Chop[y[t] /. sol4[[1]]]; 
    endtime = 8; Plot[y3[t], {t, 0, endtime}, PlotStyle -> {{Thickness[0.01], Blue}}, 
      AxesLabel -> {"t", ""}, AspectRatio -> (2/endtime), 
      PlotRange -> {-2.5, 2.5}, PlotLabel -> "Plot A"]; Plot[y2[t], {t, 0, endtime}, PlotStyle -> {{Thickness[0.01], Blue}}, 
      AxesLabel -> {"t", ""}, AspectRatio -> (2/endtime), 
      PlotRange -> {-2.5, 2.5}, PlotLabel -> "Plot B"]; Plot[y4[t], {t, 0, endtime}, PlotStyle -> {{Thickness[0.01], Blue}}, 
      AxesLabel -> {"t", ""}, AspectRatio -> (2/endtime), 
      PlotRange -> {-2.5, 2.5}, PlotLabel -> "Plot C"]; Plot[y1[t], {t, 0, endtime}, PlotStyle -> {{Thickness[0.01], Blue}}, 
      AxesLabel -> {"t", ""}, AspectRatio -> (2/endtime), 
      PlotRange -> {-2.5, 2.5}, PlotLabel -> "Plot D"];

diffeq1) ⟶ Plot........   diffeq2) ⟶ Plot........
           diffeq3) ⟶ Plot........   diffeq4) ⟶ Plot........
  How does the starter data signal that each plot must go up before it can go down?

L.3)

a) Write down the  characteristic equation for this unforced linear oscillator diffeq
       y''[t] + 6 y'[t] + 25 y[t]=0
b)  Solve the characteristic equation.
c)  What information about all solutions of 
       y''[t] + 6 y'[t] + 25 y[t]=0
do you get  by inspecting the solutions of the characteristic equation?

L.4)

a) Write down the  characteristic equation for this unforced linear oscillator diffeq
       y''[t] + 7 y'[t] + 12 y[t]=0
b) Solve the characteristic equation.
c) Explain how the solutions of the characteritic equation signal take all solutions of
       y''[t] + 7 y'[t] + 12 y[t]=0 
head to 0 without oscillation.

L.5)

a)  Write down the characteristic equation for this unforced linear oscillator diffeq
       y''[t] + 7 y'[t] + 25 y[t]=0.
b) What information about all solutions of 
       y''[t] + 7 y'[t] + 25 y[t]=0
do you get when you solve the charcteristic equation?

L.6)

Here are three plots of solutions of a damped forced oscillator 
       y''[t] + 0.7 y'[t] + 5.5 y[t]=2 Sin[3 t]  +  Cos[1.5 t] .
Two of the solutions have random starting values on y[0] and y'[0].
The other solution is yzeroinput[t] which has starting values y[0] = 0 and y'[0]=1:

b  =  0.7;
c  =  5.5;
f[t_]  =  2 Sin[3 t]  +  Cos[1.5 t];
forcedoscillator  =  
        y''[t]  +  b y'[t]  +  c y[t]  ==  f[t];
 
endtime  =  30;
Clear[y, y1, y2, y3, t];
initialdisplacement1  =  Random[Real, {4, 10}];
initialvel1  =  Random[Real, {-6, 6}];
initialdisplacement2  =  0;
initialvel2  =  0;
initialdisplacement3  =  Random[Real, {-10, -4}];
initialvel3  =  Random[Real, {-6, 6}];

\(sol1\  = \ 
        NDSolve[

          \ \ \ {forcedoscillator, 
\ \ \ \ 
            y[0]\  == \ initialdisplacement1, \(y'\)[0]\  == \ initialvel1}, 
          y[t], {t, 0, endtime}]\);
\(sol2\  = \ 
        NDSolve[

          \ \ \ {forcedoscillator, 
\ \ \ \ 
            y[0]\  == \ initialdisplacement2, \(y'\)[0]\  == \ initialvel2}, 
          y[t], {t, 0, endtime}]\);
    
\(sol3\  = \ 
        NDSolve[

          \ \ \ {forcedoscillator, 
\ \ \ \ 
            y[0]\  == \ initialdisplacement3, \(y'\)[0]\  == \ initialvel3}, 
          y[t], {t, 0, endtime}]\);
    
    
    
y1[t_]  =  y[t] /. sol1[[1]];
y2[t_]  =  y[t] /. sol2[[1]];
y3[t_]  =  y[t] /. sol3[[1]];
\(plots\  = \ 
        Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
\ \ \ \ \ \ 
          PlotStyle -> {{Thickness[0.01], Blue}, {Thickness[0.01], Red}, 

              \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
                Thickness[0.01], Black}}, 
\ \ \ \ \ \ PlotRange -> All, 
          AspectRatio -> Automatic, 
\ \ \ \ \ \ 
          AxesLabel -> {"t", ""}]\);

The three solutions begin their trip in totally different ways, but when t is large, you need a scorecard to tell them apart. Explain how you could have anticipated this in advance.

L.7)

Here are three plots of solutions of the UNDAMPED forced oscillator 
       y''[t] + 5.5 y[t]=2 Sin[3 t]  +  Cos[1.5 t] .
Two of the solutions have random starting values on y[0] and y'[0].
The other solution is yzeroinput[t] which has starting values y[0] = 0 and y'[0]=1:

b  =  0;
c  =  5.5;
f[t_]  =  2 Sin[3 t]  +  Cos[1.5 t];
forcedoscillator  =  
        y''[t]  +  b y'[t]  +  c y[t]  ==  f[t];
 
endtime  =  30;
Clear[y, y1, y2, y3, t];
initialdisplacement1  =  Random[Real, {4, 10}];
initialvel1  =  Random[Real, {-6, 6}];
initialdisplacement2  =  0;
initialvel2  =  0;
initialdisplacement3  =  Random[Real, {-10, -4}];
initialvel3  =  Random[Real, {-6, 6}];

\(sol1\  = \ 
        NDSolve[

          \ \ \ {forcedoscillator, 
\ \ \ \ 
            y[0]\  == \ initialdisplacement1, \(y'\)[0]\  == \ initialvel1}, 
          y[t], {t, 0, endtime}]\);
\(sol2\  = \ 
        NDSolve[

          \ \ \ {forcedoscillator, 
\ \ \ \ 
            y[0]\  == \ initialdisplacement2, \(y'\)[0]\  == \ initialvel2}, 
          y[t], {t, 0, endtime}]\);
    
\(sol3\  = \ 
        NDSolve[

          \ \ \ {forcedoscillator, 
\ \ \ \ 
            y[0]\  == \ initialdisplacement3, \(y'\)[0]\  == \ initialvel3}, 
          y[t], {t, 0, endtime}]\);
    
    
    
y1[t_]  =  y[t] /. sol1[[1]];
y2[t_]  =  y[t] /. sol2[[1]];
y3[t_]  =  y[t] /. sol3[[1]];
\(plots\  = \ 
        Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
\ \ \ \ \ \ 
          PlotStyle -> {{Thickness[0.01], Blue}, {Thickness[0.01], Red}, 

              \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
                Thickness[0.01], Black}}, 
\ \ \ \ \ \ PlotRange -> All, 
          AspectRatio -> Automatic, 
\ \ \ \ \ \ 
          AxesLabel -> {"t", ""}]\);

This time the solutions do not settle into a common steady state. 
Explain how you could have anticipated this in advance.

L.8)

Which do you expect to decay to 0 faster:
a)  Solutions of y''[t] + 0.3 y'[t] + 8.7 y[t]=0.
b)  Solutions of y''[t] + 0.9 y'[t] + 8.7 y[t]=0.
My response is.....................

L.9)

Many folks would say that
a) The oscillator  y''[t] + 4 y[t]=0 is undamped.
b) The oscillator  y''[t] + 2 y'[t] + 4 y[t]=0 is under damped.
c) The oscillator  y''[t] + 4 y'[t] + 4 y[t]=0 is critically damped
d) The oscillator  y''[t] + 4.2 y'[t] + 4 y[t]=0 is overdamped.
Where are these folks coming from?

L.10)

Here are plots of forced oscillators coming from 
        y''[t]+0.3 y'[t]+1.6 y[t]=f[t]
with
        y[0] = 4.0 and y'[0]=-0.7
for five choices of forcing functions f[t]:

b = 0.3; c = 1.6; initialdisplacement = 4.0; endtime = 40; initialvel = -0.7; Clear[f, t, forcedy];
f[t_] = 8 Sin[t];
forcedoscillator=y''[t]+b y'[t]+c y[t]==f[t];forcedsolution=NDSolve[{forcedoscillator,y[0] == initialdisplacement,y'[0]==initialvel},y[t],{t, 0, endtime}];forcedy[t_] = 
      y[t] /. forcedsolution[[1]]; forcedplot = 
      Plot[forcedy[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Plot 1", AspectRatio -> (1/3)]; \(Clear[f, t]; 
f[t_] = 12.3\ DiracDelta[t - 20]; 
\(b = 0.3; \)\)c = 1.6; \(Clear[K1, K2, z, z1, z2, generalsol, t]; 
zeqn = z\^2 + b\ z + c == 0; 

    \(zsols = Solve[zeqn, z]; \)\)generalsol[t_] = 
      K1 (E^(z1 t)) + K2 (E^(z2 t)) /. {
          z1 -> zsols[[1, 1, 2]], 
          z2 -> zsols[[2, 1, 2]]}; 
    Ksols=Solve[{generalsol[0] == initialdisplacement,generalsol'[0]==initialvel}];Clear[transient];
transient[t_] = 
        Chop[ComplexExpand[
            generalsol[t] /. Ksols[[1]]]]
          ;
Ksols=Solve[{generalsol[0] == 0,generalsol'[0]==1}];\(Clear[green]; 

    green[t_] = 
      Chop[ComplexExpand[
          generalsol[t] /. Ksols[[1]]]]; 
    
Clear[steadystate, x]; 

    steadystate[t_] = 
      Simplify[Chop[
          ∫\_0\%t\( f[x]\ green[t - x]\) \[DifferentialD]x]]; 

    Clear[formulay]; 
formulay[t_] = Expand[transient[t] + steadystate[t]]; 
    
\(oscplot = 
      Plot[formulay[t], {t, 0, endtime}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Plot 2", AspectRatio -> 1\/3]; \)\)b = 0.3; c = 1.6; initialdisplacement = 4.0; initialvel = -0.7; Clear[f, t, forcedy];
f[t_] = 2.6;
forcedoscillator=y''[t]+b y'[t]+c y[t]==f[t];forcedsolution=NDSolve[{forcedoscillator,y[0] == initialdisplacement,y'[0]==initialvel},y[t],{t, 0, endtime}];forcedy[t_] = 
      y[t] /. forcedsolution[[1]]; forcedplot = 
      Plot[forcedy[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Plot 3", AspectRatio -> (1/3)]; b = 0.3; c = 1.6; initialdisplacement = 4.0; initialvel = -0.7; Clear[f, t, forcedy];
f[t_] = 2.3 (E^(0.05 t)) Cos[t];
forcedoscillator=y''[t]+b y'[t]+c y[t]==f[t];forcedsolution=NDSolve[{forcedoscillator,y[0] == initialdisplacement,y'[0]==initialvel},y[t],{t, 0, endtime}];forcedy[t_] = 
      y[t] /. forcedsolution[[1]]; forcedplot = 
      Plot[forcedy[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Plot 4", AspectRatio -> (1/3)]; b = 0.3; c = 1.6; initialdisplacement = 4.0; initialvel = -0.7; Clear[f, t, forcedy];
f[t_] = 7.3 (E^(-0.4 t));
forcedoscillator=y''[t]+b y'[t]+c y[t]==f[t];forcedsolution=NDSolve[{forcedoscillator,y[0] == initialdisplacement,y'[0]==initialvel},y[t],{t, 0, endtime}];forcedy[t_] = 
      y[t] /. forcedsolution[[1]]; forcedplot = 
      Plot[forcedy[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Plot 5", AspectRatio -> (1/3)];

The five forcing functions f[t] used in these plots are: 
         f[t] = 7.3 (E^(-0.4 t))
         f[t] = 12.3 DiracDelta[t - 20]
         f[t] = 2.3 (E^(0.05 t)) Cos[t]
         f[t] = 8 Sin[t]
 and
         f[t] = 2.6
 Your job is to match the forcing functions to the plots.

L.11)

a) When you go after a formula for the oscillator coming from a forced, damped oscillator
               y''[t]+b y'[t]+c y[t]=f[t]
               (with b > 0 and c > 0)
with given starter data
               y[0] = p and y'[0]=q,
you put 
             yformula[t] = yunforced[t]  +  yzeroinput[t]
where yunforced[t] solves
             y''[t] + b y'[t] + c y[t]=0
with 
            y[0]=p and y'[0]=q
and yzeroinput[t] solves
            y''[t] + b y'[t] + c y[t]=f[t]
with 
             y[0] = 0 and y'[0]=0.

Explain these statements:
a) You are guaranteed that yunforced[t] -> 0 as t -> ∞.
(This is why some folks call yunforced[t] by the name "transient.")

b) If you are interested in looking at the steady state (long term global scale) behavior of this oscillator, then all you have to do is to plot yzeroinput[t].

c) Does the steady state (long term, global scale) behavior of this oscillator depend in any way on the specific values of the starting data? 
How do you know for sure? 

d) How do your answers change if you drop the damping term and go after a formula for the oscillator coming from a forced, undamped oscillator
               y''[t]+c y[t]=f[t]  (with c > 0)
with given starter data
               y[0] = p and y'[0]=q?

L.12)

a) Explain the Mathematica output:

ComplexExpand[(E^((-0.26 + I 3.78) t))]

b) Explain the Mathematica output:

ComplexExpand[(E^((-0.26 - I 3.78) t))]

c) Explain the Mathematica output:

TrigToExp[Sin[((π t)/4)]]

d) Explain the Mathematica output:

TrigToExp[Cos[((π t)/2)]]

L.13)

When you use the characteristic equation to get the form of the general solution of
          y''[t] + 25 y[t]=0, 
you get:

\(Clear[t]; 
b = 0; 
\(c = 25; \)\)\(Clear[K1, K2, z, z1, z2, generalsol, t]; 

    charequation = z\^2 + b\ z + c == 0; 
\(zsols = Solve[charequation, z]; 
    \)\)generalsol[t_] = 
      K1 (E^(z1 t)) + K2 (E^(z2 t)) /. {
          z1 -> zsols[[1, 1, 2]], 
          z2 -> zsols[[2, 1, 2]]}

But when you look in the old printed text used in many DiffEq classes, you read that the general solution is
           C1 Cos[5 t]  +  C2 Sin[k t].
How do you reconcile these two answers?

L.14)

a) If y1[t] solves
        y''[t] + 3 y'[t] + 9 y[t]=0
with 
        y[0] = 4 and y'[0]=-2, 
and if y2[t] solves
         y''[t] + 3 y'[t] + 9 y[t]=0.2 Sin[t]
with 
        y[0] = 0 and y'[0]=0, 
then what differential equation does
        y[t] = y1[t]  +  y2[t]
solve?  Include starter data.

b)  Somebody hands you a function y1[t] solving
         y''[t] + 4 y[t]=DiracDelta[t - 6]
with 
         y[0] = 1 and y'[0]=-1, 

Explain how you would go about coming up with a the formula for function y2[t] such that
         y[t] = y1[t]  +  y2[t]
solves 
         y''[t] + 4 y[t]=DiracDelta[t - 6]
with 
         y[0] = 1 and y'[0]=1.
      
c) If y1[t] solves
        y''[t] + 3.2 y'[t] + 9.8 y[t]=Cos[8 t]
with 
        y[0] = 4 and y'[0]=-2, 
and if y2[t] solves
         y''[t] + 3.2 y'[t] + 9.8 y[t]=0.2 Sin[t]
with 
        y[0] = 3 and y'[0]=3, 
then what differential equation does
        y[t] = y1[t]  +  y2[t]
solve?  Include starter data.

L.15)

a) Given an oscillator differential equation 
          y''[t] + b y'[t] + c y[t]=0,
what information do you get from the solutions of the characteristic equation?

b) If b > 0 and (b^2) - 4 c > 0, why do most folks say that the oscillator coming from
            y''[t] + b y'[t] + c y[t]=0,
is overdamped?

L.16)

Here's how you use Mathematica and the convolution integral method to go after an exact formula for the solution of
        y''[t] + 4 y'[t] + 13 y[t]=3.6 Sin[t]
with
         y[0] = 3 and y'[0]=-1:

\(Clear[f, t]; 
f[t_] = 3.6\ Sin[t]; 
\(b = 4; \)\)c = 13; \(Clear[K1, K2, z, z1, z2, generalsol, t]; 

    charequation = z\^2 + b\ z + c == 0; 
\(zsols = Solve[charequation, z]; 
    \)\)generalsol[t_] = 
      K1 (E^(z1 t)) + K2 (E^(z2 t)) /. {
          z1 -> zsols[[1, 1, 2]], 
          z2 -> zsols[[2, 1, 2]]}; 
    Ksols=Solve[{generalsol[0] == 3,generalsol'[0]==-1}];\(Clear[yunforced]; 

    yunforced[t_] = 
      Chop[ComplexExpand[
          generalsol[t] /. Ksols[[1]]]]
        \)

What differential equation does yunforced[t] solve?
Include starter data.

b) The next step is to calculate the unit impulse response:

Ksols=Solve[{generalsol[0] == 0,generalsol'[0]==1}];\(Clear[yinitimpulse]; 

    yunitimpulse[t_] = 
      Chop[ComplexExpand[
          generalsol[t] /. Ksols[[1]]]]
        \)

And then you convolve the unit impulse response function with f[t] to get the zero input solution:

\(Clear[yzeroinput, x]; 

    yzeroinput[t_] = 
      Expand[Chop[
          ∫\_0\%t\( f[x]\ yunitimpulse[t - x]\) \[DifferentialD]x]]
        \)

What differential equation does yzeroinput[t] solve?
Include starter data.

c) The last step is to put 
          yformula[t] = yunforced[t]  +  yzeroinput[t].

\(Clear[yformula]; 

    yformula[t_]\  = \ yunforced[t]\  + \ yzeroinput[t]\)

What differential equation does yformula[t] solve?
Include starter data.
Identify the transient and the steady state parts of the solution.

L.17)

What do folks mean when they talk about resonance?
What do folks mean when they talk about beating?

L.18)

a) Go with {a, b} different from {0, 0} and write down the choice of number k that puts the oscillator coming from
            y''[t] + 9 y[t]=a Cos[k t]  +  b Sin[k t]
into resonance.

b) Go with {a, b} different from {0, 0} and write down a choice of number k that makes the oscillator coming from
           y''[t] + 9 y[t]=a Cos[k t]  +  b Sin[k t]
beat without being in true resonance.

L.19)

a) You are given a solution y[t] of 
       y''[t]+2 y'[t]+5 y[t]=6  DiracDelta[t - 4.7].
What happens to y'[t] at t = 4.7?
b)  How does your response indicate that
       DiracDelta[t  -  4.7] ≠ 2  DiracDelta[t  -  4.7]?

L.20)

a) Go with any positive h and write down the value of  
       (∫_( b - h)^(b + h ))12 DiracDelta[t - b] 
          dt

b) Write down the values of
       (∫_0^4.99999) DiracDelta[t - 5] dt
       (∫_0^5.00001) DiracDelta[t - 5] dt
       (∫_5.00001^(100 ))  
        DiracDelta[t - 5] dt.
      
c) Here's a plot of UnitStep[t  -  2]:

\(Clear[t]; 

    Plot[UnitStep[t - 2], {t, 0, 6}, 
      PlotStyle -> {{Blue, Thickness[0.015]}}, 
      AxesLabel -> {"t", ""}, AspectRatio -> 1\/3]; \)

And here's a plot of UnitStep[t  -  4]:

\(Clear[t]; 

    Plot[UnitStep[t - 4], {t, 0, 6}, 
      PlotStyle -> {{Blue, Thickness[0.015]}}, 
      AxesLabel -> {"t", ""}, AspectRatio -> 1\/3]; \)

Digest the two plots and explain the result of this Mathematica calculation:

\(Clear[t, x]; 

    ∫\_0\%t\( Sin[x]\ DiracDelta[x - 5]\) \[DifferentialD]x\)

d) Calculate this convolution integral by hand:
           (∫_0^t) DiracDelta[x - 5] g[t - x] dx
for 
           g[t] = (E^(-t)) Cos[3 t].

L.21)

a) Just to see that you can do something with no machine help,
use the convolution integral method and hand calculation to come up with the exact formula for the solution of
         y''[t] + 5 y'[t] + 6 y[t]=3.8 (E^(-t))
with
         y[0] = 2 and y'[0]=-1.

b) Just to see that you can do even more with no machine help,
use the convolution integral method and hand calculation to come up with the exact formula for the solution of this impulse-forced oscillator:
        y''[t] + 2 y'[t] + 5 y[t]=7 DiracDelta[t - 4]
with 
        y[0] = 3 and y'[0]=0.