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

\(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

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] == 2Perpend[a_?VectorQ] := 
      Normalize[
          If[a[[1]] == 0, {1, 0, 0}, {0, 
              a[[3]], 
              -a[[2]]}]] /; 
        Length[a] == 3Base3D = 
      Table[{Cos[((2 k π)/8.) ], Sin[((2 k π)/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] == 2Vector[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[#1, tip]]&)/@
                    Partition[
                      (tip - hr diff + #1&)/@
                        ((ap hr Norm[diff] Base3D) . {temp, 
                              Normalize[temp\[Cross]diff]}), 2, 1]]}}]] /; 
        Length[from] == Length[to] == 3Arrow[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[]; \)\)

VectorField and Phase Lines

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, #1 ≥ a && #1 ≤ b&]; 
        zero1 = Sort[zero ⋃ {a, b}]; len = Length[zero]; 
        len1 = Length[zero1]; 
        signs = Table[
            Sign[fcn /. 
                var -> 
                  (1/2) (
                      zero1[[k + 1]] + 
                        zero1[[k]])], {
              k, 1, len1 - 1}]; 
        stick = {Thickness[0.007], Line[{{pos, a}, {pos, b}}]}; 
        dots = (Point[{pos, #1}]&)/@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 slen)/4)]; 
        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, #1 ≥ a && #1 ≤ b&]; 
        zero1 = Sort[zero ⋃ {a, b}]; len = Length[zero]; 
        len1 = Length[zero1]; 
        signs = Table[
            Sign[fcn /. 
                var -> 
                  (1/2) (
                      zero1[[k + 1]] + 
                        zero1[[k]])], {
              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]], #1}]&)/@
            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), (slen/2)]; 
        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, #1]&)/@Options[picture]; 
          opt3 = opt3 //. {
                (AxesLabel -> _) -> AxesLabel -> None, 
                (PlotRange -> _) -> 
                  PlotRange -> Automatic, 
                (PlotRegion -> _) -> 
                  PlotRegion -> Automatic, 
                (Ticks -> _) -> Ticks -> Automatic}; 
          picture = List@@picture; 
          picture[[-1]] = opt3; 
          picture = 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[];

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, π\/8, \(3\ π\)\/8}, 
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[0.02]}}, 
        Axes -> None, AspectRatio -> 1\/3, 
        DisplayFunction -> Identity]; \)

B.1)  Underdamped oscillators (b > 0 and small) and overdamped oscillators (b > 0 and large)

B.1.a.i)

The undamped oscillator diffeq is y''[t]+ c y[t]=0  with c > 0.
Here's a random undamped oscillator diffeq with random starter values for y[0] and y'[0]:

\(c\  = \ Random[Real, \ {3, \ 5}]; \ 

    startery\  = \ Random[Real, \ {2, \ 5}]; \ 

    starteryprime\  = \ Random[Real, \ {\(-2\), \ 2}]; \ 
Clear[t, \ y]; 

    undampedoscdiffeq\  = \ 

      \ \ \ {\(\(y'\)'\)[t]\  + \ c\ y[t]\  == \ 0, \ y[0]\  == \ startery, \ 
        
\ \ \ \ \(y'\)[0]\  == \ starteryprime}; \ 

    ColumnForm[Thread[undampedoscdiffeq]]\)

The solution plots out this way:

\(endtime = 30; 

    ndsol = NDSolve[undampedoscdiffeq, y[t], {t, 0, endtime}]; 

    Clear[yundamped]; 

    yundamped[t_] = y[t] /. ndsol[[1]]; 
    
undampedplot = 
      Plot[yundamped[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Blue}}, 
        AspectRatio -> 1\/4, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}]; \)

A sine wave.
Regular repeating oscillations forever.
How do clued-in diffeq folks reach in and damp this oscillator so that the the solutions tend to 0 as t gets large?

The undamped oscillator diffeq is
            y''[t] + c y[t]=0:

\(Clear[t, \ y]; 

    undampedoscdiffeq\  = \ 

      \ \ \ {\(\(y'\)'\)[t]\  + \ c\ y[t]\  == \ 0, \ y[0]\  == \ startery, \ 
        
\ \ \ \ \ \(y'\)[0]\  == \ starteryprime}; \ 

    ColumnForm[undampedoscdiffeq]\)

To damp this oscillator so that all the solutions tend to 0 as t gets large, you keep everything the same but you insert a damping term b y'[t] ( with b > 0) into the middle of the diffeq to get the damped oscillator diffeq
       y''[t] + b y'[t] + c y[t]=0:

\(b\  = \ 0.18; 
Clear[t, \ y]; 

    dampedoscdiffeq\  = \ 

      \ \ \ {\(\(y'\)'\)[t]\  + \ b\ \(y'\)[t]\  + \ c\ y[t]\  == \ 0, \ 
        y[0]\  == \ startery, \ 
\(y'\)[0]\  == \ starteryprime}; \ 

    ColumnForm[dampedoscdiffeq]\)

The solution plots out this way:

\(endtime = 30; 

    ndsol = NDSolve[dampedoscdiffeq, y[t], {t, 0, endtime}]; 


    Clear[ydamped]; 

    ydamped[t_] = y[t] /. ndsol[[1]]; 
    dampedplot1 = 
      Plot[ydamped[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, AspectRatio -> 1\/4, 
        PlotRange -> All, AxesLabel -> {"t", "y[t]"}]; \)

A damped sine wave which is squashed to 0 as t gets large.
See the undamped solution and the damped solution together:

Show[undampedplot, dampedplot1];

Damped and undamped versions of the same oscillator.

B.1.a.ii)

Here is the plot of the damped oscillator from part i)

Show[dampedplot1];

How do you damp it even more?

To get this plot, you inserted b y'[t] with this b:

b

To damp it even more you raise b:

higherb  =  0.5;
Clear[t, y];
moredampedoscdiffeq = {y''[t] + higherb y'[t] + c y[t] == 0,
             y[0]  ==  startery,y'[0]  ==  starteryprime};
ColumnForm[Thread[moredampedoscdiffeq]]

The new solution plots out this way:

endtime  =  30;
ndsol = NDSolve[
   moredampedoscdiffeq  ,y[t],{t, 0, endtime}];

Clear[newydamped];
ymoredamped[t_] = 
        y[t] /. ndsol[[1]];dampedplot2 = 
        Plot[ymoredamped[t], {t, 0, endtime}, 
          PlotStyle -> {{Thickness[0.01], Magenta}}, 
          AspectRatio -> (1/4), PlotRange -> All, 
          AxesLabel -> {"t", "y[t]"}];

See them all:

Show[undampedplot, dampedplot1, dampedplot2];

The higher b is, the more damping you get.

B.1.a.iii) Underdamped versus overdamped

Here's what happens when you increase the damping term b quite a bit more:

bigb  =  5.2;
Clear[t,  y];
\(overdampedoscdiffeq\  = \ 

          \ \ \ {\(\(y'\)'\)[t]\  + \ bigb\ \(y'\)[t]\  + \ c\ y[t]\  == \ 0, 
            \ y[0]\  == \ startery, \ 
\(y'\)[0]\  == \ starteryprime}\); 
ColumnForm[overdampedoscdiffeq]

The new solution plots out this way:

\(endtime = 30; 

    ndsol = NDSolve[overdampedoscdiffeq, y[t], {t, 0, endtime}]; 

    Clear[yoverdamped]; 

    yoverdamped[t_] = 
      y[t] /. ndsol[[1]]; 

    overdampedplot = 
      Plot[yoverdamped[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.015], Red}}, 
        AspectRatio -> 1\/4, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Overdamped"]; \)

Why do folks call this by the name overdamped?

They say this ocscillator is overdamped because the plot doesn't wiggle.
Too much friction.

B.2)  Forcing the linear oscillator:
                 yforced[t] = yunforced[t] + yzeroinput[t] 
        For damped linear oscillators:
                 yunforced[t] is transient and yzeroinput[t] gives the steady state behavior

B.2.a.i)

Look at this plot of the solution of:
        y''[t] + 0.4 y'[t] + 2.1 y[t]=0;
        with  y[0] = 1.5  and  y'[0]=3.0:

b = 0.4; c = 2.1; Clear[y, t];
unforcedoscillator=y''[t]+b y'[t]+c y[t]==0;
initialdisplacement = 1.5; initialvel = 3.0; endtime = 30; Clear[s, y, unforcedy, t];
unforcedsolution=NDSolve[{unforcedoscillator,y[0] == initialdisplacement,y'[0]==initialvel},y[t],{t, 0, endtime}];
\(unforcedy[t_] = 
        y[t] /. unforcedsolution[[1]]; 

      
unforcedplot = 
        Plot[unforcedy[t], {t, 0, endtime}, 
          PlotStyle -> {{Thickness[0.01], MediumVioletRed}}, 
          PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
          PlotLabel -> "Unforced Oscillator", 
          AspectRatio -> 1\/2]; \)

Here's what you get when you take the same damped oscillator and the same starter data, but at time t, you throw in an external force proportional to 
        f[t] = 2 Sin[t]:

Clear[f, t, forcedy];
f[t_] = 2 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 -> "Forced Oscillator", AspectRatio -> 1\/2]; 
    \)

Compare:

both = 
      Show[unforcedplot, forcedplot, PlotLabel -> "Both", 
        AspectRatio -> (1/4)];

Describe what you see and try to explain why you see it.

Answer:

The unforced (damped) linear oscillator comes from:

unforcedoscillator

with starter data y[0] = 1.5 and y'[0]=3.0;
The forced (damped) linear oscillator comes from:

forcedoscillator

with the same starter data.
Take another look:

Show[both];

Evidently the effect of the forcing function 
      f[t] = 2 Sin[t] 
eventually overcomes the damping term and sends the damped oscillator into a regular steady state.

B.2.a.ii) Steady state and the zero input response

Look at a new forced linear oscillator:

b = Random[Real, {0.2, 0.8}]; c = Random[Real, {2.0, 3.5}]; \(Clear[y, t]; 
initialdisplacement = 5; 
\(initialvel = 3.0; \)\)endtime = 30; Clear[s, y, notforcedy, t];
f[t_] = 2 Cos[2 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], VenetianRed}}, 
        PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Forced Oscillator", AspectRatio -> 1\/3]; 
    \)

This comes from:

forcedoscillator

with starter data y[0] = 1.5  and  y'[0]=3.0.
This forced oscillator is confused for a while and then settles into steady-state in the long run.  
How do you get your hands on plotting the long run steady state behavior?

It's simple!
 You just change the starter data to something really easy such as
       y[0] = 0 and y'[0]=0. 
And then you solve:

forcedoscillator

with starter data y[0] = 0 and y'[0]=0.

initialdisplacement = 0; initialvel = 0; endtime = 30; forcedoscillator=y''[t]+b y'[t]+c y[t]==f[t];zeroinputsolution=NDSolve[{forcedoscillator,y[0] == initialdisplacement,y'[0]==initialvel},y[t],{t, 0, endtime}];\(Clear[yzeroinput]; 

    yzeroinput[t_] = 
      y[t] /. zeroinputsolution[[1]]; 

    \(steadystateplot = 
      Plot[yzeroinput[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.02], GreenDark}}, 
        PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Steady State", AspectRatio -> 1\/3]; 
    \)\)

Compare with the plot of the original forced oscillator:

Show[steadystateplot, forcedplot, PlotLabel -> None];

In the long run they eventually settle into the same steady state oscillations. 
This will happen no matter what starter data you go with.
The zero input solution always captures the essence of the long-term action of any forced damped oscillator.

B.2.a.iii) Transient

Stay with the same oscillator and look at this plot of 
      yforced[t]  -  yzeroinput[t]:

transientplot = 
      Plot[forcedy[t] - yzeroinput[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.02], GeraniumLake}}, 
        PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Transient", AspectRatio -> (1/3)];

Why do folks call this the transient effect?
How do you plot this directly?

Answer:

It's called the transient effect because it dies out because of the effect of the damping term.
You can plot it directly by taking the original forced oscillator, keeping the same starter data, but zeroing out the forcing function f[t]:

endtime = 30; initialdisplacement = 1.5; initialvel = 3.0; Clear[s, y, yunforcedy, t];
unforcedoscillator=y''[t]+b y'[t]+c y[t]==0;
unforcedsolution=NDSolve[{unforcedoscillator,y[0] == initialdisplacement,y'[0]==initialvel},y[t],{t, 0, endtime}];\(yunforced[t_] = 
      y[t] /. unforcedsolution[[1]]; 


    directtransientplot = 
      Plot[yunforced[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, PlotRange -> All, 
        AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Transient", AspectRatio -> 1\/3]; \)

Compare:

Show[transientplot, directtransientplot];

On the money.

B.2.b.i)

Here's a new forced oscillator subjected to a constant forcing function:

b = Random[Real, {0.5, 1.4}]; c = Random[Real, {4.0, 6.0}]; initialdisplacement = -3.7; initialvelocity = 4.0; Clear[f, t, yforced];
f[t_] = 5.0;
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}];\(yforced[t_] = 
      y[t] /. forcedsolution[[1]]; 


    forcedplot = 
      Plot[yforced[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], NavyBlue}}, 
        PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Forced Oscillator", AspectRatio -> 1\/4]; 
    \)

Come up with plots of the zero input solution and the the unforced solution.
Show a plot of 
          yzeroinput[t]  +  yunforced[t]
together with the plot immediately above and describe what you see.

Here comes a plot of yunforced[t]:

Clear[y, t, transient];
unforcedoscillator=y''[t]+b y'[t]+c y[t]==0;
unforcedsolution=NDSolve[{unforcedoscillator,y[0] == initialdisplacement,y'[0]==initialvel},y[t],{t, 0, endtime}];\(yunforced[t_] = 
      y[t] /. unforcedsolution[[1]]; 


    transientplot = 
      Plot[yunforced[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], MediumBlue}}, 
        PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "unforced", AspectRatio -> 1\/3]; \)

Here comes a plot of yzeroinput[t]:

Clear[yzeroinput];
forcedoscillator=y''[t]+b y'[t]+c y[t]==f[t];
zeroinputsolution=NDSolve[{forcedoscillator,y[0] == 0,y'[0]==0},y[t],{t, 0, endtime}];yzeroinput[t_] = 
      y[t] /. zeroinputsolution[[1]]; 
    steadystateplot = 
      Plot[yzeroinput[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.02], GreenDark}}, 
        PlotRange -> All, AxesLabel -> {"t", "y[t]"}, 
        PlotLabel -> "Steady state", AspectRatio -> (1/3)];

Zap.

Here's a plot of yzeroinput[t]  +  yunforced[t]:

combinedaction = 
      Plot[yzeroinput[t] + yunforced[t], {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.02], HotPink}}, 
        PlotRange -> All, AxesLabel -> {"t", ""}, 
        AspectRatio -> (1/3)];

Compare to Mathematica's plot of the action of the forced oscillator:

Show[combinedaction, forcedplot, 
      AxesLabel -> {"t", ""}];

They are the same. It works this way everytime. 
If you want to see why, go on.

B.2.b.ii)  yforced[t] = yzeroinput[t] + yunforced[t]

In the last part it became apparent that
        yforced[t] = yzeroinput[t]  +  yunforced[t].
Explain why this happened and why you are guaranteed that this will happen anytime you are dealing with a forced linear (damped or undamped) oscillator.

Go with a forced linear oscillator with cleared forcing function and cleared coefficients:

Clear[y, f, b, c];
forcedoscillator=y''[t]+b y'[t]+c y[t]==f[t]

with cleared starter data y[0] = p and y'[0]=q. 
Facts to ponder:
→ yzeroinput[t] solves:

Clear[yzeroinput];
zeroinputeq=forcedoscillator/.{y[t] -> yzeroput[t],y'[t]->yzeroinput'[t],y''[t]->yzeroinput''[t]}

with starter data yzeroinput[0] = 0 and yzeroinput'[0]=0.

→ yunforced[t] solves:

Clear[yunforced];
unforcedeq=forcedoscillator/.{y[t] -> yunforced[t],y'[t]->yunforced'[t],y''[t]->yunforced''[t],f[t] -> 0}

with starter data yunforced[0] = p and yunforced'[0]=q. 
Adding the two equations and remembering that the derivatives of the sums are sums of derivatives, you see that when you put
        yforced[t] = yzeroinput[t]  +  yunforced[t],   
then you are guaranteed that yforced[t] solves the original forced oscillator differential equation
        yforced''[t] + b yforced'[t] + c yforced[t]=f[t]  +  0 = f[t].

Also this gives you 
        yforced[0] = 
        yzeroinput[0]  +  yunforced[0] = 0  +  p = p
and     
        yforced'[0]=yzeroinput'[0] + yunforced'[0]=0  +  q = q,
     
The upshot: yforced[t] solves

Clear[y, f, b, c];
forcedoscillator=y''[t]+b y'[t]+c y[t]==f[t]

with correct starter data y[0] = p and y'[0]=q. 

That's why you can count on 
        yforced[t] = yzeroinput[t]  +  yunforced[t]  
for any and all forced linear oscillators.

B.3)  Euler's identity:
           (E^((a + I b) t)) = (E^( at))(Cos[b t]  +  I Sin[b t]) where I = (Sqrt[-1])

B.3.a)

Look at these:

\(Clear[t]; 
ComplexExpand[E\^\(I\ t\)]\)

\(Clear[t]; 
ComplexExpand[E\^\(I\ 2\ t\)]\)

\(Clear[t]; 
ComplexExpand[E\^\(I\ π\ t\)]\)

b = Random[Real, {-3, 3}]; ComplexExpand[(E^(I b t))]

a = Random[Real, {-3, 3}]b = Random[Real, {-3, 3}]ComplexExpand[(E^((a + I b) t))]

What the heck is going on here?

Mathematica is spitting out Euler's identity.
Euler's  Identity says:
If a and b are real numbers and I  = (Sqrt[-1]), then 
             (E^((a + I b) t)) = 
        (E^(a t)) (Cos[b t] + I Sin[b t]).
Folks call this the complex exponential. It is programmed right into Mathematica, and you will want it to be programmed into your brain as well.
Check it out:

\(a = Random[Real, {\(-3\), 3}]; 
b = Random[Real, {\(-3\), 3}]; 

    ComplexExpand[E\^\(\((a + I\ b)\)\ t\)]\)(E^(a t)) (Cos[b t] + I Sin[b t])

It happens every time.

\(TextForm
    \`\(\(B  .3 . b)\)\ The\ first\ and\ second\ derivatives\ of\ 
      E\^\(\((a\  + \ I\ b)\)\ t\)\ \)\)

Check this out:

\(Clear[y, t, a, b, K]; 
y[t_] = K\ E\^\(\((a + I\ b)\)\ t\)\)

y'[t]

(y')'[t]

What's the message?
Explain why it happens.

The message is that when you go with
       y[t]  =  K (E^((a  +  I b) t))
then you are guaranteed that
       \(TextForm
      \`\(y\^′\)[t]\  = \ 
        \(\((a\  + \ I\ b)\)\ K\ E\^\(\((a\  + \ I\ b)\)\ t\)\ \  = \ 
          \((a\  + \ I\ b)\)\ y[t]\)\)
and
       y''[t]=  (((a  +  I b))^2) K (E^((a  +  I b) t))  =  
          (((a  +  I b))^2) y[t]
You take the derivative of the complex exponential the same way you were used to taking the derivative of any exponential back in calculus.

Start with
       \(TextForm
      \`y[t]\  = 
        \(K\ E\^\(\((a\  + \ I\ b)\)\ t\) = 
          K\ E\^\(a\ t\)\ \((Cos[b\ t]\  + \ I\ Sin[b\ t])\)\)\)

\(Clear[k, a, b, y, t]; 

    y[t_] = K\ E\^\(a\ t\)\ \((Cos[b\ t] + I\ Sin[b\ t])\)\)

Multiply it out:

Expand[y[t]]

Use the product rule to differentiate both parts with respect to t and multiply out, remembering that (I^2) = -1:

Expand[y'[t]]

Compare:

Expand[(a + I b) y[t]]

So when you start with   y[t] = K (E^((a  +  I b) t)) , you get  
      (y^′)[t]  =  
        K (a  +  I b) (E^((a  +  I b) t))

B.4) The characteristic equation of the unforced linear oscillator diffeq
        y''[t]+b y'[t]+c y[t]=0 
is 
         (z^2) + b z + c = 0.
If (Subscript[z, 1]) and (Subscript[z, 2]) solve the characteristic equation,  and (Subscript[K, 1]) and (Subscript[K, 2]) are constants, then 
         (Subscript[K, 1]) (E^((Subscript[z, 1]) t)) + (Subscript[K, 2]) (E^((Subscript[z, 2]) t)) 
is a solution of
            y''[t]+b y'[t]+c y[t]=0.
How to set (Subscript[K, 1]) and (Subscript[K, 2]) to handle starter data on y[0] and y'[0].

B.4.a.i) Using the solutions of the characteristic equation

Here's an unforced linear oscillator diffeq:

b = 0.45; c = 2.2; startery = -1.0; starteryprime = 2.0; Clear[t, y, Derivative];
oscdiffeq={y''[t]+b y'[t]+c y[t]==0,y[0] == startery,y'[0]==starteryprime};
ColumnForm[Thread[oscdiffeq]]

The characteristic equation for this oscillator diffeq is:

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

Use the solutions of the characteristic equation to come up with a formula for the solution of this oscillator diffeq.
Plot the solution and describe what you see.

Step1: Write down the characteristic equation:

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

Note the relation between the characteristic equation and the oscillator diffeq:

ColumnForm[Thread[oscdiffeq]]

Step 2: Solve the charactistic equation for z to set up the general solution of the diffeq

zsols = Solve[charequation, z]

Fish the solutions for z out:

z1  = zsols[[1, 1, 2]]

z2  = zsols[[2, 1, 2]]

Set up the general form of the solutions.

\(Clear[gensol, K1, K2]; 

    gensol[t_] = K1\ E^\((z1\ t)\) + K2\ E^\((z2\ t)\)\)

Step 3: Solve for the K1 and K2 that correspond to the given starting data on y[0] and y'[0].

\(ystarteq\  = \ \((gensol[0]\  == \ startery)\); 

    yprimestarteq\  = \ \((\(gensol'\)[0]\  == \ starteryprime)\); 

    Ksols\  = \ Solve[{ystarteq, yprimestarteq}, {K1, K2}]\)

Substitute these values of K1 and K2 in to get the raw form of the exact formula:

gensol[t] /. Ksols[[1]]

Make it look nice by hitting gensol[t] with the fundamental identity
        (E^((a  +  I b) t)) = 
        (E^(a t)) Cos[b t]  +  I (E^(a t)) Sin[b t]:

\(Clear[yformula]; 

    yformula[t_]\  = \ Chop[ComplexExpand[gensol[t] /. Ksols[\([1]\)]]]\)

See the solution :

\(ColumnForm[Thread[oscdiffeq]]
\)Plot[yformula[t], {t, 0, 30}, 
      PlotStyle -> {{Thickness[0.01], Red}}, 
      AspectRatio -> (1/GoldenRatio), PlotRange -> All, 
      AxesLabel -> {"t", "y[t]"}];

A damped oscillator.
If you want to see why this method works, then go on to part iv) below.
If you want another example, go on to the next part.

B.4.a.ii) A random unforced damped oscillator

Here's an a random unforced damped linear oscillator diffeq:

\(b\  = \ Random[Real, {0.1, 0.6}]; 
c\  = \ Random[Real, {1, 3}]; 

    startery\  = \ Random[Real, {\(-2\), 2}]; 

    starteryprime\  = \ Random[Real, {\(-2\), 2}]; 

Clear[t, y]; 

    linoscdiffeq\  = 
      \ {\(\(y'\)'\)[t]\  + \ b\ \(y'\)[t]\  + \ c\ y[t]\  == \ 0, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ y[0]\  == \ startery, 
        \(y'\)[0]\  == \ starteryprime}; 

    ColumnForm[Thread[linoscdiffeq]]\)

The characteristic equation for this oscillator diffeq is:

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

Use the solutions of the characteristic equation to come up with a formula for the solution of this oscillator diffeq.
Plot the solution and describe what you see.

Step1: Write down the characteristic equation:

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

Note the relation between the characteristic equation and the oscillator diffeq:

ColumnForm[Thread[linoscdiffeq]]

Step 2: Solve the charactistic equation for z to set up the general solution of the diffeq

zsols = Solve[charequation, z]

Fish the solutions for z out:

z1 = zsols[[1, 1, 2]]

z2 = zsols[[2, 1, 2]]

Set up the general form of the solutions.

\(Clear[gensol, K1, K2]; 

    gensol[t_] = K1\ E\^\(z1\ t\) + K2\ E\^\(z2\ t\)\)

Step 3: Solve for the K1 and K2 that correspond to the given starting data on y[0] and y'[0].

ystarteq = gensol[0] == startery; yprimestarteq=gensol'[0]==starteryprime;Ksols = Solve[{ystarteq, yprimestarteq}, {K1, K2}]

Substitute these values of K1 and K2 in to get the raw form of the exact formula:

gensol[t] /. Ksols[[1]]

Make it look nice by hitting gensol[t] with the fundamental identity
        (E^((a  +  I b) t)) = 
        (E^(a t)) Cos[b t]  +  I (E^(a t)) Sin[b t]:

\(Clear[yformula]; 

    yformula[t_] = 
      Chop[ComplexExpand[
          gensol[t] /. Ksols[[1]]]]\)

See the solution :

\(ColumnForm[Thread[linoscdiffeq]]
\)Plot[yformula[t], {t, 0, 30}, 
      PlotStyle -> {{Thickness[0.01], Red}}, 
      AspectRatio -> (1/GoldenRatio), PlotRange -> All, 
      AxesLabel -> {"t", "y[t]"}];

A damped oscillator.
If you want to see why this method works, then go on to part iv) below.
If you want another example, go on to the next part.

B.4.a.iii) A random unforced undamped oscillator

Here's an a random unforced undamped linear oscillator diffeq:

b  =  0;
c  =  Random[Real, {1, 3}];
startery  =  Random[Real, {-2, 2}];
starteryprime  =  Random[Real, {-2, 2}];

Clear[t, y];
\(linoscdiffeq\  = 
        \ {\(\(y'\)'\)[t]\  + \ b\ \(y'\)[t]\  + \ c\ y[t]\  == \ 0, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ y[0]\  == \ startery, 
          \(y'\)[0]\  == \ starteryprime}\);
ColumnForm[Thread[linoscdiffeq]]

The characteristic equation for this oscillator diffeq is:

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

Use the solutions of the characteristic equation to come up with a formula for the solution of this oscillator diffeq.
Plot the solution and describe what you see.

Step1: Write down the characteristic equation:

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

Note the relation between the characteristic equation and the oscillator diffeq:

ColumnForm[Thread[linoscdiffeq]]

Step 2: Solve the charactistic equation for z to set up the general solution of the diffeq

zsols = Solve[charequation, z]

Fish the solutions for z out:

z1 = zsols[[1, 1, 2]]

z2 = zsols[[2, 1, 2]]

Set up the general form of the solutions.

\(Clear[gensol, K1, K2]; 

    gensol[t_] = K1\ E\^\(z1\ t\) + K2\ E\^\(z2\ t\)\)

Step 3: Solve for the K1 and K2 that correspond to the given starting data on y[0] and y'[0].

ystarteq = gensol[0] == startery; yprimestarteq=gensol'[0]==starteryprime;Ksols = Solve[{ystarteq, yprimestarteq}, {K1, K2}]

Substitute these values of K1 and K2 in to get the raw form of the exact formula:

gensol[t] /. Ksols[[1]]

Make it look nice by hitting gensol[t] with the fundamental identity
        (E^((a  +  I b) t)) = 
        (E^(a t)) Cos[b t]  +  I (E^(a t)) Sin[b t]:

\(Clear[yformula]; 

    yformula[t_] = 
      Chop[ComplexExpand[
          gensol[t] /. Ksols[[1]]]]\)

See the solution :

\(ColumnForm[Thread[linoscdiffeq]]
\)Plot[yformula[t], {t, 0, 30}, 
      PlotStyle -> {{Thickness[0.01], Red}}, 
      AspectRatio -> (1/GoldenRatio), PlotRange -> All, 
      AxesLabel -> {"t", "y[t]"}];

A undamped unforced oscillator.
If you want to see why this method works, then go on.

B.4.a.iv)

Why does this method work?

Start with cleared coefficients b and c:

\(Clear[y, t, b, c]; 

    linoscdiffeq\  = \ 
      \((\(\(y'\)'\)[t]\  + \ b\ \(y'\)[t]\  + \ c\ y[t]\  == \ 0)\)\)

And do educated guesswork.
Go with numbers K and z and put 
        guessy[t] = K (E^(z t))
Because the diffeq is
          y''[t] + b (y^′)[t] + c y[t]=0, 
you want to select z so that 
           guessy''[t] + b guess'[t] + b guessy[t]=0
This is the same as
          
      K (z^2) (E^(z t))  +   b z K (E^(z t)) + 
          c K (E^(z t )) = 0
 This is the same as
          (z^2)  + b z  +  c  =  0 .
The upshot: 
If z1 and z2 solve   (z^2)  +  b z  +  c  =  0 , and K1 and K2 are any constants, then
then both
       K1 (E^(z1 t)) and K2 (E^(z2 t)) 
are guaranteed to solve
       y''[t] + b y'[t] + c y[t]=0.
And so is
       K1 (E^(z1 t)) +  K2 (E^(z2 t)) 
 guaranteed to solve
       y''[t] + b y'[t] + c y[t]=0.
 This leaves you two degrees of freedom to use to set K1 and K2 to solve
       K1 (E^(z1 t)) +  K2 (E^(z2 t)) =  
        given starting value on y[0]
and
       K1 z1 (E^(z1 t)) +  K2 z2  (E^(z2 t))=given starting value on y'[0].

B.4.a.v)

What happens in the disgusting (and rare) case that the characteristic equation (z^(2 )) +  b z  +  c = 0 has only one solution?

If the characteristic equation
       (z^2)  +  b z  +  c  =  0  (double root), you take the lone solution z1
  and form
       K1 (E^(z1 t))+ K2 t (E^(z1 t)).
This is guaranteed to solve
       y''[t] + b y'[t] + c y[t]=0.

 This leaves you two degrees of freedom to use to set K1 and K2 to solve to get the right starting values on y[0] and y'[0].
  Case in point:

\(b\  = \ 6; 
c\  = \ 9; 
startery\  = \ Random[Real, {1, 7}]; 

    starteryprime\  = \ Random[Real, {\(-2\), 2}]; 

Clear[t, y]; 

    linoscdiffeq\  = 
      \ {\(\(y'\)'\)[t]\  + \ b\ \(y'\)[t]\  + \ c\ y[t]\  == \ 0, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ y[0]\  == \ startery, 
        \(y'\)[0]\  == \ starteryprime}; 

    ColumnForm[Thread[linoscdiffeq]]\)

Step1: Write down the characteristic equation:

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

Note the relation between the characteristic equation and the oscillator diffeq:

ColumnForm[Thread[linoscdiffeq]]

Step 2: Solve the charactistic equation for z to set up the general solution

zsols = Solve[charequation, z]

Fish the solutions for z out:

z1 = zsols[[1, 1, 2]]

z2  = zsols[[2, 1, 2]]

Only one solution - this tells you to insert the extra t:

Set up the general form of the solutions.

Clear[gensol, K1, K2];
gensol[t_]=K1 (E^(z1 t))+K2 t (E^(z2 t))

Step 3: Solve for the K1 and K2 that correspond to the given starting data on y[0] and y'[0].

ystarteq = gensol[0] == startery; yprimestarteq=gensol'[0]==starteryprime;Ksols = Solve[{ystarteq, yprimestarteq}, {K1, K2}]

Substitute these values of K1 and K2 in to get the raw form of the exact formula:

gensol[t] /. Ksols[[1]]

Make it look nice by hitting gensol[t] with the fundamental identity
        (E^((a  +  I b) t)) = 
        (E^(a t)) Cos[b t]  +  I (E^(a t)) Sin[b t]:

\(Clear[yformula]; 

    yformula[t_] = 
      Chop[ComplexExpand[
          gensol[t] /. Ksols[[1]]]]\)

See the solution :

\(ColumnForm[Thread[linoscdiffeq]]
\)Plot[yformula[t], {t, 0, 20}, 
      PlotStyle -> {{Thickness[0.01], Red}}, 
      AspectRatio -> (1/GoldenRatio), PlotRange -> All, 
      AxesLabel -> {"t", "y[t]"}];

Damped out fast!

B.5) Trying to come up with formulas for solutions  of the forced oscillator diffeq  
y''[t]+b y'[t]+c y[t]=f[t] via convolution integrals

Here's how you use the convolution integral method to go after an exact formula for the oscillator coming from 
         y''[t] + 4.1 y'[t] + 7.2 y[t]=f[t] 
with  f[t] = 5 Sin[2  t] , y[0] = 3.0  and  y'[0]=2.2:

Step 1: Calculate yunforced[t]

This involves coming up with the formula for the unforced solution. 
This is the solution of the unforced damped oscillator
       y''[t] + 4.1 y'[t] + 7.2 y[t]=0 
with  y[0] = 3.0  and  y'[0]=2.2.
You already know how to do this:

\(Clear[f, t]; 
f[t_] = 5\ Sin[2\ t]; \((b = 4.1; )\)\)c = 7.2; ystarter = 3.0; yprimestarter = 2.2; \(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]]]]
        \)

You can see the effect of the damping term on this formula.
Put yunforced[t] aside for a moment.

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

This involves coming up with the formula of the solution of the unforced oscillator
       y''[t] + 4 y'[t] + 5 y[t]=0 
with  y[0]=0  and  y'[0]=1.
You already know how to do this:

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) yunitimpulse[t - x] f[x] 
            dx

This involves going after a formula for the zero input solution, yzeroinput[t], of the forced oscillator
       y''[t] + 4 y'[t] + 5 y[t]=f[t] 
with zeroed out starter data  y[0] = 0  and  y'[0]=0:
A miracle of calculus (to be explained later) is that you can get a formula for yzeroinput[t] by setting  
          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]]
        \)

Looks good.

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

Get your shot at an exact formula for the solution of 
       y''[t] + 4 y'[t] + 5 y[t]=f[t] 
with f[t] = 5 Sin[t],  y[0] = 3  and  y'[0]=2
by putting
       formulay[t] = yzeroinput[t]  +  yunforced[t]:

\(Clear[yformula]; 

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

Done!
See it:

Plot[yformula[t], {t, 0, 14}, 
      PlotStyle -> {{Thickness[0.01], Red}}, 
      AspectRatio -> (1/GoldenRatio), 
      AxesLabel -> {"t", "y[t]"}];

You can see that the forcing function  f[t] = 5  Sin[2  t] controls the steady state behavior.

B.5.a)

Use the convolution integral method to go after an exact formula for the oscillator coming from 
       y''[t] + 2.4 y'[t] + 6.7 y[t]=f[t]
with f[t] =  20  DiracDelta[t  -  4],  y[0] = 8.0  and  y'[0]=-2.7.

Here you go:
Copy paste and edit to get:

Step 1: Calculate yunforced[t]

This involves coming up with the formula for the unforced solution. 
This is the solution of the unforced damped oscillator
       y''[t] + 2.4 y'[t] + 6.7 y[t]=0 
with  y[0] = 8.0  and  y'[0]=-2.7.
You already know how to do this:

Clear[f, t];
f[t_]  = 30 DiracDelta[t  -  6];
b = 2.4;
c = 6.7;
ystarter = 8.0;
yprimestarter = -2.7;

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

You can see the effect of the damping term on this formula.
Put yunforced[t] aside for a moment.

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

This involves coming up with the formula of the solution of the unforced oscillator
       y''[t] + 2.4 y'[t] + 6.7 y[t]=0
with  y[0] = 0  and  y'[0]=1.
You already know how to do this:

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) yunitimpulse[t - x] f[x] 
            dx

This involves going after a formula for the zero input solution, yzeroinput[t], of the forced oscillator
       y''[t] + 2.4 y'[t] + 6.7 y[t]=f[t] 
with zeroed out starter data  y[0] = 0  and  y'[0]=0:
A miracle of calculus (to be explained later) is that you can get a formula for yzeroinput[t] by setting  
       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]]
        \)

Looks good.

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

Get your shot at an exact formula for the solution of 
       y''[t] + 2.4 y'[t] + 6.7 y[t]=f[t]
with f[t] =  20  DiracDelta[t  -  4],  y[0] = 8.0  and  y'[0]=-2.7.
by putting
       yformula[t] = yzeroinput[t]  +  yunforced[t]:

\(Clear[yformula]; 

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

Done!
See it:

Plot[yformula[t], {t, 0, 14}, 
      PlotStyle -> {{Thickness[0.01], Red}}, 
      AspectRatio -> (1/GoldenRatio), PlotRange -> All, 
      AxesLabel -> {"t", "y[t]"}];

The forcing function
       f[t] = 20  DiracDelat[t - 6] 
 goosed the oscillator at t = 6.

B.5.b)

Why do lots of old timers (those old profs who wear bowties) like to say that the convolution integral method is the method of dividing and conquering?

The goal is to come up with a formula for the oscillator coming from 
       y''[t]+b y'[t]+c y[t]=f[t] 
with given starting data 
       y[0] = p  and  y'[0]=q.
The convolution integral method divides the problem in two pieces:
yunforced[t] solves
       y''[t] + b y'[t] + c y[t]=0
with given starting data 
       y[0] = p  and  y'[0]=q.
yzeroinput[t] solves
       y''[t] + b y'[t] + c y[t]=f[t]
with starting data 
       y[0]=0  and  y'[0]=0.
The method conquers by putting 
       yformula[t]=yzeroinput[t] + yunforced[t]
which solves 
       y''[t] + b y'[t] + c y[t]=f[t] + 0=f[t]
with given starting data 
       y[0]=p + 0=p  and  y'[0]=q + 0=q.
Kinda neat and kinda sweet.

B.5.c)

What can break the back of the convolution integral method?

When you go with the convolution integral method to come up with a formula for the solution of
       y''[t] + b y'[t] + c y[t]=f[t] 
with given starting data 
       y[0] = p  and  y'[0]=q, 
the machine has to be able to do the convolution integral
       yzeroinput[t] = 
        (∫_0^x) yunitimpulse[t - x] f[x] 
            dt
where yunitimpulse[t] solves
       y''[t] + b y'[t] + c y[t]=0 
with  y[0] = 0 and y'[0]=1.
If f[t] is nasty enough, this integral can become painful or plain impossible.

B.5.d) The calculus behind the convolution integral method

What guarantees that the convolution integral method always works when the integrals can be done?

Calculus does!
There is only one issue to explain. 
The issue is to explain why 
       yzeroinput[t] = 
        (∫_0^t) yunitimpulse[t - x] f[x] 
            dx
 solves
       y''[t] + b y'[t] + c y[t]=f[t] 
with  y[0] = 0  and  y'[0]=0.
         
To this end, remember that yunitimpulse[t] solves
       y''[t] + b y'[t] + c y[t]=0
with  y[0] = 0  and y'[0]=1; 
 so that
→ yunitimpulse[0] = 0
→ yunitimpulse'[0]=1 and
→ yunitimpulse''[t] + b yunitimpulse'[t] + c yunitimpulse[t]=0.
Enter yzeroinput[t] into Mathematica:

\(Clear[f, t, yunitimpulse, yzeroinput, y, b, c, x]; 
	
    yzeroinput[t_] = 
      ∫\_0\%t\( yunitimpulse[t - x]\ f[x]\) \[DifferentialD]x\)

Let Mathematica calculate yzeroinput'[t] and yzeroinput''[t]:

yzeroinput'[t]

yzeroinput''[t]

Remembering that yunitimpulse[0] = 0 and yunitimpulse'[0]=1 , look at 

        yzeroinput''[t] + b yzeroinput'[t] + c yzeroinput[t]:

(yzeroinput''[t]  +  b yzeroinput'[t]  +  
          c yzeroinput[t])/.{yunitimpulse[0] -> 0,yunitimpulse'[0] -> 1}

This tells you that
          yzeroinput''[t] + b yzeroinput'[t] + c yzeroinput[t] 
 is the same as
       f[t]+(∫_0^t)f[x]( yunitimpulse''[t - x]+ b yunitimpulse'[t - x]+ c yunitimpulse[t - x])dx.
 
 This collapses to
       f[t] +(∫_0^t)f[x]( 0)dx = f[t].

The upshot: 
          yzeroinput''[t] + b yzeroinput'[t] + c yzeroinput[t] 
 is guaranteed to be f[t].
        
This explains why 
       yzeroinput[t] = 
        (∫_0^t) f[x] yunitimpulse[t - x] dx
is  guaranteed to solve
       y''[t] + b y'[t] + c y[t]=f[t].
 Now check whether yzeroinput[0] is guaranteed to be 0:

yzeroinput[t]

When you plug in t = 0, you get:

yzeroinput[0]

This checks.
Finally check whether yzeroinput'[0] is guaranteed to be 0:

yzeroinput'[t]

When you plug in t = 0, you get:

yzeroinput'[0]

And this is zero because  
       yunitimpulse[0] = 0.
Explanation complete and you're out of here.