November 2, 2004 --- Class 19 --- Relaxation, Runge-Kutta, Logistic
Map and Chaos (continued)
Activities:
Relaxation Time
Some people misinterpreted the definition of relation time for the
damped harmonic oscillator. You need to consider the peak height
in every half-cycle of the oscillator. It is these peak heights
that decay exponentially. How can we easily see this? One way to
to make a semi-log plot of the trajectory. However, half the time
the coordinate is negative, so we must change the sign of the negative
values. This is easily done with awk. In a file abs_traj.a we
have:
{if($2 <0) print $1,-$2 else print $1,$2}
Then if your program output has time in column 1 and position in
column 2 and is called data:
awk -f abs_traj.a data |axis y l |plot -T x
will make a nice graph. There is a series of bumps of decreasing
amplitude and it is obvious that the peaks of the bumps fall on a
straight line. The slope of that line determines the relaxation
time. The peaks occur when the velocity changes sign. If you
identify the time and position where they occur, you can easily
solve for the relation time.
An alternative procedure to show the exponential decay is to draw
on the graph the "envelope" around the curve showing the
oscillation. If the initial amplitude is 1, the envelope is
+-exp(-t/tau). We used Mathematica to make a table of values of
the exponential.
data = Table[ {t,Exp[-t/4.]},{t,0,10,.2}]
We exported the data to a file with this command:
Export["nov2.dat", data, "Table"]
After adding a blank label to the last line, we are ready to plot
all the curves:
(awk '{$2= -$2; print}' nov2.dat; cat nov2.dat data )|axis |plot -Tx
The order of the files makes use of the fact the the two envelope
curves end with blank labes and the data file does not. Also, only
the first two columns (time and position) are plotted from the data
file. You might instead use:
(awk '{$2= -$2; print}' nov2.dat; cat nov2.dat; awk -f traj.a data ) \
|axis |plot -Tx
Fourth Order Runge-Kutta
A number of students did not use the 4th order Runge-Kutta
algorithm as described in class. They used use some strange
notation where instead of having k1x, k1v, k2x, k2v, etc. they just
used k1, k2, k3 and k4. These variables were first defined in
terms of the acceleration function. For example:
k1 = dt*accel(omega_sqr, (*y_ptr) );
k2 = dt*accel(omega_sqr,(*y_ptr+0.5*k1));
The real RK algorithm uses:
k1v = dt*accel(omega_sqr, (*y_ptr) );
k1x = dt*(*v_ptr);
k2v = dt*accel(omega_sqr,(*y_ptr+0.5*k1x));
k2x = dt*(*v_ptr + 0.5*k1v);
Notice that in the third line for k2v we use the argument:
*y_ptr+0.5*k1x. On the other hand, in the k2 line the argment is
*y_ptr+0.5*k1, where k1 is equivalent to k1v. We can't have k1x
in one place and k1v in the other. Thus, the first set of two
lines is totally bogus. Another way of seeing this is by
dimensional analysis. k1 and k1v have the dimensions of a
velocity. So, you clearly can't add 0.5*k1 to a position because
that is adding a lenght to a velocity.
I ran the bogus algorithm and showed that energy conservation did
not hold very well in comparison to the properly coded 4th order
Runge-Kutta. Woe to those who don't check energy conservation when
they are itegrating the equations of motion for a conservative
system!
Chaos and the Logistic Map
Previously, we introduced the logistic map from
consideration of an insect population model.
Upon scaling the population, the population is
replaced by x which lies in the range [0,1]. The next generation
is given by
x <- 4 r x * (1-x) .
A program entitled map 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.
We have seen 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
Students were asked to try different starting values. We then
constructed a shell script to loop over different starting values.
#!/bin/csh
set r = 0.1
foreach x (0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9)
map $r $x 100 >>! temp
end
When r is sufficiently small, no matter what the starting
population, after a long time the population approaches 0.
Next, we tried r=0.30. For any starting value other than 0 or 1,
after a number of generations the population
approaches 0.166667
With r=0.6, 0.7 and 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.
In the next class, we will explore yet larger r values and try to
understand what why this is happening.