#ToDo: On these example scripts, we need to label the graphs, #and also find a way to allow the student to step through a long #script in an interactive way, so he/she can see what is going on. # # #Solution script for Math Methods, Problem Set 3 from ClimateUtilities import * import math #Problem 1 def f(x,junk): return math.sin(x) - 10.*math.exp(-x) #Instantiate a Newton solver root = newtSolve(f) #The only real challenge in this problem is to find #a good set of initial guesses. In this case, since #the exponential is positive and monotonically decreasing but the #sin oscillates and goes negative, we can infer that there is a #chance for a pair of solutions between each two zeros of the #sin function enclosing a positive loop of sin(x). (Sketch a #graph of sin(x) and 10exp(-x) if this is unclear -- the script #does that for you to help out.) Therefore, if we use the zeros #of sin(x) as initial guesses, and try each one in succession until #we have five different roots, we have a pretty good chance #of finding the first five. To help increase our chances that we #didn't miss one, we actually try some intermediate values as #well. #Plot the functions to get an idea of what is going on x = [i*20*math.pi/100. for i in range(100)] sinx = [math.sin(xx) for xx in x] expx = [10.*math.exp(-xx) for xx in x] c = Curve() c.addCurve(x) c.addCurve(sinx,'sin','sin(x)') #Gives the column a name and a label c.addCurve(expx,'expx','10*exp(-x)') w1 = plot(c) #From looking at the curve, you should see that in fact the farther #out zeros are quite close the the zeros of sin, since the exponential #decays so rapidly. #Since Newton's method blows up if the derivative of the #objective function is small near the initial guess, we avoid #those points (hence the pi/7). print 'Problem 1 solution' for guess in [i*math.pi/7 for i in range(100)]: print root(guess) #From this list we can see some of the peculiar features of Newton's #method. It can't be counted on to converge to the closest root. Still, #by examining the list, we can pick out the values we want. #Problem 2: Newton's method as a dynamical system import map1D def g(x,junk): return .5*x - .5/x #Instantiate a map explorer m = map1D.mapExplorer(g) #Let's find the pre-images of zero, so we know what points #will blow up after a given number N of iterates. We do #this by solving x=g(x) for g, yielding x = g + sqrt(g*g+1) and #x = g - sqrt(g*g+1). Note this is a multiply valued function: #more than one point will map to zero after N iterates. Therefore, #we compile an ever-lengthening list of pre-images of zero, by iterating #over the previous list print "Pre-images of zero at n iterates" preImages = [0] for n in range (4): #Note how we use list comprehension and list concatenation #to avoid writing a loop. Much prettier this way preImages = [gg+math.sqrt(gg*gg+1) for gg in preImages]+ \ [gg-math.sqrt(gg*gg+1) for gg in preImages] preImages.sort() #Sort the list (in place) print n+1,preImages #Look at the plots for the function and the first few #compositions, to get an idea of whether there are any #periodic orbits #Take a look at the iteration function itself. As expected, #no fixed points. Note the blow-up at x=0. Because we know #g(-x) = -g(x) we don't have to look at negative values separately w3 = m.plotComposition(1,[.1,.9]) w4 = m.plotComposition(1,[1.1,10]) #Period two orbits in the unit interval w5 = m.plotComposition(2,[.1,.9]) #Looks like there's an unstable period 2 orbit near .6. Let's #find it more precisely print "Period 2 orbit characteristics:",m.findPeriodicOrbits(2,[.5,.6]) #Any periodic orbits for larger x? We are avoiding the singularity w6 = m.plotComposition(2,[1.1,10]) #Looks like there aren't. Now let's look at what happens when we #let our period-2 orbit go unstable. w6 = m.plotOrbit(0.57727727727727729,50) # #For a while, the orbit alternates between .5769 and -.5769. #After 10 iterates, the instability has grown to the point #where the orbit no longer looks periodic. The general description #of the orbit is that it oscillates for a while at moderate values, #then every once in a while comes close to zero, after which it #maps to very large values. Then, the orbit relaxes exponentially #(like 1/2**n) toward small values, and the process starts #all over again. Let's look at a really long orbit, to see how big #the values get in the long term w7 = m.plotOrbit(.5,1000) #Interesting question: What is the maximum value encountered after #n interations? What is the form of the probability distribution #of the large x tail? (i.e. the probability of encountering a #value larger than X). #We've now seen most of the interesting behavior, but let's #finish the rest of the problem. First look for period-3 #orbits. To avoid problems with singular points, we need to #look separately on the sub-intervals where the composition 3 #map is non singular (guided by the pre-images of zero calculated #above) w8 = m.plotComposition(3,[.1,.35]) w9 = m.plotComposition(3,[.45,.9]) w10 = m.plotComposition(3,[1.1,2.2]) w11 = m.plotComposition(3,[2.5,10]) #Looking at the graphs, we find an unstable period-3 orbit starting from each of #the first three subintervals. Let's use findPeriodicOrbits to get accurate #values orbits3A = m.findPeriodicOrbits(3,[.1,.35]) orbits3B = m.findPeriodicOrbits(3,[.45,.9]) orbits3C = m.findPeriodicOrbits(3,[1.1,2.2]) print 'Period 3 orbit characteristics:' print orbits3A,orbits3B,orbits3C #Note that these are all unstable (the growth rates are all positive) # #Is this a single period-3 orbit, or #are there two or more distinct ones? You can look at some #iterates to find out. This detail is left to the student. #Next, lets look for a period-4 orbits. We won't exhaust ourselves #by doing an exhaustive search. Just looking in the first sub-interval, #we find w12 = m.plotComposition(4,[.01,.15]) #we find an unstable periodic orbit here, too. print 'Characteristics of one of the Period-4 Orbits:' print m.findPeriodicOrbits(4,[.01,.15]) #There are other unstable periodic orbits in the other subintervals #Now, finally, we check for sensitive dependence on initial conditions, #to look for one of the symptoms of chaos. orbit1 = m.orbit(.5,30) orbit2 = m.orbit(.500001,30) logdiff = [math.log(abs(orbit1[i]-orbit2[i])) for i in range(len(orbit1))] w13 = map1D.ezplot(logdiff) #Seems to have sensitive dependence on initial conditions #We can also use the lyap method to compute the lyapunov exponent #over n iterates. Recall that the lyapunov exponent is the exponential #rate of separation of nearby trajectories. Try out computing the #lyapunov exponent for different n and for different initial conditions. #Do you always get approximately the same answer? Try looking at #different initial values for small n and for large n. What do you notice? #You could try making some graphs of the lyapunov exponent vs. initial #x, to explore this question. print 'The lyapunov exponent is approximately',m.lyap(.5,1000) #Note that the lyapunov exponent is quite close to the growth #rate of instabilities to the periodic orbits we discovered earlier. #Finally, do a quick and dirty check of 'ergodicity' , focusing #on the unit interval. We divide this up into N subintervals, #and see which ones are visited by a long orbit N = 1000 dx = 1./N count = [0 for i in range(N)] # count is initialized to zero orbit = m.orbit(.5,1000000) #Assign points to their bins for x in orbit: i = int(x/dx) if (i>=0)&(i