8
$\begingroup$

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]

SHM transients

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]

response

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.

| improve this question | | | | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Jens Sep 2 '14 at 4:21
8
$\begingroup$
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]

Mathematica graphics

| improve this answer | | | | |
$\endgroup$
  • $\begingroup$ 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[] $\endgroup$ – Dr. belisarius Sep 1 '14 at 22:30
  • $\begingroup$ That's great. The thing I learned here was the possibility and efficiency of breaking down equations into individual terms. $\endgroup$ – NLambert Sep 2 '14 at 16:59
  • 1
    $\begingroup$ @NLambert yes, but you should understand that it's an artificial way to circumvent the Limit[] limitations (pun intended) $\endgroup$ – Dr. belisarius Sep 2 '14 at 17:03
2
$\begingroup$

It is much easier and more general.

Your equation

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

| improve this answer | | | | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.