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

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

T.1)   Amplitude and frequency of unforced oscillators

T.1.a) Amplitude and frequency of unforced damped oscillators:

You are in the lab and are in the process of coming up with a formula for this given damped oscillator with these given starter values on y[0] and y'[0]:

\(Clear[y, t]; 
b = 0.8; 
\(c = 24.0; \)\)linoscdiffeq=y''[t]+b y'[t]+c y[t]==0startery = 5.9; starteryprime = -0.4;

You copy, paste and edit from the Basics:

Clear[z];
charequation = (z^2) + b z + c == 0;
zsols = Solve[charequation, z];

Clear[gensol, C1, C2];
gensol[t_] = 
        C1 (E^(zsols[[1, 1, 2]] t)) + 
          C2 (E^(zsols[[2, 1, 2]] t
                ));

ystarteq = gensol[0] == startery;
yprimestarteq=gensol'[0]==starteryprime;\(Csols = Solve[{ystarteq, yprimestarteq}, {C1, C2}]; 

Clear[yformula]; 
    
yformula[t_] = 
      Chop[ComplexExpand[
          gensol[t] /. Csols[[1]]]]\)

Your job is to calculate the amplitude and the frequency for this damped unforced oscillator. Do it and show off your work with convincing plots.

Answer:

Look at:

yformula[t]

Read off:

\(Clear[amplitude, t]; 

    amplitude[t_] = \@\(5.9\^2 + 0.401424\^2\)\ E\^\(\(-0.4\)\ t\)\)

frequency = N[(4.88262/(2 π))]

This damped oscillator goes through about  0.8 oscillations per time unit.
See 1 complete oscillation:

endtime = (1/frequency); oscplot = 
      Plot[{yformula[t], amplitude[t], -amplitude[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], Red}, {
              Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "y[t]"}, PlotPoints -> 100, 
        AspectRatio -> (1/2)];

See 2 complete oscillations:

endtime = (2/frequency); oscplot = 
      Plot[{yformula[t], amplitude[t], -amplitude[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], Red}, {
              Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "y[t]"}, PlotPoints -> 100, 
        AspectRatio -> (1/2)];

See 4 complete oscillations:

endtime = (4/frequency); oscplot = 
      Plot[{yformula[t], amplitude[t], -amplitude[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], Red}, {
              Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "y[t]"}, PlotPoints -> 100, 
        AspectRatio -> (1/2)];

See 8 complete oscillations:

endtime = (8/frequency); oscplot = 
      Plot[{yformula[t], amplitude[t], -amplitude[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], Red}, {
              Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "y[t]"}, PlotPoints -> 100, 
        AspectRatio -> (1/2)];

Petering out.
Throw in a movie:

endtime = (8/frequency); \(Clear[oscplotter, s]; 

    oscplotter[s_] := 
      Plot[{yformula[t], amplitude[t], \(-amplitude[t]\)}, {t, 0, 
          s\/frequency}, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], Red}, {
              Thickness[0.01], Red}}, 
        PlotRange 
          -> {{0, endtime}, {\(-amplitude[0]\) - 0.1, 
              amplitude[0] + 0.1}}, AxesLabel -> {"t", "y[t]"}, 
        PlotPoints -> 100, 
        PlotLabel -> s\ "= complete oscillations", 
        AspectRatio -> 1\/3]; 
\(Table[oscplotter[s], {s, 1, 8}]; 
    \)\)

T.1.b) Amplitude and frequency of unforced undamped oscillators

Will the same thing work for undamped unforced oscillators?

Yah sure! You bettcha!
Copy, paste, edit and put it to work on this one:

\(Clear[y, t]; 
b = 0; 
\(c = 32.7; \)\)linoscdiffeq=y''[t]+b y'[t]+c y[t]==0startery = 12.7; starteryprime = 3.1;

\(Clear[z]; 
charequation = z\^2 + b\ z + c == 0; 

    zsols = Solve[charequation, z]; 

Clear[gensol, C1, C2]; 

    gensol[t_] = 
      C1\ E\^\(zsols[[1, 1, 2]]\ t\) + 
        C2\ E\^\(zsols[[2, 1, 2]]\ t\); 
    
\(ystarteq = gensol[0] == startery; \)\)yprimestarteq=gensol'[0]==starteryprime;\(Csols = Solve[{ystarteq, yprimestarteq}, {C1, C2}]; 

Clear[yformula]; 
    
yformula[t_] = 
      Chop[ComplexExpand[
          gensol[t] /. Csols[[1]]]]\)

Read off:

\(Clear[amplitude, t]; 
amplitude[t_] = \@\(12.7\^2 + 0.54211\^2\)\)

frequency = N[(5.71839/(2 π))]

This oscillator goes through about almost 1 complete oscillation in one time unit.

See 1 complete oscillation:

endtime = (1/frequency); oscplot = 
      Plot[{yformula[t], amplitude[t], -amplitude[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], Red}, {
              Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "y[t]"}, PlotPoints -> 100, 
        AspectRatio -> (1/2)];

See 2 complete oscillations:

endtime = (2/frequency); oscplot = 
      Plot[{yformula[t], amplitude[t], -amplitude[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], Red}, {
              Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "y[t]"}, PlotPoints -> 100, 
        AspectRatio -> (1/2)];

See 4 complete oscillations:

endtime = (4/frequency); oscplot = 
      Plot[{yformula[t], amplitude[t], -amplitude[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], Red}, {
              Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "y[t]"}, PlotPoints -> 100, 
        AspectRatio -> (1/2)];

See 8 complete oscillations:

endtime = (8/frequency); oscplot = 
      Plot[{yformula[t], amplitude[t], -amplitude[t]}, {t, 0, endtime}, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], Red}, {
              Thickness[0.01], Red}}, 
        AxesLabel -> {"t", "y[t]"}, PlotPoints -> 100, 
        AspectRatio -> (1/2)];

Throw in a movie:

endtime = (8/frequency); \(Clear[oscplotter, s]; 

    oscplotter[s_] := 
      Plot[{yformula[t], amplitude[t], \(-amplitude[t]\)}, {t, 0, 
          s\/frequency}, 
        PlotStyle 
          -> {{Thickness[0.015], Blue}, {Thickness[0.01], Red}, {
              Thickness[0.01], Red}}, 
        PlotRange 
          -> {{0, endtime}, {\(-amplitude[0]\) - 0.1, 
              amplitude[0] + 0.1}}, AxesLabel -> {"t", "y[t]"}, 
        PlotPoints -> 100, 
        PlotLabel -> s\ "= complete oscillations", 
        AspectRatio -> 1\/3]; 
\(Table[oscplotter[s], {s, 1, 8}]; 
    \)\)

T.2)  Impulse-forced oscillators: The physical meaning of the impulse  hit

T.2.a.i)

Showcase the convolution integral method and the Dirac delta function by coming up with the exact formula for the forced oscillator coming from 
       y''[t] + 0.9 y'[t] + 4.1 y[t]=f[t]
with  f[t] = 10 DiracDelta[t - 7],  y[0] = 2.0  and  y'[0]=7.1.
Show off your work with a good plot.
Describe what you see.

Step 1: Calculate yunforced[t]

\(Clear[f, t]; 
f[t_] = 10\ DiracDelta[t - 7]; 
\(b = 0.9; \)\)c = 4.1; ystarter = 2.0; \(yprimestarter = 7.1; 

Clear[C1, C2, z, z1, z2, generalsol, t]; 

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

Csols=Solve[{generalsol[0] == ystarter,generalsol'[0]==yprimestarter}];\(Clear[yunforced]; 

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

Step 2: Calculate the unit impulse response, yunitimpulse[t]

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

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

Step 3: Calculate 
           yzeroinput[t] = 
        (∫_0^t) yunitimpulse[t - x] f[x] dx

\(Clear[yzeroinput, x]; 

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

Step 4: Set formulay[t] = yzeroinput[t] + yunforced[t]:

\(Clear[formulay]; 

    formulay[t_] = Expand[yzeroinput[t] + yunforced[t]]\)

And a plot:

\(\(oscplot = 
      Plot[formulay[t], {t, 0, 12}, 
		
        PlotStyle -> {{Blue, Thickness[0.01]}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}]; \)\)

Just when this oscillator seemed to be settling down, the impulse hit at t = 7 goosed this oscillator.
For the fun of it, look at this plot of yformula[t] and yzeroinput[t]:

Plot[{formulay[t], yzeroinput[t]}, {t, 0, 12}, 
      PlotStyle -> {{Blue, Thickness[0.02]}, {Red, Thickness[0.01]}}, 
      PlotRange -> All, AxesLabel -> {"t", "y[t]"}];

Think about what this plot is trying to tell you.

T.2.a.ii) Impulses and Dirac Delta Functions

Take another peek at the plot of the impulse-forced oscillator coming from
       y''[t] + 0.9 y'[t] + 4.1 y[t]=f[t]
with  f[t] = 10 DiracDelta[t - 7],  y[0] = 2.0  and  y'[0]=7.1.
Show off your work with a good plot.
Describe what you see.

tail = {4, 4}; Show[oscplot, 
      Arrow[{7, formulay[7]} - tail, Tail -> tail, 
        VectorColor -> Magenta]];

Explain why the oscillator had a sudden change of direction at t = 7.

Look at the differential equation
       y''[t] + 0.9 y'[t] + 4.1 y[t]=10 DiracDelta[t - 7]
and remember that 
       DiracDelta[t - 7] = 0  for t not equal to 7. 
So, aside from the what happens at t = 7, this oscillator behaves just as an oscillator coming from
       y''[t] + 0.9 y'[t] + 4.1 y[t]=0.
Look again:

tail = {6, 4}; Show[oscplot, 
      Arrow[{7, formulay[7]} - tail, Tail -> tail, 
        VectorColor -> Magenta]];

The differential equation
       y''[t] + 0.9 y'[t] + 4.1 y[t]=10 DiracDelta[t - 7]
tells you that the forcing term 
       10 DiracDelta[t - 7] 
packs all its punch at t = 7 and abruptly redirects the oscillator from its original path to a new (unforced) path. 
This is in direct response to the impulse force that 
             10 DiracDelta[t - 7] 
lays on the oscillator at t = 7.

Show[oscplot, 
      Arrow[{7, formulay[7]} - tail, Tail -> tail, 
        VectorColor -> Magenta], 
      Graphics[Text["impulse hit", tail]]];

T.2.a.iii) Physical meaning of the impulse hit

Stay with the same setup:
       y''[t] + 0.9 y'[t] + 4.1 y[t]=f[t]
with  f[t] = 10 DiracDelta[t - 7],  y[0] = 2.0  and  y'[0]=7.1.
and look at this plot of the derivative of the solution:

velocityplot=Plot[formulay'[t],{t, 0, 12},PlotStyle -> {{Red, Thickness[0.01]}},AxesLabel -> {"t", "velocity[t]"},DisplayFunction -> Identity];

tail = {4, 4};
label = Graphics[{Text["Impulse hit", tail]}]; pointer=Arrow[{6.999,formulay'[6.999]}-tail,Tail -> tail];

setup = 
        Show[velocityplot, label, pointer, 
          DisplayFunction -> $DisplayFunction];

Use this plot to interpret the physical meaning of the impulse hit with f[t] = 10 DiracDelta[t - 7] at t = 7.

The impulse hit with
       f[t]=10 DiracDelta[t - 7]
instantaneously jumps the velocity up by 10 at the instant of the impulse hit at t = 7.

T.2.a.iv)

Go with d > 0, and go with
       y''[t] + b y'[t] + c y[t]=K DiracDelta[t  -  d]
with given starting values on y[0] and y'[0].
The experiment above sends the message that the impulse hit with
       f[t] = K DiracDelta[t - d]
instantaneously jumps  
       y'[t] up by K
at the instant of the impulse hit at t = d.
Explain why this is guaranteed.

Go with h > 0 and integrate both sides of 
       y''[t] +  b y'[t] +  c y[t]=K DiracDelta[t -  d]
from d - h to d + h to get
       (∫_(d - h)^(d + h))(y''[t]+ b y'[t]+ c y[t])dt 
        = 
        K (∫_(d - h)^(d + h  ))DiracDelta[t - d] 
              dt
        = K

This is the same as
       (∫_(d - h )^(d + h))  y''[t]dt + (∫_(d - h)^(d + h)) b y'[t]dt + (∫_(d - h)^(d + h))c y[t] 
              dt=K
As h closes in on 0,
       (∫_(d - h)^(d + h))c y[t] dt  
        →  0
and
       (∫_(d - h)^(d + h))b y'[t]dt →0.

Consequently, as h closes in on 0, 
       (∫_(d - h)^(d + h))y''[t]dt→K.
But 
       (∫_(d - h)^(d + h))y''[t]dt=y'[d + h]-y'[d - h]
So as h closes in on 0,
       y'[d + h]-y'[d - h]→ K.

The upshot: 
y'[t] jumps by K at t = d.

T.3) Convolutions involving DiracDelta[x - d] are easy to do by hand

T.3.a.i)

Look at this experiment involving a convolution integral:

\(Clear[f, g, x, t]; 
d = 2.5; 
\(f[x_] = Cos[x]; 

    g[t_] = ∫\_0\%t\( f[t - x]\ DiracDelta[x - d]\) 
          \[DifferentialD]x\)\)

And plot:

\(\(gplot\  = \ 
      Plot[g[t], {t, 0, 8}, 
		PlotStyle -> {{Thickness[0.015], Blue}}, 

        \ \ \ \ \ \ \ \ \ AspectRatio -> 1\/2, 
        AxesLabel -> {"t", "g[t]"}, PlotRange -> All, 

        \ \ \ \ \ \ \ \ \ Epilog -> {Red, Text["d", {d, \(-0.08\)}]}]; 
    \)\)

Explain what happened by explaining why
       (∫_0^t) f[t - x] DiracDelta[x - d] dx = 
        f[t - d] UnitStep[t - d].

Because DiracDelta[x - d] = 0 for x < d, you are guaranteed that

       g[t] = 
        (∫_0^t) f[t - x] DiracDelta[x - d] 
              dx = 0 for t < d.
  
 But for t > d,
       g[t] = 
        (∫_0^t) f[t - x] DiracDelta[x - d] 
              dx = f[t - d].
    
 So g[t] = f[t - d] UnitStep[t - d].

See the plot of  f[t - d] UnitStep[t - d]   :

thisplot = 
      Plot[f[t - d] UnitStep[t - d], {t, 0, 8}, 
        PlotStyle -> {{Thickness[0.01], Red}}, AspectRatio -> (1/2), 
        AxesLabel -> {"t", "g[t]"}, 
        Epilog -> {Red, Text["d", {d, -0.08}]}];

Compare with the original:

Show[gplot, thisplot];

Just as theory predicted, they are the same!

T.4)  Ramp-forced oscillators

T.4.a)

Here is a ramp forcing function and its plot:

\(Clear[ramp, g, t]; 
g[t_] = 0.5\ \((E\^\(0.7\ t\) - 1)\); 

    \(changeover = 2.7; \)\)ramp[t_] = 
      g[t] UnitStep[changeover - t] + 
        g[changeover] UnitStep[t - changeover]Plot[ramp[t], {t, 0, 9}, PlotStyle -> {{Thickness[0.01], Red}}, 
      PlotRange -> All, AxesLabel -> {"t", "ramp[t]"}, 
      AspectRatio -> (1/GoldenRatio)];

Ramp[t] runs 
→ with g[t] = 0.5 ((E^(0.7 t)) - 1) for t < 2.7
→ with g[2.7] for t > 2.7.
Here is a the plot of a damped oscillator forced by ramp[t]:

\(Clear[y, t, f]; 
f[t_] = ramp[t]; 
\(b = 1.6; \)\)c = 4.2; ystarter = 3.5; yprimestarter = 1.5; diffeq=y''[t]+b y'[t]+c y[t]==f[t];

endtime = 10;
solution=NDSolve[{diffeq,y[0] == ystarter,y'[0]==yprimestarter},y[t],{t, 0, endtime}];

approxy[t_] = 
        y[t] /. solution[[1]];
ndsplot = 
        Plot[approxy[t], {t, 0, endtime}, 
          PlotStyle -> {{Thickness[0.015], Blue}}, 
          PlotRange -> All, AspectRatio -> (1/GoldenRatio), 
          AxesLabel -> {"t", "y[t]"}];

Use the convolution integral method to produce an exact formula for this ramp-forced oscillator.

Here you go:

Step 1: Calculate yunforced[t]

\(Clear[C1, C2, z, z1, z2, generalsol, t]; 

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

Csols=Solve[{generalsol[0] == ystarter,generalsol'[0]==yprimestarter}];\(Clear[yunforced]; 

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

Step 2: Calculate the unit impulse response, yunitimpulse[t]

Csols=Solve[{generalsol[0] == 0,generalsol'[0]==1}];\(Clear[yunitimpulse]; 

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

Step 3: Calculate 
           yzeroinput[t] = 
        (∫_0^t) unitimpulse[t - x] f[x] dx

\(Clear[yzeroinput, x]; 

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

Step 4: Set formulay[t] = yzeroinput[t] + yunforced[t]:

\(Clear[formulay]; 

    formulay[t_] = Simplify[yzeroinput[t] + yunforced[t]]\)

Here is the formula plot together with the NDSolve plot:

formulaplot = 
      Plot[formulay[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.008], Red}}, 
        DisplayFunction -> Identity]; Show[ndsplot, formulaplot];

Perfecto.

T.5)  Electrical folks and the forced oscillator

T.5.a)

Lots of folks other than spring specialists are interested in the oscillator.  In fact, electrical folks really get into the oscillator in a big way. Why?

When spring specialists think of the forced oscillator differential equation
       y''[t] + b y'[t] + c y[t]=f[t], 
they think of:
-> t as time
-> y[t] as displacement from mass stretched position,
-> b as a drag measurement,
-> c as a measurement of the stiffness of the spring, and
-> f[t] as an externally applied force which varies with time. 

When the electrical specialist thinks of the the oscillator differential equation
       y''[t] + b y'[t] + c y[t]=f[t], 
they think of:
-> t as time
-> y[t] as voltage
-> b as a measurement of resistance,
-> c as a measurement of impedance and
-> f[t] as a measurement of electromotive force.
Don't worry about these words if they are unfamiliar to you.