Sept. 9, 2004 --- Class 4 --- Error Analysis for the Euler Method, Shell Scripts

Activities:
	
	First homework:  The first homework, which is just a simple exercise
		was assigned and due on Tuesday.  See the class web page
		Homework section for details.

	A simple integral first (CSM, Secs. 2.2-2.5)
		In ~sg/chap2 you will find a file euler.c that you can use
		to do a numerical integral.  It might be set up with the
		function 2*x or 3*x^2.

		We complied the euler program using the command:
		cc -o euler euler.c
		This creates an output executable called euler.  If you just
		type cc euler.c, your executable will be called a.out.

		We ran the program for different steps sizes and found that
		the final answer (and all intermediate values) vary with
		step size.  However, the value value converges as 
		the step size decreases.

		We would like to automate this procedure so that we can easily
		plot the final result as we vary the step size.  We looked
		at the shell script ~sg/chap2/ans_vs_step.csh which solves
		this problem.

		Plot the result as a function of step size.


	Error estimates
		The rest 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.

		We wrote an awk script to show that the error for the
		function f(x)=3*x^2 agrees perfectly with the expression above.
		We would like to combine the two steps of running 
		ans_vs_step.csh and running awk on the output into a single
		script.  Note that a simple solution is to pipe the output of
		ans_vs_step.csh to the desired awk command.

		However, if one wants a single shell script there is the
		challenge of dealing with the fact that ans_vs_step.csh
		used two separate commands, echo -n and euler to print
		a single line.  Here are two ways to solve the problem.

		Append output to a temporary file and then awk the file.
		This is found in ~sg/chap2/err_vs_step.csh

		#!/bin/csh
		set noclobber
		rm -f temp.ans
		foreach step ( 0.5 0.2 .1 .05 .02 .01 .005 )
        		echo -n $step "  "  >>! temp.ans
        		euler $step |tail -1|awk '{print $2}'  >>temp.ans
		end
		awk '{print $1,(8-$2),4.5*$1,4.5*$1-0.5*$1*$1}' temp.ans

		Notes: >> appends the standard output to the filename that
		follows.  >>! does it even when the file does not already
		exist.  This is necessary if you set the noclobber variable.
		In this example we first remove the temporary file temp.ans
		with the rm -f command.  The -f option forces the file to
		be removed without displaying permissions, asking questions or
		reporting errors.  So, if the file does not already exist,
		you will not have an error message.  After the end of the
		foreach loop, we run the single awk command to create the
		display we want.  You must run this command in a directory
		where you can create the file temp.ans.

		An alternative solution uses the csh () to join two commands.
		It can be found in ~sg/chap2/err2_vs_step.csh

		#!/bin/csh
		set noclobber
		foreach step ( 0.5 0.2 .1 .05 .02 .01 .005 )
       		(echo -n $step "  "; euler $step |tail -1|awk '{print $2}') | \
          		awk '{print $1,(8-$2),4.5*$1,4.5*$1-0.5*$1*$1}' 
		end

		Notes: In this case, we use two unix commands that both
		print to the standard output and we put the two commands
		within parentheses.  The two commands are separated by a
		semicolon.  We then pipe that single output 
		stream to the awk command.  It is on a separate line, but the
		line that starts with "(echo" ends with a backslash.  That tells
		the shell that the next line is a continuation.  Important:
		don't put any spaces after that backslash!