September 19, 2013 --- Class 8 --- Error of the Euler-Richardson Algorithm,
Extrapolation to zero height, Graphing the velocity ratio
Activities:
Error of the Euler-Richardson Algorithm
To find the dependence of the error in the Euler-Richardson
we used the output created by the script runem.csh that
has step sizes from 0.005 to 2.0 seconds.
We plotted the position at t=140 for various step sizes.
The dependence seemed nonlinear, but we really wanted to
find the error.
We then discussed how we can easily use awk to find the
error. plotem2.csh or plotem3.csh leave a file temp1.ax
with the position as a function of step size. It looks
like this:
[sw246-05:~/chap3] sg% mo temp1.ax
2.0 6.454932e+03
1.0 6.454655e+03
0.5 6.454585e+03
0.2 6.454566e+03
0.1 6.454563e+03
0.05 6.454562e+03
0.02 6.454562e+03
0.01 6.454562e+03
0.002 6.454562e+03
0.005 6.454562e+03
Clearly, given the limited number of digits we printed out,
the error is zero by the smallest step size, so we can use the
postion on the last line as the correct value and subtract
that from the previous values to get the error. How do
we get awk to do the subtraction on lines with larger step
size? Any easy way is to store what we read in and do the
subtraction at the end. This is complicated enough that it
deserves to be in a file of its own.
Here is ~sg/chap3/subtract_last.a:
{step[NR]=$1; y[NR]=$2}
END{const=y[NR]; for(i=1; i<=NR;++i) print step[i],y[i]-const}
The result of running this awk script on temp1.ax
is this output
2.0 0.37
1.0 0.093
0.5 0.023
0.2 0.004
0.1 0.001
0.05 0
0.02 0
0.01 0
0.002 0
0.005 0
To see if the error obeys a power law, we would like to make
a log-log plot and see it is a straight line. However,
we can't do that with the data above, because no computer
likes to calculate the log of zero. However, if we use the
head command, we can restrict the input of axis to the first
5 lines.
awk -f subtract_last.a temp1.ax |head -5 | axis x l y l |xplot
does the trick! You will see a perfectly straight line that
indicates that Euler-Richardson errors are order (Delta t)^2.
This is really great, because if means and if we reduce
the step size by two, the error is four times smaller.
That is a big saving.
Extrapolation error
The code integrated until the height is negative, but we need
the velocity when the height is zero. When I first wrote the
code, I used a constant velocity extrapolation to get back to zero
height. However, that is dumb and when print size was a
second or more could result in a significant error. It is
much smarter to recall elementary mechanics for uniform
acceleration. In this way, we compute the speed at
zero height as sqrt(v*v+2*g*y). See ~sg/p609_interpolation.pdf
for the details.
Adaptive step size algorithm
We looked at the code ~sg/chap3/adapt_step_nonlin_interp_2012.c
This code can change its step size to keep the velocity error
small. We don't have to be too careful about the step
size we set since the code will adjust.
We ran the code adapt_step_nonlin_interp_2012 to see how
it increases or decreases the step size depending on how
accurately our two estimates of the velocity agree with
each other. With different initial step sizes and print
times, we found that the step size could increase or decrease
depending upon whether we picked too large or too small
an initial step size. Sometimes the print time prevented
the step size from increasing because there would not a
integer number of steps between print times.
Root finding
Since we really want to find the initial height that
gives v_true/v_naive=0.99, we are really trying to
find the root of an equation. Let ratio(height)=
v_true/v_naive. We are trying to solve ratio(height)-0.99=0.
In the program adapt.c, the secant method is used to solve
the problem of what height must the object be dropped from.
We will talk more about root finding in the next class.
Graphing the ratio to the naive velocity
The problem is stated in terms of the ratio of the velocity in
the position dependent potential as compared with constant
acceleration g. Of course, we can find the velocity with
constant acceleration, so the code currently prints that out
on a line by itself after the velocity from integration with
the position dependent acceleration is printed. We may use
AWK to get the ratio of velocities.
We talked about how we can find the two correct values.
We use the awk file ~sg/chap3/ratio.a. It contains two lines:
/new way/{record=NR; v=$NF}
NR==(record+1) && NF==1{print v/$1}
The first line looks for the output line from
adapt_step_nonlin_interp_2012 that looks like this:
Velocity at y=0 new way: 1.389139e+03
This is our new interpolation formula for the speed when y=0.
The following line of output has a single value that is the
naive velocity, so we have awk look for the next line number
assure that there is only one value on the line. We
print the ratio of the two values.
Graphing the velocity ratio as function of the height via shell script
It is relatively easy to make a graph of the velocity ratio
as a function of the initial height. We can use a foreach
loop to set the initial height. Then, we need to manipulate
the output of adapt_step_nonlin_interp_2012 to get the
velocity ratio. This last step is done with the awk file
described above.
Below is a shell script to run adapt_step_nonlin_interp_2012
and pipe the output through awk. If you run this at the
screen, you will find all the prompts from
adapt_step_nonlin_interp_2012, which are written
to the standard error, mixed up with the output from the
awk script, which are written to the standard output.
If you redirect the standard output to a file with the ">"
symbol, only the output you want will wind up in the file.
You may also pipe the output of the script directly to
axis and xplot. The script is ~sg/chap3/ratio_vs_height.csh .
#!/bin/csh
foreach height ( 20000 40000 60000 80000 100000 120000 140000 160000)
echo -n $height " "
adapt_step_nonlin_interp_2012 <