# Driven simple harmonic oscillator — amplitude of steady state motion

I am (partly as an exercise to understand Mathematica) trying to model the response of a damped simple harmonic oscillator to a sinusoidal driving force. I can solve the differential equation with some arbitrarily chosen boundary conditions, and get a nice graph;

params = {ν1 -> 1.0, ω1 -> 10.0, F -> 4.0};
system =
{D[x1[t], {t, 2}] == -ν1 D[x1[t], t] - ω1^2 x1[t] + F Cos[ω t], x1[0] == 1, x1'[0] == 0};
soln = DSolve[system /. params, x1[t], t][[1]][[1]];
Plot[x1[t] /. soln /. ω -> 8, {t, 1, 20}, Frame -> True, Axes -> False]


But I don't care about the transients - I just care about the steady state situation. I tried using Limit to extract this;

amp = Table[Max[Limit[x1[t] /. soln, t -> ∞]], {ω, 1, 20, 1}]
ListPlot[amp]


Looks a bit peculiar to me. Also, this is incredibly slow, and doesn't work symbolically.

I thought I could do something along the lines of forcing it to take as a solution of the DE $a Sin(\omega t + \phi)$;

params = {ν1 -> 40.0, ω1 -> 10.0, F -> 10.0};
x1 = a Sin[ω t + ϕ];
system = D[x1, {t, 2}] == -ν1 D[x1, t] - ω1^2 x1 + F Cos[ω t]
amp = Solve[system /. params, a]
phase = Solve[D[a /. amp, t] == 0, ϕ][[1]][[1]]


but this just turns into a mess and doesn't give the right result either.

Is there a canonical way to tackle this sort of problem? I don't really know what I'm doing with Mathematica yet, so any explanations would be gratefully received.

• What I would consider the canonical approach to a steady-state problem is to go to the Fourier domain and convert the differential equation to an algebraic one for the amplitude. Is this something you want, or are you determined to use the differential-equation approach?
– Jens
Sep 2 '14 at 4:21

params = {ν1 -> 1, ω1 -> 10, F -> 4};
system = {D[ x1[t], {t, 2}] == -ν1 D[x1[t], t] - ω1^2 x1[t] + F Cos[ω t], x1[0] == 1, x1'[0] == 0};
soln = DSolve[system /. params, x1[t], t][[1, 1]];

(* and the steady state is*)
lim = ((List @@ (Expand@soln[[2]])) /. x_ /; (Limit[x, t -> Infinity] == 0) :> 0)
steadyState = Simplify@Together[Plus @@ lim]
(* (-4 (-100 + ω^2) Cos[t ω] + 4 ω Sin[t ω])/(10000 - 199 ω^2 + ω^4) *)


$\frac{4 \omega \sin (t \omega )-4 \left(\omega ^2-100\right) \cos (t \omega )}{\omega ^4-199 \omega ^2+10000}$

And the amplitudes, with a nice resonance peak:

Plot[NMaxValue[steadyState, t], {ω, 1, 20}, PlotRange -> All]


• Note: In the general case a more careful condition is needed in the /; clause while calculating lim. The current condition is enough for this problem since there are no non-vanishing constant terms in the Limit[] Sep 1 '14 at 22:30
• That's great. The thing I learned here was the possibility and efficiency of breaking down equations into individual terms. Sep 2 '14 at 16:59
• @NLambert yes, but you should understand that it's an artificial way to circumvent the Limit[] limitations (pun intended) Sep 2 '14 at 17:03

It is much easier and more general.

system = {D[x1[t], {t, 2}] == -ν1 D[x1[t], t] - ω1^2 x1[t] + F Cos[ω t]};


without specifying neither the initial conditions nor the parameters (except of course v1 > 0, ω1^2 > 0) is solved by

xx[t_]=x1[t]/. DSolve[System, x1[t], t][[1]];


The behaviour of the solution at large t is simply given by letting die out all eponentially decaying terms.

Hence:

y[t_] = xx[t] /. Exp[t_] -> 0


$-\frac{4 F \left(-\omega ^2 \text{Cos}[t \omega ]+\text{$\omega $1}^2 \text{Cos}[t \omega ]+\text{$\nu $1} \omega \text{Sin}[t \omega ]\right)}{\left(\text{$\nu $1}^2+2 \omega ^2-2 \text{$\omega $1}^2+\text{$\nu $1} \sqrt{\text{$\nu $1}^2-4 \text{$\omega $1}^2}\right) \left(-\text{$\nu $1}^2-2 \omega ^2+2 \text{$\omega $1}^2+\text{$\nu $1} \sqrt{\text{$\nu $1}^2-4 \text{$\omega $1}^2}\right)}$

Hope this helps, Wolfgang