Sept. 8, 2005 --- Class 5---Error Analysis for the Euler Method, Shell Scripts,
Introduction to Chapter 3
Activities:
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 compilied 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.
The shell script ~sg/chap2/ans_vs_step.csh which solves
this problem. Here is ans_vs_step.csh with extra comments:
#!/bin/csh # anything after the hash mark or sharp is a comment
# the first comment is special in that it tells the shell
# to use /bin/csh to interpret or run the commands in this file
# such a comment starts with #!
# the word foreach starts a loop. The loop ends with the "end" line
foreach step ( 0.5 0.2 .1 .05 .02 .01 .005 )
echo -n $step " " # -n tells echo not to end with a new line
euler $step |tail -1|awk '{print $2}'
end
The script above with use the tail command to look at the last
line of output. There is also a head command to look at the
beginning of a file or standard input. head -9 will look at the
first 9 lines. head +3 will look at all of the file starting
with line 3. Consult the man pages for head and tail
to learn more.
Plot the result as a function of step size.
Error estimates
Most of 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!
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.