#Solution script for Math Methods Problem Set 2 #Problem 1: Analyzing an interated map dynamical system import math from ClimateUtilities import * import map1D #Define the map function def f(x,a): return a*x*math.cos(x) #Use Newton's method to find the stability #boundaries for the fixed points. def fplus(a,j): return (math.acos(1./a) + 2.*j*math.pi)*math.sqrt(a*a-1.) - 2. def fminus(a,j): return (-math.acos(1./a) + 2.*j*math.pi)*math.sqrt(a*a-1.) - 2. rootPlus = newtSolve(fplus) rootMinus = newtSolve(fminus) #Compute and print out stability boundaries print "Stability boundaries for fixed points:" rootPlus.params = 0. print "Smallest fixed point goes unstable at a = ",rootPlus(1.001) print "The next few positive fixed points go unstable at a =:" for j in range(1,6): rootPlus.params = j rootMinus.params = j print " ",rootMinus(1.001) print " ",rootPlus(1.001) print "-----------------------------------------" #Find the position of the first few positive maxima/minima #using Newton's method def g(x,dummy): return x*math.tan(x) - 1. r = newtSolve(g) for j in range(5): if j ==0: guess = .5 else: guess = j*math.pi #Use approximate value as guess xm = r(guess) print "Maximum at x = %f, value = %f *a"%(xm,xm*math.cos(xm)) #Now we try out some ideas with mapExplorer m = map1D.mapExplorer(f) #Look at the map for a = 1.5, and a typical orbit. #In this regime, there is a stable positive fixed point. m.param = 1.5 m.plotComposition(1,[0.,math.pi/2]) m.plotOrbit(.01,20) #Note that the solution is attracted to the stable fixed point. # #We suspect that when the fixed point goes unstable, it bifurcates #to a period-2 orbit. To verify this, we look at the composition #plot just a little over the stability boundary. m.param = 2.2 m.plotComposition(2,[0.,math.pi/2.]) m.plotOrbit(.01,20) #Look for the stability boundary of the period-2 orbit, #and check for the emergence of a stable period-4 orbit #You can do this by just plotting the orbits for various #values of a, and looking for when the orbit is attracted #to a period-4 rather than a period-2 orbit. Here, instead #we use the findPeriodicOrbits method, and vary a until #the period-2 orbit goes unstable. Remember that the #findPeriodicOrbits method reports the exponential growth rate, #not the slope. Thus, a negative value is stable, and a positive #value is unstable. There is only one stable orbit, so we can #use Python's min(...) function, which finds the minimum value #in a list, and then look for when this goes positive (unstable) print "Looking for stability boundary of period-2 orbit" for a in [2.2+ .002*i for i in range(100)]: m.param = a orbits,stability=m.findPeriodicOrbits(2,[0.,math.pi/2]) print a,min(stability) #From looking at the output of this, it looks like the stability #boundary is near a = 2.393. Let's look at some composition maps #and orbits above the stability boundary, to check m.param = 2.44 # #Look at the composition plot m.plotComposition(4,[0.,math.pi/2]) #Find the unstable period-2 orbit to use as an initial condition orbits,stability = m.findPeriodicOrbits(2,[0.,math.pi/2]) xstart = orbits[1] #To find which one to take, we just have to look #at the list, and ignore the ones that we know #are actually fixed points, not periodic orbits m.plotOrbit(xstart,80) #Note how the period-2 orbit breaks down into a period-4 orbit #Plot the orbit diagram to look for the onset of a dense #set of periodic orbits. We vary the parameter a in #the region between the stability boundary of the period-2 #orbit and the upper limit of the value that guarantees #that the orbit stays in the interval. import orbitDiagram orbitDiagram.orbitDiagram(f,2.2,2.79,[0.,math.pi/2]) #Now let's look at some chaotic orbits. m.param = 2.7 orbit1 = m.orbit(.5,100) orbit2 = m.orbit(.50001,100) map1D.ezplot(orbit1,orbit2) #Now compute the log of the difference, and estimate the #Lyapunov exponent diff = [abs(orbit1[i]-orbit2[i]) for i in range(len(orbit1))] logdiff = [math.log(d) for d in diff] map1D.ezplot(logdiff[0:30]) #Just plot first 30 #Looks like the exponentially growing phase is between iterate 1 and #20. Take the slope of this segment to estimate the lyapunov exponent. slope = (logdiff[20]-logdiff[0])/20 print "The Lyapunov exponent is approximately ",slope #-------Problem 2: Computing a histogram--------------- def count(interval, data): sum = 0 for x in data: if (x < interval[1])&(x >= interval[0]): sum += 1 #shorthand for sum = sum+1 return sum #Arguments: # bins is a list of the boundary points of the bins # in ascending order # # data is a list of the data points whose histogram we # wish to compute. def histo(bins,data): counts = [] for i in range(len(bins)-1): counts.append(count([ bins[i],bins[i+1] ],data)) return counts #Now test it out: #First make some data import math n = 1000 data = [math.sin(i*math.pi/n) for i in range(n+1)] #Now define the bins nbins = 100 bins = [i/float(nbins) for i in range(nbins+1)] #Note that the float(...) function converts the integer to a float #Compute the histogram myHisto = histo(bins,data) #Now plot it out # First we compute the midpoints of the bins, for plotting purposes midpoints = [.5*(bins[i]+bins[i+1]) for i in range(len(bins)-1)] # # Now stuff it in a curve and plot it from ClimateUtilities import * c = Curve() c.addCurve(midpoints) c.addCurve(myHisto) w = plot(c) #The curve looks a little bumpy with n = 1000 and 100 bins. This #is because of the relatively small number of points per bin. #Experiment with changing n and nbin. # #Note the apparent singularity in the histogram near x = 1. Can #you explain this in terms of the shape of sin(x)? #---------------------Extra for the PowerPythonista:---------------------- # #Note that the algorithm used above is the simplest way to #program a histogram computation, but it is far from the most #efficient. The number of operations is proportional to nbin*n. #If the bins are equally spaced, the histogram can actually be #calculated with only order(n) operations, since we can immediately #find which bin a data point x belongs to using index = int((x-x1)/dx), #where x1 is the coordinate of the first bin boundary and dx is the #bin width. The int(...) function converts the float to its integer part. # #Even if the bins aren't equally spaced, the search for the right bin #can be speeded up considerably using the fact that the bin boundaries #are sorted in ascending order. If the bins are not too nonuniform, #a binary search can reduce the count to log2(nbin) for each data point. #If one knows something more about the nature of the bins, the count #can be reduced even more. # #Either of these cases requires changing the structure of the program. #Instead of looping over datapoints in count(...), we need to #write a routine called findIndex(bins,x) which returns the bin #index corresponding to datapoint x. Then, we loop over all datapoints, #incrementing the appropriate bin. The following code illustrates #how to do this. We don't write a findIndex(...) function ourselves, #because the Numeric array package includes a function called searchsorted(a,x) #which already does this for us, and does it efficiently using #compiled code. Note that searchsorted assumes that the data in the array #is sorted in ascending order. It returns index of the first element #in the array a which exceeds x, and returns len(a) (the #highest index of array plus 1) if x is larger than the last element. # #If the bins are equally spaced, the use of Numeric.searchsorted(...) #can be replaced by the arithmetic formula for the index # #The following code assumes a little familiarity with the Numeric array #module, mainly because we want to use the efficient Numeric.searchsorted(...) #function. import Numeric #The routine will work even if the input argument bins is not #converted to a Numeric array, but it will be very inefficient, #because Numeric will do the conversion on the fly for each call #to searchsorted. To get the benefit of the speedup, it is necessary #to convert bins to an array. Conversion of data is optional. To make #life easier for the user, we do the conversion in the function. def betterHisto(bins,data): #Convert bins to a Numeric array, for more efficiency bins = Numeric.array(bins) #Note: this won't affect bins in the calling program! #First make an integer array to hold the counts. Initially, #all the elements are zero. counts = Numeric.zeros(len(bins)-1,Numeric.Int) # #Now loop over the data and do the counts. Test the index # to make sure the data is not out of the range of the bins. for x in data: i = Numeric.searchsorted(bins,x) if (i>0)&(i