October 19, 2006---Class 16--- AHO, NDSolve, Chaos and Logistic Map

Activities:
	
	Anharmonic Oscillator: Analytic Considerations

	For homework you were asked to study the period of the 
	harmonic oscillator and how it depends on the force constant.
	It turns out that a lot can be learned about this problem by
	analytic means.  We discussed how to use time translation to
	get simple initial conditions.  We then considered how we
	could rescale x and t to simplify the equation of motion.
	When we had done so, we found that we could determine the
	dependence of the period on k, m and the initial amplitude,
	up to a single unknown constant.  This constant was the period
	with unit mass, k and amplitude.  We just found that period using
	Mathematica.

	Solving Differential Equations with Mathematica: Anharmonic Oscillator

	If you have a simple differential equation to solve, you might want
	to use Mathematica rather than a numerical code you develop yourself.
	To solve a differential equation numerically, we use NDSolve.  To solve
	one analytically, use DSolve.  For example:

	NDSolve[ {x''[t] == - x[t]^3, x[0] == 1, x'[0] == 0}, x, {t, 0, 10}]

	is the command required to solve the anharmonic oscillator with k/m=1.
	x''[t] stands for the second derivative at time t.  x'[0] stands for
	the first derivative, or velocity, at t=0.  Within the list that is
	the first argument to NDSolve, we have the differential equation
	itself and the initial conditions on position and velocity that are
	required to specify the solution.  The solution will be an 
	"interpolating function" for x[t] over the range
	0 <= t <= 10.

	If your next command is

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

	Mathematica will display a plot of the solution.  The command
	Evaluate[x[t] /. %] tells Mathematica to use the output of the 
	previous command as a substitution rule and then to evaluate the
	function x.  Assuming that In[1] was the command to solve the
	differntial equation, Out[1] is the rule for the interpolating function
	that defines the solution.  You can then find the solution at a 
	particular time with a command like:

	x[7.5] /. Out[1]

	I do not find the following commands intuitive, but they can be used
	to find the period:

	x /. First[Out[1]]

	FindRoot[ %[t] == 1, {t, 7.5}]

	The output of the first command is the interpolating function itself.
	In the second command, % stands for the interpolating function, we
	make it a function of t and then ask Mathematica to find the value of
	t that makes the function 1.  We have the search start near 7.5 because
	looking at the plot of the solution to the anharmonic oscillator, we
	see the the period is close to 7.5.  The result from Mathematica is
	7.4163.

	Chaos and Logistic Map Continued

	We continued looking at our population model, the logistic map.

	When r is <0.25, no matter what the starting
	population, after a long time the population approaches 0.
	
	For larger values of r, the long time behavior is a constant and
	the final value increases as r increases.

	With r=0.78, we find something very curious.  The
	population oscillates in time between two different values
	5.475744e-01 and 7.729384e-01.  You should verify that this does
	not depend on the starting value (as long as it is not 0 or 1).
	This result is called having an attractor of period 2.

	Now, we are getting the idea that there can be different long time
	behaviors.  For small r, we go to zero; for larger r, we approach a
	non-zero constant; for yet larger r, we oscillate between two
	values.

	With r=0.863, we found that the system has a periodicity of 4.
	With even larger values of r, there are more period doublings until
	the system finally becomes chaotic.

	We did a fixed point analysis to find that if |f'(x_f)| < 1, we
	have a stable fixed point, otherwise the fixed point is unstable.
	It was easy to find the fixed points from the equation

	 	x_f = 4 r x_f (1-x_f)

	There are two solutions, x_f = 0 and x_f = 1- 1/(4r)

	We easily calculated the derivatives at the fixed points.  For r <
	0.25, x_f=0 is a stable fixed point, and the other value is not
	between 0 and 1.  For 0.25 < r < 0.75, x_f = 1- 1/(4r) is a stable
	fixed point and x_f = 0 is unstable.  When r > 0.75 neither is a
	stable fixed point.  However, when we have a trajectory of period
	2, then x_(n+2) = f(f(x_)) can have a stable fixed point.  We
	explored this using Mathematica.