February 2, 2017 --- Class 8 --- Extrapolation to zero height, Adaptive
Step Size Method, Graphing velocity ratio vs. initial height
Activities:
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
I briefly explained the idea behind an adaptive step size
algorithm. We use the Euler or Euler-Cromer method with
step sizes dt, dt/2, dt/4. Knowing the algorithm is first
order, we can extrapolate to zero step size by taking
A12 = 2 R(dt/2) - R(dt)
A24 = 2 R(dt/4) - R(dt/2)
If these values don't agree to a desired precision, that is
because higher order corrections are to big, so we decide we
need to discard dt and use dt/2, dt/4, and dt/8 instead.
You should look at my 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.
In our next class, we will run the code 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, you may find that the step size increases or decreases
depending upon whether you pick too large or too small
an initial step size. Sometimes the print time prevents
the step size from increasing because there would not a
integer number of steps between print times.
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.
There are a few variations of my code and the line with the
velocity at y=0 may take slightly different forms. The next
line has the velocity based on constant acceleration g. You
may use slightly different awk commands, if you are not using
a code with the "new way" indication for the nonlinear
interpolation. However, it is important for you to understand
and realize how useful NF, $NF, and NR can be in awk.
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 <