Period of Harmonic Oscillator using Numerical Methods

[Collisionman]

Homework Statement

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:

\frac{d^{2}}{dt^{2}} = - \omega^{2}_{0}x

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 \omega_{0} = 1.

Homework Equations

Euler-Richardson algorithm:

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

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

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:

\frac{dx}{dt} = 0
\frac{dv}{dt} = 1

However, the only thing is, I don't know how to determine the period numerically from this.

 

[jedishrfu]

Aren't you supposed to use the numerically generated data to determine the period?

 

[SteamKing]

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?

 

[Collisionman]

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?

Sorry, I should have written the initial equation of motion as:

\frac{d^{2}x}{dt^{2}}=-\omega^{2}_{0}x

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.

 

[Collisionman]

Aren't you supposed to use the numerically generated data to determine the period?

I know, but I don't know how to do this.

 

[D H]

Homework Equations

Euler-Richardson algorithm:

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

No! Let me fix that for you:

1. a_{n} = \frac{F(x_n,\,v_n,\,t)} m
2. V_{mid}=V_{n}+\frac{a_{n}Δt}{2}
   
3. X_{mid}=X_{n}+\frac{V_{n}Δt}{2}
   
4. a_{mid} = \frac{F(x_{mid},\,v_{mid},\,t+Δt/2)} m
   
5. X_{n+1} = X_{n}+\frac{V_{mid}Δt}{2}
   
6. V_{n+1} = V_{n}+\frac{a_{mid}Δt}{2}

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?

 

[Collisionman]

No! Let me fix that for you:

1. a_{n} = \frac{F(x_n,\,v_n,\,t)} m
2. V_{mid}=V_{n}+\frac{a_{n}Δt}{2}
   
3. X_{mid}=X_{n}+\frac{V_{n}Δt}{2}
   
4. a_{mid} = \frac{F(x_{mid},\,v_{mid},\,t+Δt/2)} m
   
5. X_{n+1} = X_{n}+\frac{V_{mid}Δt}{2}
   
6. V_{n+1} = V_{n}+\frac{a_{mid}Δt}{2}

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?

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: \frac{d^{2}x}{dt^{2}}=-\omega^{2}_{0}x, so I subbed in the angular frequency (\omega_{0} = 1) and x(at t=0)=1 to get my initial acceleration a_{n}. Using these equations which I previously wrote, the result I get resembles a cosine wave, wouldn't this correspond with the analytical solution x(t) = A(cos(t) + \phi)?

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?

 

[D H]

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=-ω₀x.

Your equation for a_mid 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.

 

[Collisionman]

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=-ω₀x.

So, am I right or wrong in saying that a_{n} is just -ω^{2}_{0}x, with x(t=0)=0 and ω₀ = 1?

Your equation for a_mid 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).

So how would I write a_{mid}?

 

[Chestermiller]

Let's go back to your original equation:

\frac{d^2x}{dt^2}=-\omega_0^2x

Let's rewrite it as follows:

\frac{d}{dt}\left(\frac{dx}{dt}\right)=-\omega_0^2x

But

\frac{dx}{dt}=v
So
\frac{dv}{dt}=-\omega_0^2x
But \frac{dv}{dt}=a
So, we have three equations, one for a, one for dv/dt, and one for dx/dt:

a=-\omega_0^2x
\frac{dv}{dt}=a
\frac{dx}{dt}=v

These are the equations you are supposed to use recursively in conjunction with the Euler Richardson finite difference formulas.

 

[Chestermiller]

In my judgement, the version of the algorithm appropriate to this particular problem should be:

1. a_n=-\omega_0^2x_n

2. v_{mid}=v_n+\frac{a_nΔt}{2}

3. x_{mid}=x_n+\frac{v_nΔt}{2}

4. a_{mid}=-\omega_0^2x_{mid}

5. v_{n+1}=v_n+a_{mid}Δt

6. x_{n+1}=x_n+v_{mid}Δt

Chet

 

[D H]

In my judgement, the version of the algorithm appropriate to this particular problem should be: ...

5. v_{n+1}=v_n+a_{mid}Δt

6. x_{n+1}=x_n+v_{mid}Δt

Chet

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.

 

[Collisionman]

In my judgement, the version of the algorithm appropriate to this particular problem should be:

1. a_n=-\omega_0^2x_n

2. v_{mid}=v_n+\frac{a_nΔt}{2}

3. x_{mid}=x_n+\frac{v_nΔt}{2}

4. a_{mid}=-\omega_0^2x_{mid}

5. v_{n+1}=v_n+a_{mid}Δt

6. x_{n+1}=x_n+v_{mid}Δt

Chet

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?

 

[SteamKing]

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.

 

[Chestermiller]

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.

 

[Collisionman]

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.

I would used EXCEL, but I am required to determine the frequency and therefore the period using only MATLAB.

 

[Collisionman]

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.

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.

 

[Chestermiller]

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.

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.

Chet

 

[SteamKing]

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.