Physics 319
Homework #7
Oscillators

Double click on the ⌉ symbol with an arrow on the bottom at the far right to open the subsection of interest.  Double clicking again will close it back up.

Damped, Driven Harmonic Oscillator

In class, we analytically solved for the Damped, Driven, Harmonic Oscillator finding the complimentary solution and the particular solution, the complete solution being the sum of the two.  Let's numerically solve the Damped Driven Harmonic Oscillator, following the equations from our notes, where the driving term is given by (Subscript[A, d]) cos((Subscript[ω, d]) t).  We'll write things in such a way that we can make parametric plots of x, v, and ϕ.  That way we can make Phase Plots.  The parameter values are given below.  Notice that the two of the normal unknown constants
from the two solutions (amplitude of complimentary solution and phase of particular solution)  are contained in the inital conditions of x and v.

\(β = 1; 
ωo = 10; 
Ad = 1; 
ωd = 2; 

    s1 = NDSolve[{\(x'\)[t] == v[t], 
          \(v'\)[t] == 
            \(-2\)\ β*v[t]\  - \ ωo^2*x[t]\  + 
              Ad*Cos[ϕ[t]], \(ϕ'\)[t] == \ ωd, \ 
          x[0] == ϕ[0] == 0, v[0] == 0.3}, {x, v, ϕ}, {t, 0, 10}]
        \)

On a separate sheet of paper, show that the above expressions are consistent with the equation derived in class....

Now, let's plot the position versus time.

Plot[Evaluate[x[t] /.  s1], {t, 0, 10}]

Why does the graph look like this?  What happens over time to the strange behavior at the beginning?  Give detailed descriptions of the various parts of the curve and the reasons for them.  (Hint:  think about the complimentary and particular solutions and the various parameter values.)   Calculate how long until the intial oscillatory amplitude is down to 5% of the inital amplitude.  Calculate the amplitude of the final motion from the solution derived in class.  You can answer these questions on a separate sheet of paper.

Now, let's look at the phase plot.

ParametricPlot[Evaluate[{x[t], v[t]} /.  s1], {t, 0, 10},  
      PlotRange -> {{-0.04, 0.04}, {-0.31, 0.31}}]

This curve crosses itself many times.  Why doesn't it get into a stable orbit the first time it crosses itself?  What is different about this system than the others we've looked at?

	(answer here)

It is a little confusing to follow the above curve to tell for instance when it crossess itself.  We can spread this out by plotting along another axis.  Our other variable is 
ϕ(t).  

What does ϕ(t) represent?

	(answer here)
	
Let's make a 3D parametric plot.....

p1 = 
      ParametricPlot3D[Evaluate[{ϕ[t], x[t], v[t]} /.  s1], {t, 0, 10}, 
         BoxRatios -> {4, 1, 1},  PlotPoints -> 400,  
        AxesLabel -> {ϕ, x, v}]

Solve another example choosing parameters such that the complimentary solution does not oscillate.  Plot the position versus time.  (Note -- the cut and paste feature may be helpful here.)

\(\(β =. .\); 
\(ωo =. .\); 
\(Ad =. .\); 

    \(ωd =. .\); 
s2 = NDSolve[.....]\)

Plot[Evaluate[x[t] /.  s2], {t, 0, 10}]

Describe the curve and the reasons for it's shape....

	(answer here)

Fourier Series and The DDHO

We discussed in class the possibility of letting the driving term be any arbitrary function.  Fourier tells us that any periodic function can be described by a linear sum of harmonic functions.  Since our DDHO equation of motion is linear, the linear sum of harmonic functions leads to a superposition of solutions of the form of a complimentary solution and a particular solution that we have essentially aready found.

Let's drive our DDHO with a triangular wave: 

	F(t) = {({{A((1/2) + (t/τ))     for    - (τ/2) < t
                   < 0}, {A((1/2) - (t/τ))      for        0 < t
                   < (τ/2)}})
We can plot it below, to make sure it looks like a triangular wave.  Note that we've only defined it over one period, but that's all you need to find the Fourier Series.

\(\(triangle[A_, τ_, t_] := 
      \((
If[t ≤ 0, f = A*\((0.5 + t/τ)\), 
          f = A*\((0.5 - t/τ)\)]; f)\)	\)\)

Plot[triangle[4, 1, t], {t, -1, 1}]

On a separater sheet of paper find an expression for the coefficients ((Subscript[a, n]) and (Subscript[b, n])) of the Fourier Series that represents this function.  Be careful with the limits of integration for the two parts of F(t) in each of the coefficient calculations  If you need a particular integral, you can find it here.  (e.g., Integrate[t Cos[n w t],t] )....

You should have gotten a Cosine series with coefficients equalt to (A/π)(n^2)for only odd values of n.  Let's take a look at that Fourier Series.

Let A=π, so we can just add up (1/(n^2)) cos(nωt)for odd n.  Let's just plot the first two terms first with ω=π.

\(ω = π; 
Plot[....., {t, \(-1\), 1}]\)

Plot the next two terms.

\(ω = π; 
Plot[.......\ \(, {t, \(-1\), 1}\)]\)

We can get even fancier by using the Sum[] function.  Take a look at the syntax and use Sum[] plot the first 100 terms.

?? Sum

Note that you only want odd terms, so you will have to get the indexing of the Sum[] function correct.  Also, go ahead and plot over a wider range (say t=-4,4s); this Fourier Series will repeat itself.

\(ω = π; 
Plot[.....]\)

Before we put this into our Damped Driven Harmonic  Oscillator, let's put in a simpler expression.  Using the DDHO routine at the beginning of this notebook, put in a driving function of Ad (Cos[ϕ[t]]+3 Cos[3 ϕ[t]]).

\(β = 1; 
ωo = 10; 
Ad = 1; 
ωd = 2; 

    s3 = NDSolve[....., {t, 0, 10}]\)

Plot x[t] vs. t.

Plot[.....], {t, 0, 10}]

How does the amplitude here compare with the case at the beginning of the notebook?

	(answer here)

Can you pick a driving frequency that will make the amplitude even bigger?  Try it.

\(β = 1; 
ωo = 10; 
Ad = 1; 
\(\(ωd =. \)....\); 

    s4 = NDSolve[.....]\)

Plot[Evaluate[x[t] /.  s4], {t, 0, 10}]

Is this a bigger amplitude?  Why or why not?

	(answer here)
	
This looks almost sinusoidal.  Why - despite the complexity of the driving function?	

	(answer here)

Now, let's put in the triangular wave.

.....

Plot x[t].

And just for fun, let's look at the 3D parametric plot (I'm assuming you've called your last calculation s5.

p2 = 
      ParametricPlot3D[Evaluate[{ϕ[t], x[t], v[t]} /.  s5], {t, 0, 10}, 
         BoxRatios -> {4, 1, 1},  PlotPoints -> 400,  
        AxesLabel -> {ϕ, x, v}]

Why is the final motion not sinusoidal?

	(answer here)

When you are ready to print the whole thing out, please do the following to minimize paper usage
•Go back through the document and clean things up, like deleting unnecessary returns, trial runs, and pesky messages about variables being named like other variables, etc...
•Then click on Format on the Toolbar and choose Printing Style Environment -> Printout.
•Under File click Printing Settings -> Page Setup and make sure you're set for Portrait.
•Now under File click Print.  In the print dialog box choose the Layout tab and select 2 pages per sheet.
•Click Print.