1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Period of Harmonic Oscillator using Numerical Methods

  1. Oct 17, 2013 #1
    1. The problem statement, all variables and given/known data

    Numerically determine the period of oscillations for a harmonic oscillator using the Euler-Richardson algorithm. The equation of motion of the harmonic oscillator is described by the following:

    [itex]\frac{d^{2}}{dt^{2}} = - \omega^{2}_{0}x[/itex]

    The initial conditions are x(t=0)=1 and v(t=0)=0 (i.e., position at time zero = 1 and velocity at time zero = 0) and [itex]\omega_{0} = 1[/itex].

    2. Relevant equations

    Euler-Richardson algorithm:

    1. [itex]a_{n} = \frac{dv}{dt}[/itex]
    2. [itex]V_{mid}=V_{n}+\frac{a_{n}Δt}{2}[/itex]
    3. [itex]X_{mid}=X_{n}+\frac{V_{n}Δt}{2}[/itex]
    4. [itex]a_{mid}=a_{n}+\frac{a_{n}Δt}{2}[/itex]
    5. [itex]X_{n+1} = X_{n}+{V_{mid}Δt}{2}[/itex]
    6. [itex]V_{n+1} = V_{n}+{a_{mid}Δt}{2}[/itex]

    Where X, V and a are the position, velocity and acceleration, respectively.

    3. The attempt at a solution

    I have already written a Matlab programme to numerically solve the differential equation describing the motion using the above algorithm. I used:

    [itex]\frac{dx}{dt} = 0[/itex]
    [itex]\frac{dv}{dt} = 1[/itex]

    However, the only thing is, I don't know how to determine the period numerically from this.
  2. jcsd
  3. Oct 17, 2013 #2


    Staff: Mentor

    Aren't you supposed to use the numerically generated data to determine the period?
  4. Oct 17, 2013 #3


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    I'm confused by your choice of initial conditions. The problem statement said to use x = 1 and v = 0, both at t = 0.
    You appear to be using v = 0 and a = 1. Why the change?

    When you generated your values of position versus time, did any periodic behavior manifest itself?
  5. Oct 18, 2013 #4
    Sorry, I should have written the initial equation of motion as:


    Basically, I forgot the 'x' on the LHS when I was writing the OP. This will mean that a = - 1 (I forgot to add the 'minus' sign), because x (t=0) = 1 and w = 0. The velocity remains v (t=0) = 0.

    Also, yes, periodic motion did manifest, and corresponded well with cos(t) for a time-step of 0.001.
    Last edited by a moderator: Oct 18, 2013
  6. Oct 18, 2013 #5
    I know, but I don't know how to do this.
  7. Oct 18, 2013 #6

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    No! Let me fix that for you:
    1. [itex]a_{n} = \frac{F(x_n,\,v_n,\,t)} m[/itex]
    2. [itex]V_{mid}=V_{n}+\frac{a_{n}Δt}{2}[/itex]
    3. [itex]X_{mid}=X_{n}+\frac{V_{n}Δt}{2}[/itex]
    4. [itex]a_{mid} = \frac{F(x_{mid},\,v_{mid},\,t+Δt/2)} m[/itex]
    5. [itex]X_{n+1} = X_{n}+\frac{V_{mid}Δt}{2}[/itex]
    6. [itex]V_{n+1} = V_{n}+\frac{a_{mid}Δt}{2}[/itex]

    Note the minor corrections to items 1, 5, and 6 and the *huge* correction to item 4.

    With regard to items 1 and 4, what is the force F in your harmonic oscillator?
  8. Oct 18, 2013 #7
    The way I initially wrote 5 and 6 here were typos, sorry. I used item 1 because that is basically what the initial acceleration is, as stated in the question: [itex]\frac{d^{2}x}{dt^{2}}=-\omega^{2}_{0}x[/itex], so I subbed in the angular frequency ([itex]\omega_{0} = 1[/itex]) and [itex]x[/itex](at [itex]t=0[/itex])[itex]=1[/itex] to get my initial acceleration [itex]a_{n}[/itex]. Using these equations which I previously wrote, the result I get resembles a cosine wave, wouldn't this correspond with the analytical solution [itex]x(t) = A(cos(t) + \phi)[/itex]?

    As for the force F in the revised algorithm you stated, I think it's something in terms of gravity and mass.

    As for getting the period of oscillations, as I initially wanted to find in my OP, would I be correct in saying that I'd need to use some sort of a fourier transform function in Matlab to get the fundamental frequency(ies) and then using that to compute the period with Period = 1/Frequency?
    Last edited: Oct 18, 2013
  9. Oct 18, 2013 #8

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Ignore gravity. The equations I posted are the generic form for Euler-Richardson integration. Force divided by mass *is* acceleration, and acceleration is the time derivative of velocity. You don't need force if you know that dv/dt=-ω0x.

    Your equation for amid in the opening post is incorrect. One look at the units should tell you that it cannot be correct. The first term on the right hand side has units of acceleration but the second has units of velocity. You can't add acceleration and velocity any more than you can add length and time. Units have to be commensurable in when adding, comparing, or assigning (e.g., F=ma is an assignment).

    Perhaps even more importantly, the algorithm posted in the original post misses the key reason for using Euler-Richardson or any other multistage numerical integration technique: Obtaining a better estimate of *the* intermediate derivative that perfectly advances state.
  10. Oct 18, 2013 #9
    So, am I right or wrong in saying that [itex]a_{n}[/itex] is just [itex]-ω^{2}_{0}x[/itex], with x(t=0)=0 and ω0 = 1?

    So how would I write [itex]a_{mid}[/itex]?
  11. Oct 18, 2013 #10
    Let's go back to your original equation:


    Let's rewrite it as follows:



    But [tex]\frac{dv}{dt}=a[/tex]
    So, we have three equations, one for a, one for dv/dt, and one for dx/dt:


    These are the equations you are supposed to use recursively in conjunction with the Euler Richardson finite difference formulas.
    Last edited: Oct 18, 2013
  12. Oct 18, 2013 #11
    In my judgement, the version of the algorithm appropriate to this particular problem should be:

    1. [itex]a_n=-\omega_0^2x_n[/itex]

    2. [itex]v_{mid}=v_n+\frac{a_nΔt}{2}[/itex]

    3. [itex]x_{mid}=x_n+\frac{v_nΔt}{2}[/itex]

    4. [itex]a_{mid}=-\omega_0^2x_{mid}[/itex]

    5. [itex]v_{n+1}=v_n+a_{mid}Δt[/itex]

    6. [itex]x_{n+1}=x_n+v_{mid}Δt[/itex]

  13. Oct 18, 2013 #12

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Ouch. What I posted is wrong. I only focused on steps 1-4. That's all for nought if the final steps are incorrect. Yours are the correct updates.
  14. Oct 18, 2013 #13
    I just realized this a second ago. Thank you, Chestermiller and D H, this is great help.

    As per my original question, does anything know how I find the period of the oscillations? Do I just use a fourier transform of some kind in MATLAB to find the frequency and then use that to find the period using T = 1/f?
    Last edited: Oct 18, 2013
  15. Oct 18, 2013 #14


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    You don't have to do a FT to find the period of oscillation. Your two eyeballs should do just fine.
    Remember, if the displacement is periodic, the values will repeat through each period. What you have to do is use your algorithm to generate a sufficient number of displacement values until you notice that they start to repeat. Plotting these values versus time should make the periodic nature of the displacement clear. Using either format, you should be able to determine the period of oscillation.
  16. Oct 19, 2013 #15
    This calculation is very easy to set up on an EXCEL spread sheet. You can then easily plot the results using the EXCEL graphics capabilities.
  17. Oct 20, 2013 #16
    I would used EXCEL, but I am required to determine the frequency and therefore the period using only MATLAB.
  18. Oct 20, 2013 #17
    I think I need to determine it more accurately than just using my eyeballs. Later on in the problem sheet I've given, I am required to include damping forces, which will make it harder to determine period using the method you've described. I am required to modify the MATLAB program I've written to determine the period, so I guess I'll have to use some sort of FT function in MATLAB, unfortunately.
  19. Oct 20, 2013 #18
    I don't have much experience with MATLAB, but, if MATLAB has a curve fitting capability that can fit a user-supplied functional form, you can use MATLAB to fit the generated results of the calculations to a damped harmonic functional form so it can determine the period, the damping, and the phase shift automatically. At least, that is what I would try to do.

  20. Oct 21, 2013 #19


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Even with all kinds of damping applied, harmonic motion is still harmonic motion. All the damping is going to do is change the amplitude of the motion. You still can use a table of motion data which you generate from your numerical method to figure out the period. You just need to extract the time between motion peaks.

    I would be surprised if MATLAB couldn't at least plot your motion data.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted