February 28, 2017---Class 15---Chaos and Logistic Map, Fixed Points,
Finding Bifurcation Points Finding Superstable Trajectories, GSL
Activities:
Chaos and the Logistic Map
The logistic map is define by:
x <- 4 r x * (1-x) .
Recall, a program ~sg/map.c allows you to follow the time history of a
population. Map takes three command line arguments: r, the
parameter that appears above in the definition of the map; x_0, the
initial value of the population; and iterations, the total number
of generations to calculate. The executable is ~sg/bin/map.
In the last class, we saw what happens after a long time if we use r=0.1.
The population dies out if we start with x=0.1
map 0.1 0.1 100
When r is sufficiently small, no matter what the starting
population, after a long time the population approaches 0.
With r a little greater than 0.25, after many generations the
population approaches a nonzero constant for any starting value
other than 0 or 1.
With larger r values, we see that the final value is increasing.
With r=0.78, we find something very curious. The
population oscillates in time between two different values
5.475744e-01 and 7.729384e-01. You should verify that this does
not depend on the starting value (as long as it is not 0 or 1).
This result is called having an attractor of period 2.
Now, we are getting the idea that there can be different long time
behaviors. For small r, we go to zero; for larger r, we approach a
non-zero constant; for yet larger r, we oscillate between two
values.
With r=0.863, we found that the system has a periodicity of 4.
We would like to explore the long time behavior by looking at the
set of x values on each trajectory. To do this, we want to print
out a certain number of values on the trajectory starting after
a number of generations. A program called fixed_pts.c can do this.
It will do this for a range of r values. It arguments are
rmin rmax dr ntransient nprint
rmin is the minimum value of r to be considered, rmax is the maximum
and dr is the increment in r. ntransient is the number of initial
generations before we start printing x_n. nprint is the number of
generations to print after the transient generations.
With this program we were able to produce pictures of period doubling
and chaos as in Figure 6.2 on page 145.
We would like to better understand how period doubling comes about.
Fixed Points
We are looking at the mapping
x_{n+1} = f(x_n)
and when r <= 0.75 we find that the long term behavior is a
constant. We call this a fixed point. Let
x_{n} = x_f + \epsilon_n
where x_f is the fixed point and \epsilon_n is how far we are from
the fixed point at the n th generation.
By Taylor expanding around x_f, we found that if |f'(x_f)| < 1, we
have a stable fixed point, otherwise the fixed point is unstable.
It was easy to find the fixed points from the equation
x_f = 4 r x_f (1-x_f)
There are two solutions, x_f = 0 and x_f = 1- 1/(4r)
We easily calculated the derivatives at the fixed points. For r <
0.25, x_f=0 is a stable fixed point, and the other value is not
between 0 and 1. For 0.25 < r < 0.75, x_f = 1- 1/(4r) is a stable
fixed point and x_f = 0 is unstable. When r > 0.75 neither is a
stable fixed point. However, when we have a trajectory of period
2, then x_(n+2) = f(f(x_)) can have a stable fixed point. We
will explore this in the next class.
Finding Bifurcation Points
We are interested in accurately finding the bifurcation points
where we go from stable attractors of period 2^n to stable
attractors of period 2^(n+1). This can be done using the code
in fixed_pts.c in my home directory, which is compiled and in my
bin directory. Recall that we used this code to prepare a
diagram like Fig. 6.2 on page 134 of CSM. We also know that at
the bifurcation points, the convergence is very slow. We can
use this to iteratively run the code with narrower windows and
more initial iterations discarded to find the bifurcation
points.
For instance, if you run:
fixed_pts .7 .999 0.001 1000 128 |axis m 0 |xplot
you will get a nice diagram like Fig 6.2. If you look at the
output of fixed_pts, you will see that for each value of r,
there are 128 lines with the x values that result after 1000
iterations are discarded. Away from the r values that
correspond to bifurcation points, convergence is fast and many
of the x values are repeated as we loop over the points in
stable attractors of period smaller than 128. However, near the
bifurcation points, convergence can be slow and we will find
more than the expected number of unique x values. How can we
use this to find the bifurcation points?
Unix has a utility called sort that is very useful. The option
-u, which stands for unique, will throw out repeated lines. For
example, in the output of
fixed_pts .7 .999 0.001 1000 128 | sort -u | more
you will find output like this:
7.430000e-01 6.635262e-01
7.440000e-01 6.639785e-01
7.450000e-01 6.644295e-01
7.460000e-01 6.648794e-01
7.470000e-01 6.653278e-01
7.470000e-01 6.653279e-01
7.470000e-01 6.653280e-01
7.470000e-01 6.653281e-01
7.470000e-01 6.653282e-01
7.480000e-01 6.657665e-01
7.480000e-01 6.657666e-01
7.480000e-01 6.657667e-01
7.480000e-01 6.657669e-01
7.480000e-01 6.657670e-01
7.480000e-01 6.657671e-01
7.480000e-01 6.657673e-01
7.480000e-01 6.657674e-01
7.480000e-01 6.657675e-01
Notice that r= 0.743 to 0.746 only appear once. For those
values of r, convergence is complete (to the printed accuracy)
after 1000 iterations. However, starting with .747, converence
is so slow that we do not just have 1 or 2 values corresponding
to a stable attractor of period 1 or 2.
We can use awk to count the number of lines for each value of r.
Here is a little awk script pts.a from my home directory:
BEGIN {old_r=0; count=0;}
{if($1 != old_r) {print old_r,count; old_r=$1; count=1;} else ++count;}
END {print old_r,count;}
Combine the commands:
fixed_pts .7 .999 0.001 1000 128 | sort -u | awk -f pts.a |mo
Part of the output now looks like this:
7.440000e-01 1
7.450000e-01 1
7.460000e-01 1
7.470000e-01 5
7.480000e-01 106
7.490000e-01 128
7.500000e-01 128
7.510000e-01 41
7.520000e-01 2
7.530000e-01 2
7.540000e-01 2
7.550000e-01 2
7.560000e-01 2
and another part looks like this:
8.590000e-01 2
8.600000e-01 2
8.610000e-01 2
8.620000e-01 128
8.630000e-01 5
8.640000e-01 4
8.650000e-01 4
8.660000e-01 4
In this way we can find ranges of r that correspond to the
transitions between stable attractors of different periods. You
may then bracket these regions, decrease the granularity in r
and increase the number of initial iterations that are thrown
away. It is important to note that the first value of r for
which the number of unique values increases will likely be a
little below the bifurcaton point.
This process can be tedious. In the text, you will find a table
of the bifurcation points to 6 digit accuracy: Table 6.1 page 153.
Can you improve on this table?
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.
My may verify this by using map with my value of r and the one in
the book. You should find that after the expected period with my value
x we returne to 0.5, but with the value in the book, it is slightly off.
If you want full credit on the homework problem involving superstable
fixed points, you would do well to use as many as you can.