Problem set #3 PHY256 Problems with Solutions¶
Problem 1. Baker map's strange attractor is a fractal, (precision)¶
The Baker map is a prototype for 2D chaotic dynamical systems.
On the unit inverval ($x \in [0,1]$ and $y \in [0,1]$) and for $0<c<1/2$
For c=1/2 the map is area preserving, otherwise a distribution of points covers a smaller and smaller region as the map is iterated. Below I show a code that plots a few thousand iterations of the map for c=1/3. As the map is iterated, points converge to a set of points that is a fractal and is called a strange attractor.
solution¶
The fractal rescales with factor fac = 3; power=1 If the precision is less than about 6 I find the routine converges to fixed points and you don't see the strange attractor at al. Convergence is extremely sensitive to the initial condition.
#Baker map excample using the decimal package to increase precision of computations
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# I increase precision by using the Decimal package
# https://docs.python.org/2/library/decimal.html
from decimal import *
# you can change number of digits of precision here
getcontext().prec =30
print ("digits of precision ", getcontext().prec) # number of digits of precision
d2 = Decimal(2) # using package decimal we need to give constants this way
d1 = Decimal(1)
dhalf = d1/d2
# baker map function returns x,y with coefficient c should be in 0<c<1/2
# this routine will not work if you pass it regular numbers,
# you need to pass decimals
def baker(x,y,c):
xout = d1 + c*(x-d1) # using decimal constants
yout = d1 + d2*(y-d1)
if (y<=dhalf):
xout=c*x
yout=d2*y
#print(xout,yout) # for testing
return xout,yout
# return two vectors of nit iterations of the Baker map
# xi,yi is initial conditions
# c is coefficient for Baker map
def vecbaker(xi,yi,c,nit):
xout =np.array([xi])
yout =np.array([yi])
x=xi
y=yi
for i in range(nit):
x,y = baker(x,y,c)
xout=np.append(xout,x)
yout=np.append(yout,y)
return xout,yout
c=d1/Decimal(3.0) # map parameter, with c<1/2 we can have a strange attractor
x0=dhalf # initial condition
y0=Decimal(2.0).sqrt()/Decimal(2) # an irrational number
print("initial condition:(",x0,",",y0,")")
# irrational initial condition in y required for non-periodic behavior
nit=60000 # numbers of iterations for the map
xv,yv = vecbaker(x0,y0,c,nit) # get some vectors of stuff to plot
plt.figure(figsize=(8,8))
plt.xlabel("x", fontsize=22)
plt.ylabel("y", fontsize=22)
plt.plot(xv,yv,'ro', ms=1) # plot the iterations of the map
# stuff to rescale plotting range, need to assign fac and power
#solution giving these factors
fac=3.0 # solution here! changing factor
power=2.0
win=1.0/fac**power
plt.xlim([0,win])
plt.ylim([0,win]);
fac=1.0 # original plot!
power=1.0
win=1.0/fac**power
plt.figure(figsize=(8,8))
plt.xlabel("x", fontsize=22)
plt.ylabel("y", fontsize=22)
plt.plot(xv,yv,'bo', ms=1) # plot the iterations of the map
plt.xlim([0,win])
plt.ylim([0,win]);
getcontext().prec =6
print ("digits of precision ", getcontext().prec) # number of digits of precision
nit=60000 # numbers of iterations for the map
xv,yv = vecbaker(x0,y0,c,nit) # get some vectors of stuff to plot
plt.figure(figsize=(8,8))
plt.xlabel("x", fontsize=22)
plt.ylabel("y", fontsize=22)
plt.plot(xv,yv,'ro', ms=1); # plot the iterations of the map
# stuff to rescale plotting range, need to assign fac and power
#solution giving these factors
fac=3.0
power=2.0
win=1.0/fac**power
plt.xlim([0,win])
plt.ylim([0,win]);
# looks bad, the plot does not look the same as above
# now change initial condition
getcontext().prec =30
print ("digits of precision ", getcontext().prec) # number of digits of precision
nit=10000 # numbers of iterations for the map
x0=Decimal(1)/Decimal(2) # initial condition
y0=Decimal(1)/Decimal(8) # a rational number
#this gives a periodic orbit
xv,yv = vecbaker(x0,y0,c,nit) # get some vectors of stuff to plot
plt.figure(figsize=(8,8))
plt.xlabel("x", fontsize=22)
plt.ylabel("y", fontsize=22)
plt.plot(xv,yv,'ro', ms=1); # plot the iterations of the map
# stuff to rescale plotting range
fac=1.0
power=1.0
win=1.0/fac**power
plt.xlim([0,win])
plt.ylim([0,win]);
Problem 2. Baker map's Lyapunov exponent¶
A symptom of chaotic behavior is exponential divergence of nearby trajectories. The maximum Lyapunov exponent describes how fast trajectories diverge. In a dynamical system this gives a timescale for how fast information of the initial conditions is lost.
We can record trajectories with two different initial conditions, where the two different initial positions are very close together, only separated by a small distance of $\epsilon$.
We can chose $\epsilon$ to be near the precision limit. For example, $\epsilon =$ 1E-41 if our numbers of digits of precision is set to 42.
Using the Baker map, compute vectors \begin{eqnarray} (x^i, y^i) &=& f^i(x_0,y_0)\\ (x_*^i, y_*^i) &=& f^i(x_0,y_0 + \epsilon) \end{eqnarray} For $i \in [0,100]$ (you can increase the range later on).
Compute the distance between each pair of iterates \begin{equation} d^i = \sqrt{(x^i - x_*^i)^2 + (y^i - y_*^i)^2} \end{equation}
If the difference depends exponentially on the number of steps then $$ d^i = d_0 e^{\lambda i}$$ Ideally the maximum Lyapunov exponent $$ \lambda = \lim_{i \to \infty} \frac{1}{i} \log \left| \frac{d_i}{d_0} \right|$$ It's a maximum exponent because the system has two dimensions and trajectories are only stretching in one direction.
Make a plot of $\log d^i$ vs $i$. While the idea is to take the limit to infinity, in practice if $i$ is too large, the distance levels off because the difference between trajectories reaches the boundary of the unit square. So you need to compute the exponent from the region in the plot where the slope is linear.
Use your plot to estimate the Lypunov exponent $\lambda$.
The Lyapunov exponent for the Baker map is $\ln 2$. How close is your measured exponent to that value?
Notes: Assuming you are working with the code from the previous problem set, everytime you do arithmetic with a constant it must be a decimal. https://docs.python.org/2/library/decimal.html np.log() doesn't work on decimal arrays but np.log10() does. If you increase the number of iterations you should see your plot level off where the log is near zero corresponding to the difference being as large as the unit interval (can't get any larger than that!)
#Solutions
#Baker map excample using the decimal package to increase precision of computations
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# I increase precision by using the Decimal package
# https://docs.python.org/2/library/decimal.html
from decimal import *
print (getcontext().prec)
# you can change this
getcontext().prec = 42 # number of digits of precision
# needs to be high!
d2 = Decimal(2.0) # using decimal we need to give constants this way
d1 = Decimal(1.0)
dhalf = d1/d2
# baker map function returns x,y
# this routine will not work if you pass it combinations of regular numbers and decimals,
# you need to pass decimals
def baker(x,y,c):
xout = d1 + c*(x-d1) # using our decimal constants
yout = d1 + d2*(y-d1)
if (y<=dhalf):
xout=c*x
yout=d2*y
return xout,yout
# return two vectors of nit iterations of the Baker map
# xi,yi is initial conditions
# c is coefficient for Baker map
def vecbaker(xi,yi,c,nit):
xout =[xi]
yout =[yi]
x=xi
y=yi
for i in range(0,nit):
x,y = baker(x,y,c)
xout=np.append(xout,x)
yout=np.append(yout,y)
return xout,yout # vectors
c=d1/Decimal(3.0) # map parameter, with c<1/2 we should have a strange attractor
x0=dhalf # initial condition
y0=Decimal(2.0).sqrt()/Decimal(2) # irrational initial condition required for non-periodic behavior
nit=200 # numbers of iterations for the map
eps = Decimal(1e-41) # a small number, the initial distance between trajectories
leps = eps.ln()
# first trajectory
xv,yv = vecbaker(x0,y0,c,nit) # get some vectors of stuff to plot
# second trajectory
xv2,yv2 = vecbaker(x0,y0+eps,c,nit) # comparison with slighly different initial conditions
diffx = xv - xv2
diffy = yv - yv2
# compute the distances between trajectories
diffr = np.sqrt(diffx*diffx + diffy*diffy)
ldiff = diffr*Decimal(0.0);
# computing the log of the difference array
for i in range(nit): # doing this by hand because np.log bombs on Decimals
ldiff[i] = diffr[i].ln()
plt.figure(figsize=(8,8))
plt.xlabel("iterate i", fontsize=22)
plt.ylabel("$\ln(d^i$)", fontsize=22)
plt.plot(ldiff,'ro', ms=5, label = 'distance') # plot the log of distance between trajectories
xlin = np.arange(1,150,1)
plt.plot(xlin*np.log(2)-96,'b-', label='predicted') # blue line shows that log2 is a pretty good match to slope
iss = 101 # a point to use to measure the slope
slope = (ldiff[iss] - leps)/Decimal(iss)
#leps is the log of the initial distance apart
print("lambda=",slope) #Lyapunov exponent computed numerically
dln2 = d2.ln() # predicted Lyapunov exponent
diff_slope = slope - dln2 # difference between predicted and computed
print("lambda - ln2=",diff_slope) #how far are we from what is expected
plt.legend(loc = 'upper left')
# result is that my measurement gives a number very very close
# to the predicted ln2
