#Program le-spectrum.py calculates the spectrum of LEs for the Henon map # Adapted from # http://sprott.physics.wisc.edu/chaos/wolfmap.bas def derivs(x, xnew, n): return tinkerbell(x, xnew, n) def Henon(x, xnew, n): a=1.4 b=0.3 # Nonlinear Henon map equations: xnew[1] = 1-a*x[1]*x[1]+b*x[2] xnew[2] = x[1] # Linearized Henon map equations: xnew[3] = -2*a*x[1]*x[3]+b*x[5] xnew[4] = -2*a*x[1]*x[4]+b*x[6] xnew[5] = x[3] xnew[6] = x[4] return [x, xnew] # Seems to lead to an error after a few thousand iterations def DelayedLogisticMap(x, xnew, n): a=2.27 # Nonlinear Delayed Logistic Map equations: xnew[1] = a*x[1] - a*x[1]*x[2] xnew[2] = x[1] # Linearized Delayed Logistic Map equations: xnew[3] = a*x[3]-a*x[3]*x[2] - a*x[1]*x[5] xnew[4] = a*x[4]-a*x[4]*x[2] - a*x[1]*x[6] xnew[5] = x[3] xnew[6] = x[4] return [x, xnew] def burger(x, xnew, n): a = .75 b = 1.75 # Nonlinear xnew[1] = a*x[1] - x[2]**2 xnew[2] = b*x[2] + x[1]*x[2] # Linearized xnew[3] = a*x[3] - 2 * x[2]*x[5] xnew[4] = a*x[4] - 2 * x[2]*x[6] xnew[5] = x[2]*x[3] + b*x[5] + x[1]*x[5] xnew[6] = x[2]*x[4] + b*x[6] + x[1]*x[6] return [x, xnew] def tinkerbell(x, xnew, n): a = 0.9 b = -0.6 c = 2.0 d = 0.5 # Nonlinear xnew[1] = x[1]**2 - x[2]**2 + a*x[1] + b*x[2] # x xnew[2] = 2*x[1]*x[2] + c * x[1] + d * x[2] # y # Linearized # To linearize these equations # Find the jacobian of the system # Then multiply the Jn by the 1 x 2 vector [delta x; delta y] to get # a 1 x 2 vector [delta x+1; delta y+1] # xnew[3] = delta x(n + 1), xnew[5] = delta y(n+1) xnew[3] = 2 * x[1] * x[3] + a * x[3] - 2 * x[2] * x[5] + b * x[5] # delta x xnew[4] = 2 * x[1] * x[4] + a * x[4] - 2 * x[2] * x[6] + b * x[6] # delta x xnew[5] = 2 * x[2] * x[3] + c * x[3] + 2 * x[1] * x[5] + d * x[5] # delta y xnew[6] = 2 * x[2] * x[4] + c * x[4] + 2 * x[1] * x[6] + d * x[6] # delta y return [x, xnew] n=2 # number of variables in nonlinear map nn=n*(n+1) # total number of variables (nonlinear + linear) m = 0 x = [] xnew = [] v = [] ltot = [] znorm = [] gsc = [] i = 0 while i <= nn: x.append(0) xnew.append(0) v.append(0) i = i + 1 i = 0 while i <= n: ltot.append(0) znorm.append(0) gsc.append(0) i = i + 1 irate=10 # integration steps per reorthonormalization io= 1000 # number of iterations between printouts cmax = 100000 # number of iterations to perform # initial conditions for nonlinear maps # must be within the basin of attraction i = 1 while i <= n: v[i] = 0.001 i = i + 1 i = n+1 while i <= nn: # initial conditions for linearized maps v[i]=0 # Don't mess with these; they are problem independent! i = i + 1 i = 1 while i <= n: v[(n+1)*i]=1 ltot[i]=0 i = i + 1 print v t=0 w = 0 while (w < cmax): j = 1 while j <= irate: i = 1 while i <= nn: x[i]=v[i] i = i + 1 [x, xnew] = derivs(x, xnew, n) i = 1 while i <= nn: v[i]=xnew[i] i = i + 1 t = t + 1 j = j + 1 #construct new orthonormal basis by gram-schmidt: znorm[1]=0 #normalize first vector j = 1 while j <= n: znorm[1]=znorm[1]+v[n*j+1]**2 j = j + 1 znorm[1] = math.sqrt(znorm[1]) j = 1 while j <= n: v[n*j+1]=v[n*j+1]/znorm[1] j = j + 1 #generate new orthonormal set: j = 2 while j <= n: k = 1 while k <= j-1: gsc[k]=0 l = 1 while l <= n: gsc[k]=gsc[k]+v[n*l+j]*v[n*l+k] l = l + 1 k = k + 1 k = 1 while k <= n: #construct a new vector l = 1 while l <= j-1: v[n*k+j]=v[n*k+j]-gsc[l]*v[n*k+l] l = l + 1 k = k + 1 znorm[j]=0 #calculate the vector's norm k = 1 while k <= n: #construct a new vector znorm[j]=znorm[j]+v[n*k+j]**2 k = k + 1 znorm[j]=math.sqrt(znorm[j]) k = 1 while k <= n: #normalize the new vector v[n*k+j] = v[n*k+j] / znorm[j] k = k + 1 j = j + 1 k = 1 while k <= n: #update running vector magnitudes if znorm[k] > 0: ltot[k] = ltot[k] + math.log(znorm[k]) k = k + 1 m = m + 1 if m % io == 0: #normalize exponent and print every io iterations print "\tLEs ="; lsum=0 kmax=0 k = 1 while k <= n: le = ltot[k] / t lsum = lsum + le if lsum > 0: lsum0 = lsum kmax = k print le k = k + 1 if ltot[1]>0: dky = kmax - t * lsum0 / float(ltot[kmax+1]) else: dky = 0 print " Kaplan-Yorke dimension =", dky w = w + 1