September 11, 2007 --- Class 5---Error Analysis for the Euler Method (II),
		Introduction to Chapter 3, Energy (non)conservation

Information:

Activities:

	Error estimates
		Most of of this class was devoted to understanding the
		error in the Euler algorithm for integrating a differential
		equation.  We did a careful estimate of the error for the
		simple case 
			dy/dx = f(x).

		This is simple because the right hand side depends on x,
		not y(x).  So this is just doing an integral.

		If we use a step size dx, the error in a bin with left hand
		edge at x is:

	(1/2!) (dx)^2 f'(x) + (1/3!) (dx)^3 f''(x) + higher order terms

		In our first example of the Euler program, f(x)=2x.  Thus,
		f'(x)=2 and f''(x)=0, so the second term above is 0.  Why
		is the error order dx, not (dx)^2?  That is because the
		above expression is for each bin.  In our example, we
		integrate x from 1 to 2, so the number of bins n=(1/dx).
		Thus, the total error is dx for this case.

		In our second example, f(x)=3x^2. f'(x)=6x and f''(x)=6,
		so the second error term does not vanish.  If you very
		carefully evalute the two error terms you will find that
		the first term is 
			3 (dx)^2 * n (1 + (n-1)*dx/2 )
		where n=1/dx.  Using the fact that n*dx=1, this first
		term can be written as
			3 dx (3/2 -dx/2) or 9/2 *dx - 3/2 * (dx)^2.
		Actually, only this term requires care.  The f'' term is
		easily seen to be 
			(dx)^2.
		Adding the two terms, our final result is
			9/2 *dx - 1/2 * (dx)^2.

		Before we did the careful estimate, we plotted the error
		divided by dt and found that the coefficient of the second
		term is negative, not positive and we would get from just
		looking at the second term in the Taylor expansion.

	Introduction to Chapter 3

	Solving Newton's Law of Motion

		We started Chapter 3 and discussed how to turn a 2nd order
		differential equation into a series of coupled first order
		differential equations.  We also discussed why position,
		velocity and acceleration are so important, and the next
		derivative, called the jerk, is not well known.

		The program fall_euler_ansi.c in ~sg contains code to 
		integrate the equations for a falling body.  
		It is set up for the Euler algorithm.  We discussed the 
		code in detail.

		If the lines for updating v and y are
		reversed, you get the Euler-Cromer algorithm.  
		You will apply the Euler-Richardson algorithm, and 
		explore the step size dependence of the trajectory of a
		falling body.  We talked about energy conservation and
		how to write an awk expression to calculate the energy.  

	Energy (non)conservation

		We studied the energy of a falling body as caculated by the
		Euler algorithm.  We looked at the energy at t = 1 sec.
		We also created a shell script to automate the procedure:

#!/bin/csh
# this is the best shell script I have ever written
foreach step ( 0.025 0.05 0.1 0.2)
fall_euler_ansi <<EOF | awk '{if($1 == 1.0) print '$step',0.5*$3*$3+9.8*$2}'
    $step
    .2
EOF
end

		Notes: When a command expects its input from the standard input
		but you want to run it from a shell script, the << sign is
		used to redirect the standard input so that it comes from the
		shell script.  To signal the end of the redirected input,
		we need a string.  By convention, I always use EOF, which is
		short for End Of File.  Note that EOF must occur right at the
		beginning of the line.  Don't indent it!
		In this case both the input and output of the fall command are
		redirected.  The output is being piped to an awk command that
		prints the step size and the energy if the first entry on the
		line is 1.0.  Note use of == to test for equality.  This is
		is test of numerical equality.  It is not a test of string 
		equivalence.  Also note that there are now two sets of
		matching single quotes in the awk command.  This means that
		the shell variable $step is not within the single quote.
		Thus, the shell does variable substition and replaces $step
		by its current value determined by the foreach loop.  If you
		put in any extra spaces around $step, the awk command will
		not work because awk takes its command from a single
		string, and the extra space breaks the string that appears
		above into two strings.  These things can be tricky, so if
		you have problems consult me.  I have said it before, but I
		will repeat myself.  The reason I like to meet in a computer 
		cluster is because I know from years of experience that very 
		simple typing changes on your part, can make a big difference.  
		I can sometimes save you hours of pondering what is wrong 
		by spotting a very tiny error.  Don't hesistate to email me 
		if you are having trouble outside of class.

		Getting back to the numerics, we found that for fixed time, the
		energy is a linear function of the step size.  Of course, the
		energy should be constant.  We also found for a single
		step size that the energy is linearly increasing with time.

		It was very easy to modify the C code to change to the 
		Euler-Cromer algorithm.  It was then trivial to rerun the
		shell script to produce the data for this algorithm.  I edited
		the data file to include labels for axis.  We found that both
		algorithms converge to the right value as the step size
		approaches 0, but their errors have the opposite sign.
		Here is the data.  (Euler-Cromer follows Euler).

# lx "step size" ly Energy
# cm 2
.025 173 "Euler"
# cm 3
.025 175 "Euler-Cromer"
# cm 2
0.025 197.201
0.05 198.401
0.1 200.802
0.2 205.604  " t=1"
# cm 3
0.025 194.8
0.05 193.599
0.1 191.198
0.2 186.396  " t=1"
# cm 2
0.025 198.401
0.05 200.802
0.1 205.604
0.2 215.208  " t=2"
# cm 3
0.025 193.599
0.05 191.198
0.1 186.396
0.2 176.792  " t=2"

		The Euler-Richardson algorithm, which first updates x and v
		to the middle of the time interval and then recalculates the
		acceleration, is a much better algorithm.  We will talk more
		about it in the next class.