October 26, 2004---Class 17---Jackknife Method Via Mathematica

Activities:
	
	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.

	Comments on Homework

	We considered how to use awk and Mathematica for problem 3.2.

	To create a table of numerical derivates and accelerations from
	the data in table 3.1, we prepare a file with time (>=0)
	and position as in the table.  We then create an awk file
	central_diff.a with these commands:

{t[NR]=$1; x[NR]=$2}
END {for(i=2; i<(NR-1); ++i) {v[i]=(x[i+1]-x[i-1])/(t[i+1]-t[i-1]);
a[i]=(x[i+1]-2*x[i]+x[i-1])/(t[i]-t[i-1])/(t[i]-t[i-1]); print t[i],x[i],v[i],a[i];}}

	The first line reads the data from each line of the input file and
	stores vectors of time and position.  "END" matches the end of
	the input file and then implements the finite difference formulae.
	The result is a table with time, position, velocity and acceleration
	on each line.

	We would like to determine if a linear or quadratic dependence of
	the drag force is more appropriate.  We should have either
				|a| = g ( 1 - v/v_t)
	or
				|a| = g ( 1 - (v/v_t)^2 ) .
	If we plot a vs. v, we find the data is quite noisy.  We can use
	Mathematica to do a least squares fit to the data.

Prepare a Table of values of velocity and acceleration:
data = { {2.25, 8}, {3.05, 8}, {3.725, 5.5}, {4.3, 6}, {4.8, 4}, {5.2, 4} }

Use Mathematica's Fit routine to fit to a linear or quadratic function of v:

                                                                                
In[15]:= Fit[ data, {1,x}, x ]
 
Out[15]= 11.8219 - 1.51903 x
 
In[16]:= Fit[ data, {1,x^2}, x ]
 
                             2
Out[16]= 9.18996 - 0.202868 x

	Unfortunately, we have no goodness of fit, but we see that the
	quadratic fit is much closer to having the constant part be 9.8
	as the acceleration should approach g for v=0.  So, it would appear
	that quadratic is favored by this analysis.  Of course, fitting the
	trajectory and adjusting v_t for both force laws is also a good
	idea and can be done for the second part of the problem.

	Another approach is to subtract 9.8 m/s^2 from the acceleration and
	then fit using a single function x or x^2.

Do[ data[[i,2]]=data[[i,2]]-9.8, {i,1,6}]

In[3]:= data

Out[3]= {{2.25, -1.8}, {3.05, -1.8}, {3.725, -4.3}, {4.3, -3.8}, {4.8, -5.8}, 
 
>    {5.2, -5.8}}

In[4]:= Fit[ data, {x}, x ]

Out[4]= -1.03189 x

In[5]:= Fit[data, {x^2}, x ]

                   2
Out[5]= -0.233811 x

In[8]:= g=9.8

Out[8]= 9.8

In[9]:= Solve[g/vt==1.03189, vt]

Out[9]= {{vt -> 9.49714}}
 
In[10]:= Solve[g/vt^2==0.233811,vt]

Out[10]= {{vt -> 6.47412}, {vt -> -6.47412}}

	Another fine approach to finding the value of v_t is to fit the 
	observations of the position as a function of time by adjusting v_t
	to most closely fit the last time point available.


	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.

	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.