Math 331.1 •  Fall 2005

Damped harmonic oscillator 
with sinusoidal forcing

Copyright © 2005 Murray Eisenberg.  All rights reserved.

The ODE

The ODE: 	y''  +  4 y'  +  5 y  =  sin(2 t)

Clear[y,  yp,  y1,  y2]

A particular solution and the general solution

General solution of the associated H ODE:

y[t]  /.  
      First @ DSolve[y''[t]  +  4 y'[t]  +  
              5  y[t]  ==  0,  y[t],  t]

yh[t_]  :=  
      k1 Exp[-2 t] Cos[t]  +  k2 Exp[-2 t] Sin[t]

Particular solution of the NH ODE:

yp[t_] := -(8/65) Cos[2  t]  +  (1/65) Sin[2  t]

Check:

yp''[t]  +  4  yp'[t]  +  5  yp[t]  ==  
        Sin[2  t]   //  Simplify

solutionp = Plot[{yp[t],  yp'[t]}, {t,  0,  π}, AspectRatio  ->  1, AxesLabel  ->  
              TraditionalForm  /@  {HoldForm[t],  HoldForm[y]}, PlotStyle  ->  {Red,  {Red,  
                  Dashing[{0.035,  0.025}]}}];

General solution of the NH ODE:

y[t_] := yh[t] + yp[t]
y[
      t]

Solutions for some sets of initial conditions

Solution for y(0)  =  0,  y' (0)  =  0:

Solve[{y[0]  ==  0,  y'[0]  ==  0},  {k1,  k2}]

y1[t_] :=  (8/65) Exp[-2 t] Cos[t]  +  
        (14/65) Exp[-2 t] Sin[t]  +  yp[t]

y1[t]

That's the solution we calculated by hand in class.

Here are plots of that solution, first of y vs. t (with the derivative y' also shown, dashed)…

solution1 = Plot[{y1[t],  y1'[t]}, {t,  0,  π}, AspectRatio  ->  1, AxesLabel  ->  
              TraditionalForm  /@  {HoldForm[t],  HoldForm[y]}, PlotStyle  ->  {Orange,  {Orange,  
                  Dashing[{0.035,  0.025}]}}];

…and now in the phase plane:

trajectory1 = ParametricPlot[Evaluate[{y1[t],  y1'[t]}], {t,  0,  π}, AspectRatio  ->  Automatic, AxesLabel  ->  
              TraditionalForm  /@  {HoldForm[y],  HoldForm[v]}, PlotStyle  ->  Orange];

Solution for y(0)  =  1/10,  y' (0)  =  2/10:

Solve[{y[0]  ==  (1/10),  y'[0]  ==  (2/10)},  {k1,  
        k2}]

y2[t_] :=  (29/130) Exp[-2 t] Cos[t]  +  
        (8/13) Exp[-2 t] Sin[t]  +  yp[t]

solution2 = Plot[{y2[t],  y2'[t]}, {t,  0,  π}, AspectRatio  ->  1, AxesLabel  ->  
              TraditionalForm  /@  {HoldForm[t],  HoldForm[y]}, PlotStyle  ->  {Magenta,  {Magenta,  
                  Dashing[{0.035,  0.025}]}}];

trajectory2 = ParametricPlot[Evaluate[{y2[t],  y2'[t]}], {t,  0,  π}, AspectRatio  ->  Automatic, AxesLabel  ->  {"y",  "v"}, PlotStyle  ->  Magenta];

Show[{solution1, solution2}];

Show[{trajectory1, trajectory2}];

Observe:  For this system, trajectories can intersect!

Question:  How is that possible?  Doesn't that contradict the Uniqueness Theorem?

Solution for y(0)  =  1/2,  y' (0)  =  0:

Solve[{y[0]  ==  (1/2),  y'[0]  ==  0},  {k1,  
        k2}]

y3[t_] :=  (81/130) Exp[-2 t] Cos[t]  +  
        (79/65) Exp[-2 t] Sin[t]  +  yp[t]

solution3 = Plot[{y3[t],  y3'[t]}, {t,  0,  2  π}, AspectRatio  ->  1, AxesLabel  ->  
              TraditionalForm  /@  {HoldForm[t],  HoldForm[y]}, PlotStyle  ->  {Blue,  {Blue,  
                  Dashing[{0.035,  0.025}]}}];

trajectory3 = ParametricPlot[Evaluate[{y3[t],  y3'[t]}], {t,  0,  π}, AspectRatio  ->  Automatic, AxesLabel  ->  {"y",  "v"}, PlotStyle  ->  Blue];

All three solutions:

Show[{solution1, solution2, solution3}];

Show[{trajectory1, trajectory2, trajectory3}];

Observe:  All these solutions—viewed either as functions y of t or as trajectories—eventually get attracted together.

Question:  Why does that happen?  And to what are they attracted?

Behavior for large t

moretrajectory1 = ParametricPlot[{y1[t],  y1'[t]}, {t,  0, 2 π}, AspectRatio  ->  Automatic, AxesLabel  ->  {"y",  "v"}, 
PlotStyle  ->  Orange, DisplayFunction  ->  Identity];
moretrajectory2 = ParametricPlot[{y2[t],  y2'[t]}, {t,  0,  2  π}, AspectRatio  ->  Automatic, AxesLabel  ->  {"y",  "v"}, PlotStyle  ->  Magenta, DisplayFunction  ->  Identity];
moretrajectory3 = ParametricPlot[{y3[t],  y3'[t]}, {t,  0,  2  π}, AspectRatio  ->  Automatic, AxesLabel  ->  {"y",  "v"}, PlotStyle  ->  Blue, DisplayFunction  ->  Identity];
Show[{moretrajectory1, moretrajectory2, moretrajectory3}, DisplayFunction  ->  $DisplayFunction];

Recall the phase-plane trajectory for the particular solution (Subscript[y, p]):

trajectoryp = ParametricPlot[{yp[t],  yp'[t]}, {t,  0,  π}, AspectRatio  ->  Automatic, AxesLabel  ->  {"y",  "v"}, PlotStyle  ->  Red];

Show[{moretrajectory1, moretrajectory2, moretrajectory3, trajectoryp}, DisplayFunction  ->  $DisplayFunction];

Plot[{y1[t], y2[t], y3[t], yp[t]}, {t,  0,  2  π}, AspectRatio  ->  1, AxesLabel  ->  
            TraditionalForm  /@  {HoldForm[t],  HoldForm[y]}, PlotStyle  ->  {Orange,  Magenta,  Blue,  Red}];

Observe:  As t increases, the solution curves (as functions y of t) tend toward the particular solution curve (Subscript[y, p]).  And in the phase plane, the trajectories get attracted to the trajectory for the particular solution (Subscript[y, p]).

Question:  Why?

Periodicity of the particular solution

Phase plane trajectory:

ParametricPlot[{yp[t],  yp'[t]}, {t,  0,  π}, AspectRatio  ->  Automatic, AxesLabel  ->  
            TraditionalForm  /@  {HoldForm[y],  HoldForm[v]}, PlotStyle  ->  Red];

Plot[yp[t], {t,  0,  3  π}, PlotStyle  ->  Red];

Observe:  The particular solution (Subscript[y, p]) is periodic—its solution curve in the phase plane is a closed oval.

Question:  Why is that so?  Why is this difference of trig functions apparently itself sinusoidal?

ϕ  =  θ  -  (π/2)  /.  θ  ->  
          ArcTan[-8]

ynew[t_]  :=  (1/(Sqrt[65])) Cos[2  t  +  ϕ]

yp[t]  ==  ynew[t]  //  Simplify

timeTicks  =  
        Transpose[{N[-ϕ Range[-1,  (3/2),  (1/2)]],  
            TraditionalForm  /@  {HoldForm[ϕ],  
                HoldForm[(ϕ/2)],  HoldForm[0],  
                HoldForm[-(ϕ/2)],  HoldForm[-ϕ],  
                HoldForm[-((3  ϕ)/2)]}}];
Plot[{(1/(Sqrt[65])) Cos[2  t], ynew[t]}, {t,  -(π/2),  ((3  π)/2)},  PlotStyle  ->  {Green,  Red}, Ticks  ->  {timeTicks,  Automatic}, TextStyle  ->  {FontSize  ->  16}];

Observe:  The particular solution (Subscript[y, p]) is "ahead" (or "behind") the sinusoidal function cos 2  t by (ϕ/2).  In other words, the "lag" (also called the "delay") between two adjacent peaks of (Subscript[y, p])(t) and cos 2  t is (ϕ/2).

Question:  Why is that lag (ϕ/2) rather than ϕ?