############################################################ # # In this template, Henon map is implemented in a class because # we need to change the parameters for future homeworks and # it is much more convenient to use class here than global # parameters. # If you are not familiar with class object in python, # have a look at the tutorial here: # https://docs.python.org/2/tutorial/classes.html # # please set case = 1 to validate your implementation first. # Then go to case2, case3 and case4. # ############################################################ import numpy as np import matplotlib.pyplot as plt from numpy.random import rand from scipy.interpolate import splrep, splev from scipy.optimize import fsolve class Henon: """ Class Henon contains functions for both forward and backward Henon map iteration. """ def __init__(self, a=6, b=-1): """ initialization function which will be called every time you create an object instance. In this case, it initializes parameter a and b, which are the parameters of Henon map. """ self.a = a self.b = b def oneIter(self, stateVec): """ forward iterate for one step. stateVec: the current state. dimension : [1 x 2] numpy.array return: the next state. dimension : [1 x 2] """ x = stateVec[0]; y = stateVec[1]; stateNext = return stateNext def multiIter(self, stateVec, NumOfIter): """ forward iterate for multiple steps. stateVec: the current state. dimension : [1 x 2] numpy.array NumOfIter: number of iterations return: the current state and the furture 'NumOfIter' states. dimension [NumOfIter+1 x 2] Hint : numpy.vstack() """ state = return state def Jacob(self, stateVec): """ The Jacobian for forward map at state point 'stateVec'. stateVec: the current state. dimension : [1 x 2] numpy.array """ x = stateVec[0]; y = stateVec[1]; jacobian = return jacobian def oneBackIter(self, stateVec): """ backward iterate for one step. stateVec: the current state. dimension : [1 x 2] numpy.array return: the previous state. dimension : [1 x 2] """ x = stateVec[0]; y = stateVec[1]; statePrev = return statePrev def multiBackIter(self, stateVec, NumOfIter): """ backward iterate for multiple steps. stateVec: the current state. dimension : [1 x 2] numpy.array NumOfIter: number of backward iterations return: the current state and the pervious 'NumOfIter' states. dimension [NumOfIter+1 x 2] """ state = stateVec tmpState = stateVec for i in range(NumOfIter): tmp = self.oneBackIter(tmpState) tmpState = tmp state = np.vstack((state, tmpState)) return state if __name__ == '__main__': case = 4 if case == 1: """ Validate your implimentation of Henon map. Note here we use a=1.4, b=0.3 in this valication case. For other cases in this homework, we use a=6, b=-1. Actually, these two sets of parameters are both important since we will come back to this model when discussing invariant measure. """ henon = Henon(1.4, 0.3) # creake a Henon instance states = henon.multiIter(rand(2), 1000) # forward iterations states_bac = henon.multiBackIter(states[-1,:], 10) # backward iterations # check your implementation of forward map fig = plt.figure(figsize=(6,6)) ax = fig.add_subplot(111) ax.scatter(states[:,0], states[:, 1], edgecolor='none') ax.set_xlabel('x') ax.set_ylabel('y') ax.set_title('(a)') plt.show() # check the correctness of backward map. The first 10 states_bac should # be the last 10 states in reverse order. # Note, backward map is very unstable for a=1.4, b=0.3, # so we only iterate backward for 10 times. print states[-10:, :] print '======' print states_bac[:10,:] print '======' # check the Jacobian matrix at (0.1, 0.2). # for a = 1.4, b = 0.3, the output should be # [[-0.28 0.3 ] # [ 1. 0. ]] print henon.Jacob(np.array([0.1, 0.2])) if case == 2: """ Try to obtain the stable and unstable manifold for equilibrium '0' in Henon map with parameter a=6, b=-1. Plotting unstable/stable manifold is a difficult task in general. We need to moniter a lot of variables, like the distances between points along the manifold, the angle formed by adjacent 3 points on the manifold, etc. However, for Henon map, a simple algorithm is enough for demonstration purpose. The algorithm works as follows: Unstable manifold: start from a point close to equilibrium '0', in the direction of the expanding eigenvector: ' eq0 + r0 * V_e '. Here 'eq0' is the equilibrium, 'r0' is a small number, 'V_e' is the expanding eigen direction of equilibrium '0'. The image of this point after one forward iteration should be very close to ' eq0 + r0 * \Lambda_e * V_e', where '\Lambda_e' is the expanding multiplier. We can confidently think that these two points are sitting on the unstable manifold of 'eq0' since 'r0' is very small. Now, we interpolate linearly between these two points and get total 'N' points. If we iterate these 'N' points forward for 'NumOfIter' times, we get the unstable manifold within some length. The crutial part of this method is that when these 'N' state points are being iterated, the spacings between them get larger and larger, so we need to use a relative large value 'N' to ensure that these state points are not too far away from each other. The following formula is used to determine 'N' ( please figure out its meaning and convince yourself ): ( (\Lamba_e - 1) * r0 / N ) * (\Lambda_e)^NumOfIter = tol (*) Here, 'tol' is the tolerance distance between adjacent points in the unstable manifold. Stable manifold: If we reverse the dynamics, then the stable direction becomes the ustable direction, so we can use the same method as above to obtain stable manifold. """ henon = Henon() # use the default parameters: a=6, b=-1 # get the two equilbria of this map. equilibrium '0' should have smaller x coordinate. eq0 = np.array([, ]) # equilibrium '0' eq1 = np.array([, ]) # equilibrium '1' # get the expanding multiplier and eigenvectors at equilibrium '0' Lamba_e = # expanding multiplier Ev = # expanding eigenvector NumOfIter = 5 # number of iterations used to get stable/unstable manifold tol = 0.1 # tolerance distance between adjacent points in the manifold r0 = 0.0001 # small length N = # implement the formula (*). Note 'N' should be an integer. delta_r = (Lamba_e-1)*r0 / N # initial spacing between points in the manifold # generate the unstable manifold. Note we do not use Henon.multiIter() here # since we want to keep the ordering of the points along the manifold. uManifold = eq0 states = np.zeros([N,2]) for i in range(N): # get the initial N points states[i,:] = eq0 + (r0 + delta_r*i)*Ev uManifold = np.vstack((uManifold, states)) for i in range(NumOfIter): for j in range(N): # update these N points along the manifold states[j,:] = henon.oneIter(states[j,:]); uManifold = np.vstack((uManifold, states)) # Please fill out this part to generate the stable manifold. # Check whether the stable manifold is symmetric with unstable manifold with # diagonal line y = x sManifold = # get the spline interpolation of stable manifold. Note: unstable manifold are # double-valued, so we only interploate stable manifold, and this is # enough since unstable manifold and stable manifold is symmetric with line y = x. tck = splrep(sManifold[:,0], sManifold[:,1], s=0) # use scipy.optimize.fsolve() to obtain intersection points B, C, D # hint: use scipy.interpolate.splev() and the fact that stable and unstable # are symmetric with y = x C = D = B = # save the variables needed for case3 # if you are using ipython enviroment, you could just keep the varibles in the # current session. np.savez('case2', B=B, C=C, D=D, eq0=eq0, eq1=eq1, sManifold=sManifold, uManifold=uManifold, tck=tck) # plot the unstable, stable manifold, points B, C, D, equilibria '0' and '1'. fig = plt.figure(figsize=(6,6)) ax = fig.add_subplot(111) ax.plot(uManifold[:,0], uManifold[:, 1], 'r-', lw=2, label=r'$W_u$') ax.plot(sManifold[:,0], sManifold[:, 1], 'c-', lw=2, label=r'$W_s$') ax.scatter(eq0[0],eq0[1]) ax.scatter(eq1[0],eq1[1]) ax.text(C[0], C[1], 'C') ax.text(D[0], D[1], 'D') ax.text(B[0], B[1], 'B') ax.text(eq0[0], eq0[1], '0') ax.text(eq1[0], eq1[1], '1') plt.legend() plt.show() if case == 3: """ Try to establish the first level partition of the non-wandering set in Henon map. You are going to iterate region 0BCD forward and backward for one step. """ henon = Henon() # use the default parameters: a=6, b=-1 # load needed variables from case2. If you have kept the variables # in the working space, then just comment out the following few lines. case2 = np.load('case2.npz') B = case2['B']; C = case2['C']; D = case2['D']; eq0 = case2['eq0']; eq1 = case2['eq1']; tck = case2['tck']; uManifold = case2['uManifold']; sManifold = case2['sManifold']; # We first make a sampling of region 0BCD. # It works like this: # we are sure that box [-0.8, 0.8] x [-0.8, 0.8] # can cover region OBCD, but not every point inside this box is in 0BCD, # so for point (x, y), how to determine whether it is in 0BCD ? # The criteria is # y < f(x) and x < f(y) # Here f() is the interpolation function of stable manifold. # It is easy to see that 'y < f(x)' enforces point (x, y) below the # stable manifold, but the fact that 'x < f(y)' enforces point (x, y) to be # at the left side of the unstable manifold is a little tricky. The answer # is that stable manifold and unstable manifold are symmetric with y = x # Anyway, this part is implemented for you. M = np.array([]).reshape(0,2) # region 0BCD x = np.arange(-0.8, 0.8, 0.01) y = splev(x, tck) for i in range(np.size(x)): for j in range(np.size(x)): if x[i] < y[j] and x[j] < y[i]: state = np.array([x[i], x[j]]) M = np.vstack( (M, state) ) # please plot out region M to convince yourself that you get region 0BCD # Now iterate forward and backward the points in region 0BCD for one step Mf1 = # forward iteration points Mb1 = # backward iteration points # plot out Mf1 and Mb1 fig = plt.figure(figsize=(6,6)) ax = fig.add_subplot(111) ax.plot(Mb1[:,0], Mb1[:,1], 'g.') ax.plot(Mf1[:,0], Mf1[:,1], 'm.') ax.plot(uManifold[:,0], uManifold[:, 1], 'r') ax.plot(sManifold[:,0], sManifold[:, 1], 'c') plt.show() # In order to see the pre-images of the boarders of Mf1 and Mb1, please # try to plot the images and per-images of 4 edges of region 0BCD. # hint: use the interpolation function of stable manifold if case == 4: """ We go further into the partition of state space in this case. In case3 you have figure out what the pre-images of the border of first forward and backward iteration, so we do not need to sample the region again, iteration of the boarder is enough. In this case we iterate forward and backward for two steps """ henon = Henon() # use the default parameters: a=6, b=-1 # load needed variables from case2 case2 = np.load('case2.npz') B = case2['B']; C = case2['C']; D = case2['D']; eq0 = case2['eq0']; eq1 = case2['eq1']; tck = case2['tck']; uManifold = case2['uManifold']; sManifold = case2['sManifold']; # initialize the first/second forward/backward iteration of the boarder Mf1 = np.array([]).reshape(0,2) # the first forward iteration of the boarder you got in case 3 Mf2 = np.array([]).reshape(0,2) # ... second forward .... Mb1 = np.array([]).reshape(0,2) # ....first backward .... Mb2 = np.array([]).reshape(0,2) # ....second backward .... # implement your code here to get Mf1, Mf2, Mb1, Mb2 # hint: use the interpolation function of stable manifold # plot out your result. fig = plt.figure(figsize=(6,6)) ax = fig.add_subplot(111) ax.plot(uManifold[:,0], uManifold[:, 1], 'r') ax.plot(sManifold[:,0], sManifold[:, 1], 'c') ax.plot(Mf1[:,0], Mf1[:,1], 'm') ax.plot(Mf2[:,0], Mf2[:,1], 'g') ax.plot(Mb1[:,0], Mb1[:,1], 'b') ax.plot(Mb2[:,0], Mb2[:,1], 'y') ax.set_title('(d)') plt.show() # find a point in the top left region (the region which is closest to point D) # as the initial condition to find a periodic period with period 4 # hint: use fsolve() guess = np.array([-0.4, 0.5]) x = # the initial condition you get from this guess print henon.multiIter(x, 4) # check whether it is periodic # if you like, you can figure out the symbolic representation # of this periodic orbit.