Lab 7

Name:
Student#:

1. The Simple Harmonic Oscillator

The equation for the harmonic oscillator (e.g. a mass attached to a spring) can be written as,
	(((d^2)x)/(d (t^2))) + (ω^2)x=0 	(1)
where ω is the angular frequency. If the oscillator has a mass m=1.5kg and is characterized by a spring constant k=22.0 N/m what is the frequency of oscillation for this system?

Remove["Global`*"]

k=

m=

ω=

f=

Answer:

Although (1) can be solved analytically using the command DSolve we will also explore a numerical solution using the command NDSolve. But first, solve this equation analytically with DSolve considering that the mass-spring system is initially stretched by 1 m and released from rest.

XE[t_]=x[t]/.DSolve[

Now, plot the solution from t=0 to t=3 s.  We will need this plot later, so we will store it in a variable plot1.

plot1=Plot[XE[t], {t, 0, 3}, AxesLabel -> {"t (s)", "x (m)"}]

Now solve this same problem numerically using the command NDSolve from  t=0 to t=3 s.

XN[t_]=x[t]/.NDSolve[

And plot the results

As well as the error

Plot[XN[t] - XE[t], {t, 0, 3}, AxesLabel -> {"t (s)", "Error in using NDSolve"}]

2. Interpolation

The InterpolatingFunction is a function that interpolates the values of x at any desired time from a discrete list of pairs {t, x}. For example, create the following list called data with the Table command, i.e. type
	data=Table[{t,(t^2)},{t,0,5}]

The values for x are only defined at some discrete time intervals, i.e. at t=0, t=1s,t=2 s etc. But what if we need the value of x at say t=1.376s? This can be done by interpolation, i.e. by fitting polynomial curves to the data. To interpolate, type the following command
	xt = Interpolation[data]

And now you can request the value at any time interval between t=0 and t=5s. Try 
	xt[1.376]

One advantage of interpolating the data is that you can use the command Plot instead of ListPlot to plot the data. Try:
	Plot[xt[t],{t,0,5},AxesLabel->{"t","x"}]

3. Interpolated Midpoint Scheme Solution

Now we will code our own solution to the harmonic oscillator equation using the midpoint scheme.  We begin by rewriting the second order equation as two first order equations.
((d v)/(d t))= - (ω^2)x=fv(x,v,t)
((d x)/(d t))= v = fx(x,v,t)

Let's define our right hand side functions

fv[x_,v_,t_]:=BLANK
fx[x_,v_,t_]:=v

And our parameters, and initial conditions

x[1] = BLANK;
v[1] = BLANK;
t[1] = BLANK;
T = 3;
h = 0.01;
nmax = Round[T/h];

Now the main loop of the Midpoint scheme (fill in the blanks)

Do[t[n+1]=n h+t[1];
sx=fx[x[n],v[n],t[n]];
sv= fv[  BLANKS ];
xtmp=(h/2)sx+BLANK;
vtmp=(h/2)sv+v[n];
sx=fx[BLANKS];
sv=fv[xtmp,vtmp,t[n]+h/2];
x[n+1]=BLANKS;
v[n+1]=h sv+v[n],
{n,1,nmax}]

And finally the interpolation

xM = Interpolation[Table[{t[n], x[n]}, {n, 1, BLANK}]];
vM = Interpolation[Table[{t[n], v[n]}, {n, 1, nmax + 1}]];

Plotting our results for x(t) and v(t)

Plot[xM[t], {t, 0, T}, AxesLabel -> {"t", "x"}]
Plot[vM[t], {t, 0, T}, AxesLabel -> {"t", "v"}]

As well as the error in x(t)

Plot[xM[t] - XE[t], {t, 0, 3}, AxesLabel -> {"t (s)", "Error in using Midpoint Scheme"}]

4. Damped Harmonic Oscillator

As seen in class, the equation for the linearly damped oscillator is:
   	(((d^2)x)/(d (t^2))) +γ((d x)/(d t)) + (ω^2)x=0
Solve this equation numerically for x(t) from t=0 to t=5 s using NDSolve given  m=1.5kg,  k=22.0 N/m,  γ=1.6 (s^(-1))and the initial conditions that the spring is initially stretched by 1 m and released from rest. Plot the solution.

plot2=Plot[

The frequency of the damped oscillator is not the same as the undamped oscillator althought they are generally close. The frequency of the damped oscillator depends upon γ and not just k and m as for the undamped oscillator. You can see that by plotting the motion of the damped and undamped oscillators over each other. To plot to graphs on a same figure use the command Show (type ?? Show and use the Help menu to learn about the command).

Show[plot1, plot2]

You should be able to see that the crests and throughs of the signals do not precisely align with each other, meaning that the frequencies aren't precisely the same.

5. Overdamped Harmonic

If you make γ very large then the oscillator becomes overdamped and you simply see the decay of the displacement, and no oscillation at all. Solve the same problem as before with γ=10 (s^(-1))and then see what happens.

6. Forced Oscillator

As has been or will be seen in class, the forced oscillator is described with

	(((d^2)x)/(d (t^2))) +γ((d x)/(d t)) + (ω^2)x=(Subscript[F, 0])cos ((Subscript[ω, f])t)
	
Solve this equation for x(t), and plot the solution, given m=1.5 kg, k=22 N/m, γ=0.5 (s^(-1)), (Subscript[F, 0])=1 N/kg and provided that the oscillator is initially at  equilibrium and at rest. Provide 3 solutions corresponding to forcing frequencies (a) (Subscript[ω, f])=ω/2 (b) (Subscript[ω, f])=ω and (c) (Subscript[ω, f])=2ω. Find the solution for the first 20 seconds using NDSolve.
(a)

(b)

(c)

Which of these 3 cases produces the largest amplitude?

Answer:

(When the forcing frequency is close to the natural frequency, large amplitudes are produced.  This is known as resonance.)