Physics 319
Homework #8
Non-Linear

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.

Simple Pendulum

We have analyzed the motion of a simple pendulum several times now, but why let that stop us from looking at the system yet again, and again....So far we have used the small angle approximation leading to Sin[θ]≈θ, but we now want to keep the Sin[θ] in place.

On a separate sheet of paper show that the potential is given by U(θ)=mgl(1-Cos[θ]).
Now use the gradient ((Overscript[∇, ⇀])U=\(\(∂U\/∂r\) e\&^\_r + 
              \(1\/r\) \(∂U\/∂θ\) 
                e\&^\_θ\)) to find the force F(θ).
With this expression for the force, use the small angle approximation to find the approximated potential (This should be analogous to the case for a spring).

Now you're ready  to plot the two expressions for the potential to see how they compare.

potential[m_, g_, l_, θ_] :=  
      m g l (1 - Cos[θ])

p1 = 
      Plot[potential[1, 9.8, 0.1, θ], {θ, -π, π}]

p2 = 
      Plot[0.5*1 (9.8) 0.1*θ^2, {θ, -π, π}]

Show[p1, p2]

The Sin[θ] causes a non-linearity in the potential.  That is, a non-parabolic potential implies a non-linear restoring force.  This also means there is a non-linearity in the equation of motion, 
\(TraditionalForm
      \`\(θ(t)\)\&¨ + \(ω\_o\%2\) sin[θ(t)] = 
        0\).This equation is in general not analytically solvable without approximatons, although the period can be found with the use of elliptic integrals (See Marion and Thornton, Sec. 4.4).  We will look at the angular position versus time, using numerical techniques.

Let's continue to compare the approximate solution to the real solution.  So let's find θ(t).   We'll use the following initial conditions θ(0)=(45^o)and \(TraditionalForm\`\(θ(0)\)\&. = 0\) with an angular frequency of ω = (Sqrt[(g/l)]), using the above values for g and l.  From your notes you can write down the solution for the approximate case as just the SHO (θ(t)=A Cos[ω t - δ]).  Find the coefficients given the initial conditions and plot the function for several periods.

p3 =. ....

For the exact solution, let's use Mathematica's numerical solver, NDSolve[].

\(g = 9.8; 
l = 0.1; 

    s1 = NDSolve[{\(θ'\)[t] == ω[t], \ 
          \(ω'\)[t] + \((g/l)\) Sin[θ[t]] == 0, 
          θ[0] == π/2, ω[0] == 0}, {θ, ω}, {
          t, 0, 6}]\)

p4 = Plot[Evaluate[θ[t] /. s1], {t, 0, 2}]

Show[p3, p4]

How do the two graphs compare?  What can you say about the approximate solution?  Use the above procedure to make plots for different values of the inital angle, θ[0]=?, to estimate at what initial angle the approximation holds. 

	(answer here)
	
	Do your work making new graphs below by clicking in-between this cell and the next to start and new input cell....

So...for about what initial angle does the approximation look reasonable?  For how many periods?


	(answer here)

What happens to the exact waveform as time goes on?  What does this imply about the exact period?  Is this consistent with the trend of Eq. 4.29 in Marion and Thronton?

	(answer here)

We can track the separation of these two plots by taking the ratio of the values of the functions at some distant period.  Let's look at the 10th period and see how different the solutions are as a function of inital angle.  I've written this like a little program.  Initial parameter values are given up front, as well as initializing three number Lists (arrays for you programers).  The Table[] function will make what Mathematica calls a List.  I've used the Table commad here to just fill our three lists with zeros.  In the Do Loop we iterate through the calculations, each time through making the intial angle, (Subscript[θ, o]), a little bigger.  In the loop we access each new element in the List (array) and put the new calculation in the new position in the list.

\(\(g = 9.8; 
l = 0.1; 
τ = 2*Pi/Sqrt[g/l]; 

    exact = Table[0, {100}]; 
approx = Table[0, {100}]; 

    ratio = Table[0, {100}]; 
θ\_o = 0.0001; 

    Do[
				θ\_o = θ\_o + 0.01; 
				
      s4\  = NDSolve[{\(θ'\)[t] == ω[t], \ 
            \(ω'\)[t] + \((g/l)\) Sin[θ[t]] == 0, 
            θ[0] == θ\_o, ω[0] == 0}, {θ, 
            ω}, {t, 10*τ}]; 
				
      exact[\([i]\)] = Evaluate[θ[10*τ] /. s4]; 
				
      approx[\([i]\)] = θ\_o*Cos[Sqrt[9.8/0.1]\ *10*τ]; 
				
      ratio[\([i]\)] = exact[\([i]\)]/approx[\([i]\)]
		, {i, 1, 100}]

    \)\)

If you want to take a look a the lists of numbers, just type the name of the list without any brackets in a new cell and execute it.

Now we can plot the diff function.  Note that I use ListPlot[] since diff is a list of numbers.  The Table[] function often puts in extra {}'s so Flatten[] is used to take them out.  Anyway, here it is!

ListPlot[Flatten[ratio]]

Why does the curve start to turn around (It might help to look at a ListPlot[] of exact and approx -- Don't forget to Flatten[])?

	(answer here)

What is the approximate value of the initial angle when the exact and approximate solutions differ by about 10%? (You'll have to conver the horizontal axis from what is currently the number in the list of (Subscript[θ, o]) values to actual (Subscript[θ, o]) values.)

	(answer here)

To  finish up with the simple pendulum, let's make a phase plot.  We expect to see a nice oscillator.....

ParametricPlot[Evaluate[{θ[t], ω[t]} /.  s1], {t, 0, 6}]

A final question.  If you were armed only with the approximate solution, would you feel confident predicting the position of the pendulum at long times?  Why?  For what initial angles?

	(answer here)

Fun with numbers!

When you watched the movie "Chaos, Fractals and Dynamics: Conputer Experiments in Mathematics", with Robert Devaney, you saw how the process of iteration can lead to some strange results.  In this section of the notebook, we'll perform so of the calculations he showed us.

First, take your calculator and on a separate piece of paper, successively take the square root of a number.  Do you see the values converge to a single value?

	(answer here)

Now let's let Mathematica do the iterations for us.  Below we have a function that will iterate the (Sqrt[□]) function.  Again, we're using the Table function to create a list into which we will place the values found at each iteration.  The input to the function is the initial value and the number of iterations.

\(sroot[xo_, N_] := 
      \((
		slist = Table[0, {N}]; 
		x = xo; 
	
        Do[{x = Sqrt[x]; \ slist[\([i]\)] = x}, {i, 1, N}]
	)\)\)

Find the first 10 iterations for some starting value.

sroot[3, 10]

Now to look at the result, we just execute the list, slist.

slist

We can of course plot this using ListPlot[].  Note that I've increased the size of the points to make it more readable.

\(sroot[16, 20]; 

    ListPlot[Abs[slist], Prolog -> AbsolutePointSize[5], PlotRange -> {0, 5}]
      \)

Plat with this function.  Does this (Sqrt[□]) function always converge?  To what value?

	(answer here)

Now let's look at the Sine function (c*sin[x]).

\(sine[xo_, c_, N_] := 
      \((
		ylist = Table[0, {N}]; 
		x = xo; 
	
        Do[{x = c*Sin[x]; \ ylist[\([i]\)] = x}, {i, 1, N}]
	)\)\)

Choose an intial x-value and play with the plots for c ranging from 0.5-3, stepping by 0.2.  You DO NOT need to keep all of these graphs!!  Just keep changing the c-value in the same graph setup below and then describe what you see.

\(sine[1, 2.8, 30]; 
ListPlot[ylist, Prolog -> AbsolutePointSize[5]]\)

What seems to be happening?

	(answer here)

Play  with it some more right around 2.50-2.75.  What is happening to the number of solutions that are possible?  As you increase the number of iterations by 1, 2, 3 (i.e., as a function of the horizontal axis) what happens to the solution?

	(answer here)

Finally, let's look at the quadratic case (c*x*(1-x)).

\(quadratic[xo_, c_, N_] := 
      \((
		xlist = Table[0, {N}]; 
		x = xo; 
	
        Do[{x = c*x*\((1 - x)\); \ xlist[\([i]\)] = x}, {i, 1, N}]
	)\)\)

Again, play with the plot, choosing a starting value and then leaving that fixed.  Adjust the c-value from 0.5 to 4.0.  Again, just keep using the same setup below -- there's no need to have a thousand graphs here.

\(quadratic[0.3, 3.45, 100]; 

    ListPlot[xlist, Prolog -> AbsolutePointSize[4]]\)

How many solutions are there for c=1, c=2, c=3, c=4?

	(answer here)
	
Go back and play some more in the c=3-4 range.  Where does this function make the transistion from 2 to four solutions? eight solutions?

	(answer here)

For c<2, do you think you could predict the result after many iterations (think long time for a physical system)?  What about c=3?  What about c=3.5?  What about c>4? 

	(answer here)

What word would you use to describe the nature of the solutions for c>4?

	(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.