Physics 319
Homework #9
Non-Linear Dynamics II

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.

Modified van der Pol Equation

We discussed in class a version of the van der Pol equation.  Here we will look at it in a bit more detail.
The equation is:
	\(TraditionalForm\`x\&¨ - \(μ(
              x\_o\%2 - x\^2 - x\& . \^2\/β\^2)\) 
            x\& .  + \(ω\_o\%2\) x = 0\)
Rewriting as two first-order differential equations, we have
	\(TraditionalForm\`x\& .  = v\)
\(TraditionalForm\`v\& .  - \(μ(
              x\_o\%2 - x\^2 - v\^2\/β\^2)\) v + \(ω\_o\%2\) x = 
        0\)
On a separate sheet of paper show that if we let μ=1, (Subscript[x, o]) = 1, (Subscript[ω, o]) = 1, β=1, the following is a solution for x: x(t)=cos(t).  Do this by just plugging in for x and v and verifying that the equation goes to 0.


Let's numerically solve this so we can make some plots.  Let μ=1, (Subscript[x, o]) = 1, (Subscript[ω, o]) = 1, β=1.  We want to plot this for different starting values of x (xinit), so I've included a parameter xinit in the input.

vanderpol[μ_, xo_, ωo_, β_, 
        xinit_] := (
s1 = 
          NDSolve[{x'[t] == 
                v[t],  v'[
                  t] == μ (xo^2 - x[t]^2 - v[t]^2/β^2) 
                    v[t] - ωo^2*x[t], x[0] == xinit, 
              v[0] == 0}, {x, v}, {t, 0, 80}, 
            MaxSteps -> 4000]
)

Plot x[t]  for a few different values of the starting position, xinit, -- some with xinit<(Subscript[x, o]) -- some with xinit>(Subscript[x, o]).

vanderpol[1, 1, 1, 1, 0.1];
Plot[Evaluate[x[t] /. s1], {t, 0, 80}]

What is happening to the amplitude as time goes on if xinit<(Subscript[x, o])?

	(answer here)

What happens if xinit>(Subscript[x, o])?

	(answer here)
	
What do you expect the phase plot to look like for these to cases?	
	
	(answer here)

Now make a Phase Plot for a couple cases.  Play around with different values of μ and xinit.  Try and get the convergence to happen slowly, so you see lots of spirals and a clear convergence to a circle.  Because the motion is attracted to this region of phase space, and eventually occupies only this region of phase space, this circle is the attractor for this equation.  Show plots for xinit<(Subscript[x, o]) and xinit>(Subscript[x, o])

vanderpol[0.1, 1, 1, 1, 0.5];
ParametricPlot[Evaluate[{x[t], v[t]} /. s1], {t, 0, 80},  
      AspectRatio -> 1]

Bifurcation

In this section we will look at bifurcation, or the period doubling route to chaos.  We will use the "logistic" equation which is the "quadratic" function from HW8.  Let's look at that again briefly before looking at bifurcation.

Here is the function from HW8 [x=c(1-x)x].  Since we use this equation iteratively, it might be better to write it as
	(Subscript[x, n + 1]) = c(1 - (Subscript[x, n])) (Subscript[x, n]).
That is, each next value of x depends on the previous.  This is the feedback nature of the non-linear equations we've been studying.

We saw last time that the value of x after many iterations (large n) may or may not converge to a given number depending on the value of c!

Play with this function again for a couple minutes putting in different values of c, especially around c<2, c=3, c=3.5, c>4.

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

quadratic[0.3, 3.56, 100];
ListPlot[xlist, Prolog -> AbsolutePointSize[3]]

For large n-values, what happens to the number of solutions as the c-value is increased?  

	(answer here)

Compare the results of c=3.52 and c=3.56.  What happend to the number of solutions?

	(answer here)
	
It would be interesting to plot (Subscript[x, n])at a big value of n, and see how this varies versus c.  This is called a bifurcation diagram.  To do this we need to redine our function so that we don't create a table of values, but only the last one for a given value of n.  Check it out.

\(quadone[xo_, c_, N_] := \((
		x = xo; 
	
        Do[x = c*x*\((1 - x)\), {N}]; x
	)\)\)

Put in some values and see what the result is.

So, you just get the value of (Subscript[x, n]), not a list of all the values leading up to it.  Let's now plot (Subscript[x, n])for lots of different c-values.

ListPlot[Table[quadone[ .3, c, 400], {c, 1, 4,  .003}], 
      Prolog -> AbsolutePointSize[1]]

Pretty cool, but there are jumps in the solution.  It looks like we're missing something.  

What's going on?  (It may be helpful to look back at our original quadratic[] function. Plot quadratic[] again for c=3.52 and describe how the solution varies in the graph as you look at increasing n-values.)

	(answer here)

So, to make a complete bifurcation diagram, we need to overlay a bunch of plots, each with n=n+1, compared to the previous.  Let's zoom in a little (shorten the range of c) and get started making the plots.  Notice that the horizontal axis is not plotting c-values but just the number in the list.

p1 = 
      ListPlot[Table[quadone[ .3, c, 400], {c, 2.8, 4,  .003}], 
        Prolog -> AbsolutePointSize[1]]

p2 = 
      ListPlot[Table[quadone[ .3, c, 401], {c, 2.8, 4,  .003}], 
        Prolog -> AbsolutePointSize[1]]

How and why are the graphs different?

	(answer here)

Continue this process, making plots of the next (Subscript[x, n]).....Use the Show[] command to overlay the plots (I think four to eight plots is enough to see what's going on -- remember to name them, e.g., p1, p2, p3, p4 so you can use Show[p1,p2,p3,p4] to overlay them.)  Your result should look like Fig 4-22 of the text. 

To see more detail, just keep increasing N  by 1 and creating a new plot.

How many bifurcations can you clearly see?  How many plots would you need to see 2 bifurcations? 3 bifurcations? 6 bifurcations?

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