November 1, 2005 --- Class 17 --- Project, AHO, Optimal Step Size,
Advanced Root Finding
Project
Each student should send me an email by 5 p.m. on Friday, Nov. 4
with at least a one paragraph description of the proposed project.
I will respond either by approving the project, suggesting
modifications or requesting additional information.
The project is due on Tuesday, Dec. 14 at 5 p.m., as stated in the
syllabus. You project report will be graded on this basis:
Statement 10
Eq Solved 5
Num. Method 5
Code 5
Results 10
Discussion 10
Critique 5
------------------
Total 50
Note the emphasis on the statement of the problem, analysis of results
and discussion of what has been learned. Also please note that your
report should include a critical discussion of how you might do it
better (or at least differently) and how it could be extended.
Anharmonic Oscillator: Analytic Considerations
For homework you were asked to study the period of the
harmonic oscillator and how it depends on the force constant.
It turns out that a lot can be learned about this problem by
analytic means. We discussed how to use time translation to
get simple initial conditions. We then considered how we
could rescale x and t to simplify the equation of motion.
When we had done so, we found that we could determine the
dependence of the period on k, m and the initial amplitude,
up to a single unknown constant. This constant was the period
with unit mass, k and amplitude. We found that period using
Mathematica a few classes ago when we studied the NDSolve function.
Optimal Step Size for Numerical Integration
On one of the homework problems you were asked to find the optimal
step size when integrating a differential equation. We have discussed
how the error should decrease as the step size decreases for various
algorithms. However, if the step size is too small, the error is
dominated by round-off error. The optimal step size depends on
whether your code is single or double precision. In my home directory
you may find two files: ener_vs_t_stepsize.ax and
ener_vs_t_stepsize_float2.ax. The first is for double precision and
we see the change in energy is a smooth function of time. As the
step size decreases from 10^-2 to 10^-4 the change in energy drops
as expected for the Euler-Richardson algorithm. However, for smaller
step sizes, Delta E is not longer smooth, and it actually starts
to increase when the step size falls below 10^-6. So, the optimal
is probably 10^-5. For a single precision code, the optimal step
size is about 10^-3.
To find the step size that conserves energy to 0.1 % as requested
in problem 5.1d, it is useful to examine energy conservation for
various step sizes and plot at horizontal line at the desired value.
File DeltaE_vs_step.ax, when ploted as a log-log plot is quite
informative. From it you can see what step size is required for the
two values of omega you were asked to study.
Superstable fixed points
In the last class, we talked about superstable fixed points.
That means that x = 0.5 is a fixed point of
f^(n) (x,r) = x .
We want to find out which values of r have 0.5 as a fixed point.
We we are trying to find the root of this equation:
f^(n) (0.5, r) = 0.5
The program findroot.c uses routines from Numerical Recipes to
accomplish this. The routine is called zbrent and may be found in
Chapter 9 on Root Finding. This may be used on the homework on
problem 6.6. We discussed how to create a polynomial interpolating
function to go through n data points. Quadratic interpolation is
used in the Brent algorithm, but is used to interpolate x and function
of y. The desired x value is the one for which y=0. This results in
a complicated, but linear function, for x the next guess at the root.
The zbrent routine is useful when you can bracket the root and only
want to evaluate the function whose root you want, not its derivative.
You can read about this routine by looking in the file ~sg/ch9-3.pdf.
We downloaded this file from a web site at Cornell that has the text
of Numerical Recipes. You may easily find this site using Google.