Forced Harmonic Oscillator

(ω^2) y(x)+2 b y'(x)+y''(x)=0

\(y(t) = 
          e\^\(\(-b\)\ t\)\ {
              c1\ e\^\(\(-Ω\)\ t\) + 
                c2\ e\^\(Ω\ t\)}

        Ω = \@\(b\^2 - ω\^2\)\)

Choose the constants c1 and c2 so that oscillator
starts out with the initial conditions y[0]=1, and
y'[0]=0,find c1=(-b+Ω)/2Ω; c2=(b+Ω)/2Ω;

Case I overdamped b>ω

ydamp[t_,b_,ω_]=(Exp[-b t](c1 Exp[-Omega t]+c2 Exp[Omega t])) //. {Omega->Sqrt[b^2-ω^2],c1->(-b+Omega)/2/Omega,c2->(b+Omega)/2/Omega}

Check initial conditions at t=0

ydamp[t,b,ω] /. t->0

velocity=Simplify[D[ydamp[t,b,ω],t]] /. t->0

Plot[ydamp[t,1.1,1],{t,0,10}]

Case II Oscillator ω > b

yosc[t_,b_,ω_]=ExpToTrig[(Exp[-b t](c1 Exp[-Omega t]+c2 Exp[Omega t])) //. {Omega->I Sqrt[-b^2+ω^2],c1->(-b+Omega)/2/Omega,c2->(b+Omega)/2/Omega}]

Plot[yosc[t,0.1,20],{t,0,10},PlotRange->All]

Lets add on a forcing term given by Sin[ω' t] on r.h.s of
differential equation

The particular solution is

yp[t_, b_, ω_, ωd_] = 
      ((-2 b ωd Cos[ωd t] + 
          (-ω(d^2) + (ω^2)) Sin[ωd t])/(4 
            (b^2) ω(d^2) + (((ω(d^2) - (ω^2)))^2)))

The general solution will be a sum of the
characteristic solution and this particular solution. In order to satisfy the initial conditions we must find new
c1 and c2 coefficients, however, for this example lets
just use the c1 and c2 coefficients that we used for the
unforced example, knowing that this implies a different set
of initial conditions. (We also consider the underdamped
case only).

y[t_,b_,ω_,ωd_]=yosc[t,b,ω]+yp[t,b,ω,ωd]

Plot[y[t,1,2,ωd] /. ωd->1,{t,0,20}]

Note the transient decaying away, followed by
a response to the driving oscillator

Plot[{y[t,1.0,10,1],Sin[t]},{t,0,20},PlotRange->{-0.1,0.1}]

Plot[{y[t,0.5,1.0,1],Sin[t]},{t,0,20}]

Plot[yp[10.0,0.1,ω,1],{ω,0.2,2.0},PlotRange->All]

Plot[yp[10.0,0.1,1,ω],{ω,0.2,2.0},PlotRange->All]