September 13, 2007 --- Class 6 --- More shell scripts; Euler-Richardson Algorithm;
Multiple Plots; A Body Falling in a Position Dependent Potential
Activities:
More shell scripts:
I forgot to go over some material related to testing the error
from integrating dy/dx = 3 x^2. We have seen how we can test
the exact error, using an awk script.
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!
Motivation for Euler-Richardson Algorithm
We discussed how a numerical integral done via a Taylor
series expansion about the midpoint of the interval results
in an integration formula with no (dt)^2 term. This was
used to motivate the Euler-Richardson algorithm that uses
the acceleration and velocity at the midpoint of the time
interval to update velocity and position through the whole
time interval.
You will get to implement the Euler-Richardson algorithm for
homework. There is a program in the textbook.
Another use for parentheses--- multiple axis graphs on one page
In ~sg/three_graphs.csh is a shell script that uses axis to
print three graphs on a single page.
#
(awk '{print $1,$2}' fall.dat |axis h .3 lx t ly y ; \
awk '{print $1,$3}' fall.dat |axis h .3 u .33 ly v -s ; \
awk '{print $1,$4}' fall.dat |axis h .3 u .66 ly a -s ) \
| plot -T x -s
In this case there are three awk and axis commands in one set of
parentheses. The awk commands pick out different columns from
the same data file. The axis commands make each graph shorter
with the -h command and graphs 2 and 3 are raised 1/3 and 2/3
of the page height with the -u option. Note that the second
two graphs must use the -s option to save the page. Also you
must add the -s option to plot, or it will complain or start
a second plot window. (I have seen both on different systems.)
The parentheses assure that all of the output from each
axis command is sent as a single stream to plot, so there is
only one xplot window with all three graphs.
Discussion of an object falling in a position dependent potential
We wish to solve is problem 3.4 of CSM. We drop a ball
from high above the surface of the earth (neglecting air
resistance!). From how great a height must we drop it so
that it's velocity when it hits the ground is 1% different
from what it would be assuming constant acceleration.
There are several potential sources of error:
1) We know from the coffee cooling problem that when we do
numerical integration there can be a step-size error. So
we must be careful about this.
2) When we integrate the equations in time, we won't
necessarily stop integrating exactly when the object
reaches height 0. We may have to interpolate to find the
velocity when the height is zero.
3) We also must find what height results in a 1%
difference. Do we guess initial heights? How do we
accurately find the height that corresponds to 1%?
Step size error
We have three approaches to the step size error issue.
One possibility is to vary the step size and see how the
answer varies. The accuracy on the velocity at height=0
must be much greater than 1%. The accuracy will depend on
height and step size.
We can also try to find a better method than Euler or
Euler-Cromer. These are each first order methods. The
Euler-Richardson method is one order higher. Try running
this code with different step sizes. Note the accuracy of
the energy and the position. (Interestingly, the
difference is negligible for the velocity.)
Another scheme involves using the Euler-Cromer algorithm
with different step sizes to extrapolate to zero step size.
This was described on the board and is coded in
adapt_step.c. This code will increase or decrease the step
size to achieve a desired level of accuracy during the
integration.
We examined the code in adapt_step.c in detail.