October 22, 2013 --- Class 16 --- AHO, Finding Superstable Trajectories
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.
You were expected to find the period of oscillation for a number
values of k and to make a plot of period vs k. Having done that,
you were expected to determine the equation that corresponds to the
plot.
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.
Here is some material about using Mathematica was not done in class:
Solving Differential Equations with Mathematica: Anharmonic Oscillator
If you have a simple differential equation to solve, you might want
to use Mathematica rather than a numerical code you develop yourself.
To solve a differential equation numerically, we use NDSolve. To solve
one analytically, use DSolve. For example:
s=NDSolve[ {x''[t] == - x[t]^3, x[0] == 1, x'[0] == 0}, x, {t, 0, 10}]
is the command required to solve the anharmonic oscillator with k/m=1.
x''[t] stands for the second derivative at time t. x'[0] stands for
the first derivative, or velocity, at t=0. Within the list that is
the first argument to NDSolve, we have the differential equation
itself and the initial conditions on position and velocity that are
required to specify the solution. The solution s will be an
"interpolating function" for x[t] over the range
0 <= t <= 10.
If your next command is
Plot[ Evaluate[x[t] /. s], {t, 0, 10}]
Mathematica will display a plot of the solution. The command
Evaluate[x[t] /. s] tells Mathematica to use the interpolating
function as a substitution rule and then to evaluate the
function x. You can find the solution at a
particular time with a command like:
x[7.5] /. s
I do not find the following commands intuitive, but they can be used
to find the period:
x /. First[s]
FindRoot[ %[t] == 1, {t, 7.5}]
The output of the first command is the interpolating function itself.
In the second command, % stands for the interpolating function, we
make it a function of t and then ask Mathematica to find the value of
t that makes the function 1. We have the search start near 7.5 because
looking at the plot of the solution to the anharmonic oscillator, we
see the the period is close to 7.5. The result from Mathematica is
7.4163.
Superstable trajectories
In our discussion of fixed point iteration, we found that if
x_{n+1} = f(x_n)
the rate of convergence near a fixed point is dependent on the
derivative of f at the fixed point. In our case, we are interested
not only in f(x)=4r x (1-x), but in the iterated functions f^(n)(x).
Because of the symmetry of f(x) around x=1/2, the iterated functions
all have a derivative of zero at x=1/2. When one of the fixed
points is at x=1/2, we have what we call a superstable trajectory.
We used Mathematica to find the values of r that lead to
to superstable trajectories, i.e., the trajectories that include
x=0.5.
Solve[ f[1/2, r] == 1/2, r]
will solve for the trajectory of period 1. You can define the
iterated functions, f2[x_,r_] := f[f[x,r],r] and
f4[x_,r_] := f2[f2[x,r],r], etc. and similarly solve for values of
r. See ~sg/Documents/March01_2012.nb for a Mathematica notebook.
I discussed the program findroot.c that uses a routine from
Numerical Recipes to solve these equations numerically.
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. Quadratic interpolation is
used in the Brent algorithm, but is used to interpolate x as a 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.
I downloaded this file from a web site at Cornell that has the text
of Numerical Recipes, at least it did a few years ago.
Running my code, I find values of s_k different from those appearing
in the book in Table 6.2. My values are correct as you should be
able to verify by iterating the fixed point equation starting with
x=0.5 and see who gets closer to 0.5 after the appropriate number
of iterations.
Root Finding with Gnu Scientific Library
The GNU Scientific Library also has root finding routines and there is
no restriction on the use of their library. To find out about
what is available from GNU, go to:
http://www.gnu.org/software
or for the GNU Scientific Library, go here:
http://www.gnu.org/software/gsl/
Documentation for the one-dimensional root finding routines
are found here:
http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html
In ~sg/gslroots you will find code for the example in the documentation
and for my code to solve for the values of r that result in
superstable trajectories. The code is more verbose that that using
the NR routines in that for each solve details of the convergence
are printed out. When high enough accuracy is demanded, the results
are seen to be the same as found with the NR based code. As pointed
out before, my values of r are more accurate that found in Table 6.2.
This was verified in one case by using my value of r and the one in
the book. We found that after the expected period with my value
x returned to 0.5, but with the value in the book, it was slightly off.