October 17, 2006---Class 15---Jackknife Mean, Jackknife Method Via Mathematica,
				Introduction to Chaos and Logistic Map

Activities:
	
	Jackknife Mean
	
	In addition to using the jackknife method for bias reduction, you
	can use it to calculate the sample mean.  The sample mean is
	related to variance of the jackknife samples.  Because these
	samples are so correlated (only one variable changes in each
	sample), the variance of the jackknife values must be multiplied by
	(n-1):

	     (\sigma_J)^2 = (n-1)/n \sum (S'_i - S)^2

	where S'_i is the jackknife statistic with the i'th measurement
	removed and S is the statistic on the entire set of n measurements.

	Notice that in the simple case that the statistic is the mean value
	of the sample, the naive variance would be given by:

	\sigma^2 = 1/[n(n-1)] \sum (x_i - x_avg)^2

	In this case, it is easy to show S'_i - S = (x_avg - x_i)/(n-1).
	Here is the deriviation (I did not complete in class):

	(n-1) S'_i + x_i  = sum of all x_i = n(S) = n x_avg = (n-1)S + x_avg

	since S and x_avg are the same in this example.  Bring the terms
	involving S to the left, the terms involving x to the right and
	divide by (n-1) to complete the derivation.

	Jackknife Method Via Mathematica

	Mathematica can be used to generate samples of random numbers and
	to apply the jackknife technique.  This is a good opportunity to
	introduce a new mathematica construct the "module", which is like
	a subroutine.  In ~sg/jackn.math  is the code to implement
	this method.

	First we need a way to create our set of random variables.  This is
	done via the newvars command:

	newvars[num_] := Do[z[i] = Random[], {i, num}]

	Once this is defined, you may create as many new random numbers
	in the vector z by calling newvars with the number you desire
	as the argument.

	A Module in mathematica takes two arguments.  The first is a list
	of variables local to the module, and the second is an expression,
	i.e., a sequence of statements.

	Using a module statement we define jackknife[x_, num_Integer],
	where x is vector of random numbers and num is the sample size.
	Here is the code:

jackknife[x_, num_Integer] := 
  Module[{mysum, mysumpr, stat, statpr, estimate, fctr, i}, 
   fctr = N[(num - 1)/num]; mysum = Sum[x[i], {i, 1, num}]; 
    stat = num/mysum; estimate = num*stat; 
    Do[mysumpr = mysum - x[i]; statpr = (num - 1)/mysumpr; 
      estimate -= fctr*statpr, {i, 1, num}]; {stat, estimate}]

	Successive calls to newvars and jackknife are equivalent to what
	was done in the C code in the last class.

October 28, 2004 --- Class 18 --- Energy Conservation; Logistic Map and Chaos

	Chaos and the Logistic Map

	We introduced the logistic map from 
	consideration of an insect population model.  

		P_{n+1} = P_n (A - B P_n)

	Upon scaling the population, the population is
	replaced by x which lies in the range [0,1].  The next generation
	is given by

	x <- 4 r x * (1-x) .

	A program entitled map allows you to follow the time history of a
	population.  Map takes three command line arguments: r, the
	parameter that appears above in the definition of the map; x_0, the
	initial value of the population; and iterations, the total number
	of generations to calculate.

	In the next class, we shall see what the long term future of the
	population is.
November 2, 2004 --- Class 19 --- Relaxation, Runge-Kutta, Logistic 
			Map and Chaos (continued)
	
Activities:
	
	Relaxation Time

	Some people misinterpreted the definition of relation time for the
	damped harmonic oscillator.  You need to consider the peak height
	in every half-cycle of the oscillator.  It is these peak heights
	that decay exponentially.  How can we easily see this?  One way to
	to make a semi-log plot of the trajectory.  However, half the time
	the coordinate is negative, so we must change the sign of the negative
	values.  This is easily done with awk.  In a file abs_traj.a we
	have:

	{if($2 <0) print $1,-$2 else print $1,$2}

	Then if your program output has time in column 1 and position in
	column 2 and is called data:

	awk -f abs_traj.a data |axis y l |plot -T x

	will make a nice graph.  There is a series of bumps of decreasing
	amplitude and it is obvious that the peaks of the bumps fall on a
	straight line.  The slope of that line determines the relaxation
	time.  The peaks occur when the velocity changes sign.  If you
	identify the time and position where they occur,  you can easily
	solve for the relation time.

	An alternative procedure to show the exponential decay is to draw
	on the graph the "envelope" around the curve showing the
	oscillation.  If the initial amplitude is 1, the envelope is 
	+-exp(-t/tau).  We used Mathematica to make a table of values of
	the exponential.  

	data = Table[ {t,Exp[-t/4.]},{t,0,10,.2}]
	
	We exported the data to a file with this command:

	Export["nov2.dat", data, "Table"]

	After adding a blank label to the last line, we are ready to plot 
	all the curves:
(awk '{$2= -$2; print}' nov2.dat; cat nov2.dat data )|axis |plot -Tx

	The order of the files makes use of the fact the the two envelope
	curves end with blank labes and the data file does not.  Also, only
	the first two columns (time and position) are plotted from the data
	file.  You might instead use:
(awk '{$2= -$2; print}' nov2.dat; cat nov2.dat; awk -f traj.a data ) \
    |axis |plot -Tx

	Fourth Order Runge-Kutta

	A number of students did not use the 4th order Runge-Kutta
	algorithm as described in class.  They used use some strange
	notation where instead of having k1x, k1v, k2x, k2v, etc. they just
	used k1, k2, k3 and k4.  These variables were first defined in
	terms of the acceleration function.  For example:

	k1 = dt*accel(omega_sqr, (*y_ptr) );
	k2 = dt*accel(omega_sqr,(*y_ptr+0.5*k1));

	The real RK algorithm uses:

	k1v = dt*accel(omega_sqr, (*y_ptr) );
	k1x = dt*(*v_ptr);
	k2v = dt*accel(omega_sqr,(*y_ptr+0.5*k1x));
	k2x = dt*(*v_ptr + 0.5*k1v);

	Notice that in the third line for k2v we use the argument:
	*y_ptr+0.5*k1x.  On the other hand, in the k2 line the argment is
	*y_ptr+0.5*k1, where k1 is equivalent to k1v.  We can't have k1x
	in one place and k1v in the other.  Thus, the first set of two
	lines is totally bogus.  Another way of seeing this is by
	dimensional analysis.  k1 and k1v have the dimensions of a
	velocity.  So, you clearly can't add 0.5*k1 to a position because
	that is adding a lenght to a velocity.

	I ran the bogus algorithm and showed that energy conservation did
	not hold very well in comparison to the properly coded 4th order
	Runge-Kutta.  Woe to those who don't check energy conservation when
	they are itegrating the equations of motion for a conservative
	system!

	Chaos and the Logistic Map

	Previously, we introduced the logistic map from 
	consideration of an insect population model.  
	Upon scaling the population, the population is
	replaced by x which lies in the range [0,1].  The next generation
	is given by

	x <- 4 r x * (1-x) .

	A program entitled map allows you to follow the time history of a
	population.  Map takes three command line arguments: r, the
	parameter that appears above in the definition of the map; x_0, the
	initial value of the population; and iterations, the total number
	of generations to calculate.

	We have seen what happens after a long time if we use r=0.1.
	The population dies out if we start with x=0.1

	map 0.1 0.1 100

	Students were asked to try different starting values.  We then
	constructed a shell script to loop over different starting values.

	#!/bin/csh
	set r = 0.1

	foreach x (0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9)

	    map $r $x 100 >>! temp

	end

	When r is sufficiently small, no matter what the starting
	population, after a long time the population approaches 0.

	With r a little greater than 0.25, after many generations the
	population approaches a nonzero constant for any starting value 
	other than 0 or 1.

	With larger r values, we see that the final value is increasing.

	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.

	In the next class, we will explore yet larger r values and try to
	understand why this is happening.