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!