January 31, 2017 --- Class 7 --- Fun with Pointers, A Body Falling in
a Position Dependent Potential, Error of the Euler-Richardson Algorithm
Activities:
Fun with Pointers (This was not covered in class. If desired we can
cover in a future class.)
We looked at a code I prepared before class
~sg/play_with_pointer.c .
It uses some features that may not be familiar to everyone.
sizeof() is a function to report how many bytes an object
requires. We defined variables of several types and found
that integers and floats take four bytes.
Doubles take 8 bytes, as do pointers. The length of a pointer
can vary according to the operating system (OS). For a 32-bit
OS, a pointer would only require 4 bytes.
Incrementing pointers increases the address by the correct
number of bytes to point to the next object the pointer can
point to. That is incrementing a pointer to an integer will
increase the address by four bytes and incrementing a pointer
a double will increase the address by eight bytes.
The code also defines two increment routines, one that has a
pointer argument and one that has an integer. Calling the latter
does not change the value of the argument in the main program,
but calling the former allows us to change the value that that
pointer argument points to.
Malloc can be used to allocate memory to store an array.
Discussion of an object falling in a position dependent potential
We wish to solve Problem 3.7 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 will be described in the next class. It 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.
Code
We looked at the code for the Euler algorithm in some detail.
We also saw the small changes required to implement the
Euler-Richardson algorithm. We compared energy conservation
with a fixed initial height. We really should automate the
process of running the code.
Aside on shell scripting
We spent quite a bit of time writing shell scripts to
run the different algorithms with various step sizes. Script
writing is a very important skill and this is something that
everyone should master. You may make your own decision about
which shell you will use. I use tcsh because it is what I
learned first. There are some people who denigrate tcsh and
think one must use sh or bash. I will create a link on the
class page to an article about this. Perhaps it will help
you decide which shell to become most familiar with.
Having a shell script lets you run multiple cases quickly. If
you master this skill it can save you a lot of time and a lot
of typing. Furthermore, you can try a new procedure
interactively, use history to see what commands you gave and
use that as a guide to create a shell script to automate
the steps you want to repeat. Having the script will remind
you, months or years later, how you did the analysis. (On the
other hand, you can also redirect your history command into
a README or other file to maintain a record of the steps you
took.
Below is an example of a script ~sg/chap3/runem.csh
#!/bin/csh
set height=100000
set alg = euler
set alg = euler_richardson
foreach step ( 2.0 1.0 0.5 0.2 0.1 0.05 0.02 0.01 0.002 0.005)
$alg <${alg}_p4_$step
$height
$step
4.0
EOF
end
With the shell script, we were able to run multiple step
sizes and euler, euler_cromer and euler_richardson algorithms.
Here is a similar script runem.bash in bash syntax:
#!/bin/bash
height=100000
alg=euler_richardson
alg=euler
for step in 2.0 1.0 0.5 0.2 0.1 0.05 0.02 0.01 0.002 0.005 0.001 0.0005 0.0002 0.0001
do
./$alg <${alg}_p4_$step.frombash
$height
$step
4.0
EOF
done
Note that the spaces have been removed before and after th
equals sign in the definitions of alg.
To plot the result as a function of step size, we created
another shell script.
Here is an example, ~sg/chap3/plotem.csh
#!/bin/csh
set noclobber
set height=100000
#set alg = euler_cromer
set alg = euler
#set alg = euler_richardson
rm -f temp1.ax temp2.ax
foreach step ( 0.2 0.1 0.05 0.02 0.01 0.002 0.005)
echo -n $step " " >>! temp1.ax
awk '$1==140{print $2}' ${alg}_$step >>temp1.ax
end
axis lx step ly position >! temp2.ax
awk '$1==140{print $3}' ${alg}_$step >>temp2.ax
end
axis >! temp1.ax
awk '$1==140{print $2}' ${alg}_p4_$step >>temp1.ax
end
axis r 0.1 h .9 w .9 lt $alg lx step ly position >! temp2.ax
awk '$1==140{print $3}' ${alg}_p4_$step >>temp2.ax
end
axis >! temp3.ax
awk '$1==140{print $5}' ${alg}_p4_$step >>temp3.ax
end
axis