I have the following code written to simulate a driven damped oscillator

Clear[x, v, K];
K = 5.0; 
m = 0.75; 
b = 0.05; 
Fo = 0.75; 
Ω = 2.80;
x[0] = 0.; v[0] = 0.;

FinalTime = 150; 
Δt = 0.005; 
Steps = Round[FinalTime/ Δt];

Do[t = n*Δt;
  F[n] = - K*x[n] - b*v[n] + Fo*Sin[Ω*t];
  v[n + 1] = v[n] + (F[n]/m)*Δt;
  x[n + 1] = x[n] + v[n + 1]*Δt,
  {n, 0, Steps}];

xdata = Table[{n*Δt, x[n]}, {n, 0, Steps}];

ListPlot[xdata, AxesLabel -> {"t", "x[t]"},
 PlotLabel -> "Driving Frequency = 2.80 ",
 PlotStyle -> PointSize[0.0005]]

I am a newbie to Mathematica, so understand that first; the class this is for is more a mechanics lab than Mathematica course (supposedly). Anyhow. I get with this a nice graph of the oscillator, and I can see where it eventually settles into steady-state.

I wanted to know how to do a couple of things.

  1. Move the start of the graph from zero to the right -- say start it at 100 -- so that the steady-state amplitude and frequency are easier to see.

  2. Show where the maxima are (part of this exercise is the relation between riving frequency and amplitude)

  3. I'd really like to display more than one driving frequency (the capital omega) on the same graph, showing that the eventual steady-state frequency depends on the mass only.

I hope it's clear what I want to do. The reason is partly that while I can take a measure of the steady-state frequency by just measuring it from the screen I wanted a better way to do that.

Anyhow, any help is appreciated. If I have approached this in entirely the wrong way let me know. And I hope this hasn't been too redundant -- i couldn't find anything on here that seemed to match exactly what i needed to do and this is one of those situations where I am learning this on the fly.

  • 1
    $\begingroup$ Don't use K as a variable, it is a reserved system symbol and if you assign a value to it, it'll cause problems sooner or later. Generally, if you avoid using symbols starting with a capital, then there won't be an accidental conflict. $\endgroup$
    – Szabolcs
    Mar 20 '14 at 15:32
  • $\begingroup$ Do you need to write your own ODE solver? NDSolve is less trouble and has much better algorithms. $\endgroup$
    – Szabolcs
    Mar 20 '14 at 15:34
  • $\begingroup$ don't need to write an ODE solver per se. And I am a beginner, remember -- i had never touched mathematica before three weeks ago. Point is, I was just hoping to have a display-type thing more than a solver, if you see what I mean -- just a little point that says "maximum here!" or something like. Does that make sense? and thx for the note about capital letters... $\endgroup$
    – Jesse
    Mar 20 '14 at 15:46
  • $\begingroup$ Is it the left end of the graph or the origin that you wish to move? If the former, do so by changing the interval over which you plot; if the latter, use option AxesOrigin inside Plot (the 4th example on the Documentation Center page ref/AxesOrigin will help you see what to do. $\endgroup$
    – murray
    Mar 20 '14 at 17:20
  • $\begingroup$ i want the time (t) to essentially start "later." so instead of starting at t=0 it would be as though I plugged in t=whatever. $\endgroup$
    – Jesse
    Mar 20 '14 at 17:23
   x[t] /. NDSolve[{m x''[t] + b x'[t] + k x[t] == 
       Fo Sin[\[CapitalOmega] t], x[0] == x0, x'[0] == v0}, 
     x, {t, tmin, tmin + deltaT}]], {t, tmin, tmin + deltaT}], {{m, 
   0.75}, 0.5, 2}, {{b, 0.05}, 0, 1}, {{k, 5}, 1, 10}, {{Fo, 0.75}, 0,
   2}, {{x0, 0}, -2, 2}, {{v0, 0}, -1, 1}, {{tmin, 0}, 0, 
  100}, {{deltaT, 10}, 10, 1000}]

enter image description here

  • $\begingroup$ thanks, so say I wanted to adjust $\CapitalOmega$. And I wanted the left side of the graph to start at t=whatever, sort of like zooming in but just so that I can move to the part where it goes steady-state. OR: I want to put more than one omega on there at once, to see the effect of changing it on the steady state frequency... $\endgroup$
    – Jesse
    Mar 20 '14 at 17:27

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.