Physics 213 9-2-03 QUESTIONS TO THINK ABOUT - The bifurcation diagram for the logistic map that we discussed in class makes a remarkable prediction: that an infinity of bifurcations can occur over a finite range of values for the parameter. Do you think it would be possible to see experimentally an infinity of bifurcations? If not, why not? - Can you find an example for the logistic map for which there are two distinct attractors, i.e., two different non-transient behaviors that are both stable? If this is possible, how would you show this on the bifurcation diagram. ANY SIMPLE SOLUTIONS? FIXED POINTS AND CYCLES - You should NEVER trust any numerical computer code until you have confirmed its behavior for some analytically known solutions. Since plotlogistic is our first Mathematica code (at least for this course), how do we know anything is working correctly? Are there some solutions that we can solve for analytically to compare with the plotlogistic plots? Can we understand at all why there are different kinds of non-transient solutions? - Answer is yes, there are many different kinds of solutions that we can work out mathematically (rigorously) and then try to compare with the numerical calculation or just try to understand conceptually and mathematically. The simplest class of solutions are ones that are constant for all time, x[i]=constant. Since x[i] is unchanged under the map function g(x), constant solutions are called FIXED POINTS since their values are fixed under the dynamics. They are specifically called fixed points, not just constant solutions, since in the space of possible states (phase space as we will call it), a constant solution is just a geotmetric point. Finding fixed points of a map requires one to find the zeros of a nonlinear function. Thus if we have a general map x[i+1] = g( x[i] ) , and we try to find a constant solution x[i] = x*, then the value x* must satisfy: x* = g( x*) , i.e., x* is a root (zero) of the function x - g(x) = 0. - Let's see how this works for the logistic equation with g(x) = r x (1 - x). Then any fixed point of the dynamics must satisfy: x* = g(x*) = r x* (1 - x*) , which has two solutions, x* = 0 extinction x* = 1 - 1 / r . finite constant population The second fixed point x* is negative and therefore unphysical if the positive parameter r is less than one. Actually, this fixed point is inaccessible and so doesn't really exist since if we start with a nonnegative initial condition x[0] in [0,1], and r in the interval [0,4], we can never approach a negative number (if we ignore the possibility of floating point errors as noted above). So the second fixed point doesn't exist until r >= 1 . Condition for nonzero positive fixed point. So our discussion has already shown something new, something special happens at r=1, a new solution becomes available. The appearance or disapperance of a solution as a parameter is varied (here r) is one of the simplest examples of a bifurcation of the dynamical system. WARNING: FINDING FIXED POINTS CAN BE VERY VERY HARD - There is a simple geometric way to identify where fixed points occur for a one-dimensional map. A constant solution x[i]=x* of the map x[i+1] = g( x[i] ) must satisfy x* = g( x* ) . We can think of this as the intersection of the 45 degree line y=x with the plot of the updating rule y=g(x), i.e., any place where the function g(x) hits the 45 line is a fixed point. Further, the logistic map is simple enough that we can find the two fixed points analytically since everyone knows how to find the roots of a quadratic equation. - I don't emphasize this geometric picture too much because it is highly misleading about how easy it is to find fixed points, especially for cases we discuss later in the course that involve two- or higher-dimensional state vectors. Even for general real functions g(x) of a real variable x, we can NOT find any solutions analytically and one must proceed numerically to find APPROXIMATE solutions. Numerical methods for finding zeros all work by generating a sequence of better and better approximations that converge in the infinite limit to a zero. We will discuss this point later on when we discuss Newton's method but, for now, you should read about the Mathematica function FindRoot for finding approximate roots by numerical methods. This function requires that you have an initial guess as to the position of a root, already a difficult problem. A typical invocation would look like this: FindRoot[ Cos[x] == x, {x, 1} ] Mathematica can find symbolic analytical answers, e.g., see what happens with: Solve[ x == r x (1 - x), x] , but there are few functions g(x) for which this will work. CHALLENGE: Use FindRoot as an alternative way to find the two fixed points of the logistic equation and make sure that you understand the way that Mathematica returns the answers as a list of substitutions. Since FindRoot locates approximate numerical values, figure out how accurately are the answers that FindRoot locates. Later in the course, we will need to find fixed points of maps with many variables. This will require us to solve systems of nonlinear equations of the form F(X)=0 where F is a N-dimensional vector function of a N-dimensional vector X and the right side is a vector of N zeros. (For some numerical problems involving partial differential equations, N may be as large a 10,000,000!) This can be VERY hard; there are no known systematic ways to find zeros of higher-dimensional nonlinear sets of equations. Even the best numerical methods become subtle and unpredictable for such large systems of nonlinear equations. At this point, laboratory experiment and simulations can be invaluable for suggesting where to look for solutions. CALCULATING CYCLES OF A MAP - Another class of simple solutions is the class of purely periodic ones, i.e, sequences that repeat themselves exactly after some period T. Such sequences are called "cycles" since they lead to cyclic behavior. Since maps take integer time steps, any period of a cycle must itself be an integer and so we talk about 2-cycles, 3-cycles or 100-cycles that repeat exactly after 2,3, and 100 time steps respectively. CHALLENGE: I meaningful to think about a fixed point as a cycle? If so, what would its period T be? Is any 2 cycle automatically a 4-cycle and an 8-cycle? - Does the logistic equation have any cycles and can we compute them? Let's start by looking for 2-cycles, which would be two numbers p and q that alternate under the dynamics: p = g(q) , q = g(p) , and so on. This will generate the endless sequence characteristic of a 2-cycle: p, q, p, q, p, q, p, q, ... Note that if we substitute the first equation, p=g(q), into the second equation, we get a fixed point problem for q: q = g(g(q)) , and similarly if we substitute the second equation into the first, we get the same fixed point problem for p: p = g(g(p)) . This is a general result: any member of a 2-cycle must be a solution of the function g(x) composed with itself twice: x = g(g(x)) . GENERAL EQUATION FOR A 2-cycle WARNING: Although any member of a 2-cycle satisfies this equation, not all solutions of this equation is a member of a 2-cycle. E.g., any fixed point x=g(x) automatically satisfies this equation. So when you solve the 2-cycle equation, you need to be careful to throw out the fixed points. - For the logistic map, g(g(x)) is a 4th-order polynomial and to find all the 2-cycles, we need to find the solutions to the following equation y = r^2 y (1 - y) (1 - r(1 - y)y) . A fundamental mathematical theorem tells us that a polynomial of degree n has exactly n roots (some of which may be complex-valued) so for this 4th-order polynomial, we expect exactly four roots. One root is obviously y=0 and the other three roots then satisfy a cubic polynomial. It is somewhat tedious to continue by hand, but straightforward using Mathematica, which can solve many problems symbolically and analytically. Try typing g[ x_ ] := r x (1 - x) (* Define logistic map *) Solve[ y == g[ g[y] ] , y ] (* Solve for member of 2-cycle *) This second statement says in English "solve the equation y=g(g(y0) for the variable y". Mathematica will produce an answer in the form of a list of 4 lists: { {y -> 0}, {y -> 1-1/r}, {y -> something messy}, {y -> something else that is messy} } MMA INTERJECTION: You should read the section of the online Mathematica book (available under the Help menu), which explains the notation "->" which is what MMa calls a "rule". One applies a rule to a given symbolic expression using the "replacement" notation "/." which basically means replace all symbols to the left of the rule by the symbols on the right of the rule. Thus if we want to see what the value is of the 4th solution, we would type the following: y /. answer[[4]][[1]] which means to pick off the first element of the list answer[[4]], which itself is a list occupying the 4th position of the list answer. The notation answer[[4,1]] can also be used to cut down on the number of brackets that you need to type. Please try typing these expressions to make sure that you understand this notation. Two of these solutions are the fixed points that we found previously and so the other two numbers must be the members of a 2-cycle. We conclude: THERE IS THEN AT MOST ONE 2-CYCLE OF THE LOGISTIC MAP FOR ANY CHOICE OF THE PARAMETER R. If you study the algebraic expression of the messy things, you will see that they both involve the expression: Sqrt[-3 - 2r + r^2] = Sqrt[(1 + r)(r - 3)] These expressions are then complex-valued unless r >= 3 and so we can make two further predictions (this is becoming fun, we can make predictions): The logistic map has no real-valued (physical) 2-cycles unless the parameter r is increased to at least the value 3. There is only ONE 2-cycle for all values of the parameter r >= 3. You can verify that these analytic predictions seem to be correctly by using the MMa fn plotlogistic for different values of r and different initial conditions. CHALLENGE: Figure out in Mathematica how to plot the two values p(r) and q(r) (that make up the 2-cycle) together as a function of the parameter r and confirm that no solutions exist until r >= 3. CHALLENGE: Write down the equation that any member of a 3-cycle must obey. Using FindRoot, determine the smallest value of r for which a 3-cycle will exist (all three members have real values). Is there only one 3-cycle for any value of r, as was the case for a 2-cycle? Similarly, work out the smallest value of r for which a 4-cycle will exist. Interesting question: can you identify a pattern to the order in which cycles of different periods come into existence as the parameter r is increased from 3 to larger values? We see that a 2-cycle turns into a 4-cycle which turns into an 8-cycle? When does a 3-cycle first appear? A 5-cycle? A p-cycle for p a prime number? COMMENT: GALOIS THEORY AND FINDING ANALYTICAL ROOTS - For the logistic equation, finding cycles of higher and higher period require finding roots of higher- and higher-order polynomials. There is a deep and famous (and old) result due to Evariste Galois who showed (at the age of 17!) that it is generally impossible to find the roots of a polynomial of degree higher than 4 as algebraic functions of the coefficients. Thus no formula similar to the quadratic formula exists in general for polys of order 5 or greater. This pure mathematical result has extremely important practical consequences for us since it implies that we have no choice but to use numerical methods for finding roots of higher-order polynomials. So you better learn how to use FindRoot as a general strategy of finding fixed points and cycles! - Mention Routh's theorem: possible to say something about stability without knowing the roots explicitly so Galois's theorem is not quite so negative. BIFURCATIONS - Bifurcations are qualitative changes in the dynamics. When a fixed point "changes into" a 2-cycle or a 2-cycle "changes into" a 4-cycle or a fixed point replaced by chaos, these represent bifurcations as some parameter is varied since the non-transient dynamics changes. What is really happening is a loss of stability. The fixed point that changes into a 2-cycle is actually still present but is unstable in the sense that most initial values that start sufficiently close to the fixed point will move away from the fixed point to some other solution. If one chooses an initial value without special knowledge (without knowing the exact value of the fixed point, here x* = 1 - 1/r), then that initial condition will not decays towards the fixed point and you will not observe the fixed point as you iterate the logistic map. For most laboratory experiments and numerical simulations, just one parameter is varied at a time (e.g., the temperature of a system or a voltage across the system). This parameter is often called the "BIFURCATION PARAMETER" since it can lead to bifurcations in the dynamics. For the logistic map, the parameter "r" would be the bifurcation parameter. A reminder from an earlier lecture: recall that there are two important kinds of bifurcations: continuous ones (the lingo is "supercritical" or "forward" bifurcations) in which the new behavior grows smoothly out of the old behavior, and discontinuous ones (the lingo is "subcritical" or "backwards" bifurcations). We will be able to classify most of the common bifurcations. LINEAR STABILITY - When certain kinds of dynamics have been identified, e.g., fixed points or 2-cycles, it is often useful to ask whether these dynamics are STABLE or UNSTABLE. We will give a careful mathematical definition of these concepts later on in the course, but a casual meaning of STABLE dynamics is that ARBITRARY sufficiently small (infinitesimal) perturbations of the dynamics lead to the SAME nontransient dynamics. A dynamics is unstable if there is at least one small perturbation that causes the dynamics to change. Since most physical systems are subject to noise, e.g., from molecular collisions, if there is at least one perturbation that leads away from a fixed point or other dynamical behavior, nature often finds it rapidly because of noise. This concept of stability is familiar to many of you in day to day life. Thus a marble placed at the top of an inverted bowl will roll down the side at the slightest touch so the top of the bowl is an unstable fixed point of the motion of the marble. A marble placed at the bottom of a normal bowl will return to the bottom (and remain stationary) if kicked a little bit and if you wait for friction to cause the transient to die out. Clearly one should consider only small perturbations since a large perturbation could kick the marble out of the bowl entirely, which is a new situation not possible to analyze in terms of the geometry of the bowl. For a simple map like the logistic map, a fixed point x* is stable if we can perturb the fixed point by an arbitrary small amount, say choosing an initial condition: x[0] = x* + eps , where eps is a small number (either positive or negative), and if the resulting sequence of numbers x[0], x[1], x[2], ... converges back to the fixed point, i.e., x[i] -> x* as i goes to infinity. On the other hand, if there is at least ONE arbitrarily small value of eps such that x[0] = x* + eps generates a sequence that leads away from the value x*, then the fixed point is unstable. - It turns out that there is a simple precise mathematical criterion for when a fixed point is stable or unstable as some parameter p is varied. We will discuss the general vector case later in the course (for maps or flows involving N variables), but can just give the criterion for a fixed point x* of a general map x[n+1] = g(x[n]): linear stability of x* if | g'(x*) | < 1 , linear instability of x* if | g'(x*) | > 1 , indeterminate linear stability if | g'(x*) | = 1 . FORESHADOWING: For a vector map X[n+1] = G(X[n]), where X[n] is a N-dimensional vector and G(X) is a N-dimensional vector field, the stability of a vector fixed point X* that satisfies X*=G(X*) is determined by whether all N eigenvalues e[i] of the NxN Jacobian matrix dG/dX[X*] (evaluated at the fixed point X*) have magnitude strictly less than one: |e[i]| < 1 for all i between 1 and N. We will derive this important result later in the course. This general condition reduces correctly to the above criterion for a one-dimensional map since the 1x1 Jacobian matrix is, indeed, the number g'(x*) which is then also the eigenvalue of the matrix. EXAMPLES OF LINEAR STABILITY - Let's first discuss some applications of the concept of linear stability and then later derive the above criteria. - When does the fixed point x=0 become unstable for the logistic map? g'(x) = r - 2 r x which has the magnitude | g'(0) | = r , for x*=0. Since we assume from the start that r >= 0, the extinction fixed point is linearly stable when r < 1, unstable for r > 1, and marginally or neutrally stable for r=1. Since biologists are mainly interested in cases when species are not going extinct, researchers often just consider the logistic equation for r >= 1, i.e., for the parameter range 1 <= r <= 4 . Recall that we want r <= 4 else the quadratic function r*x*(1-x) can produce a value x[n+1] larger than 1 for some value x[n] in [0,1], leading to a negative value x[n+2] and a diverging sequence. - When does the only other fixed point x*=1-1/r become unstable? We calculate: | g'(1-1/r) | = | 2 - r | > 1 , which occurs for r < 1 , which is when the extinction fixed point is stable r > 3 . The fixed point 1-1/r is marginally or neutrally stable when r=1 or r=3 and a nonlinear analysis is required for further insight, that takes into account the next higher-order terms in the Taylor-series expansion about the fixed point. - Let's calculate the stability of a 2-cycle p,q,p,q,p,... of the map x[n+1]=g(x[n]) , a member p of which obeys the fixed-point equation p = f(p) , where f(x) is the function given by g(g(x)). Then the criterion for instability of a 2-cycle is that: 1 < | f'(p) | = | g'(g(p)) g'(p) | = | g'(q) g'(p) | , i.e., the product of the derivatives of g(x) evaluated at the two members of the 2-cycle must have magnitude less than one. Let's check this out for the logistic equation x[n+1] = r x[n] (1 - x[n]) , using Mathematica and ask: when does a 2-cycle become linearly unstable, so that we could expect some new kind of dynamics? The MMa statements: g[ x_] := r x (1 - x) ; answer = Solve[ x == g[g[x]], x ] ; answer (* show the values of answer *) p = x /. answer[[3,1]] ; (* one member of cycle *) q = x /. answer[[4,1]] ; (* other member of cycle *) Plot[ {p, q} , {r,3,4} ] ; (* Plot 2-cycle values *) (* |d/dx g(g(x))| = |g'(p) g'(q)| *) exp1 = Abs[ D[ g[g[x]], x ] /. x -> p ] // Simplify ; Plot[ exp1 - 1 , {r,3,4} ] (* see where slope is near one *) FindRoot[ exp1 == 1, {r, {3.4, 3.5}} ] which yields the approximate numerical value r -> 3.44949 as the value of neutral stability of the 2-cycle. We can confirm that this value is close to a root by evaluation: exp1 - 1 /. r -> 3.44949 which gives the value 1.3 10^-6, fairly close to zero. Strogatz, on pages 360-361 of his book, shows that one can obtain the answer analytically which turns out to be 1+Sqrt[6], confirming the numerical result above (1+Sqrt[6] ~=3.449489743...). I have shown you how to get the result numerically because that is the best that you can do in most cases. Still, when analytical methods are possible it is VERY important to use them since you then have precise knowledge of what is going on and often gain extra intuition. We then expect the 2-cycle to become linearly stable for r < 3.44949 and linearly unstable for r > 3.44949 which you can confirm is indeed the case by iterating the logistic map near these values, e.g., plotlogistic[ 3.40, 0.45, 50 ] (* yes, a 2-cycle! *) plotlogistic[ 3.50, 0.45, 50 ] (* a 4-cycle *) where plotlogistic is the function we used in a previous lecture. Note that for r=3.44949, p and q have the approximate values 0.43996 and 0.849938 so the two above plots with initial values 0.45 are starting off close to a 2-cycle. WARNING: when you are so weakly unstable (close to a bifurcation point), the growth rates of perturbations are small (multiplicative factor g'(x*) is close to one in magnitude) and you have to wait longer for transients to die out! CHALLENGE: Repeat this analysis and determine the value of the parameter r at which a 3-cycle first appears, how many 3-cycles exist for a given value of r, for what value of r the 3-cycles all become unstable, and what kind of dynamics happens at that point? Some of this is discussed in Strogatz, pages 361-363, but the results are fairly easy to establish numerically using some of the above ideas. CONCEPTUAL CHALLENGE: Can one talk about chaos becoming stable or unstable as some parameter is varied? Is stability only well-defined for fixed points or cycles? IMPORTANCE OF LINEAR STABILITY THEORY - Linear stability theory, the analysis of the growth or decay of infinitesimal perturbations around some given dynamical state, is an extremely important and practical part of nonlinear dynamics, much more important than some of the more exotic topics that we will come across later in the course such as strange attractors and fractals. Linear stability theory is used heavily by scientists and engineers in all disciplines to predict when some familiar behavior may change into something undesirable and what some of the properties may be of the new behavior. However please keep in mind that a linear stability analysis only tells us for what parameter values a given state becomes unstable, it is a LOCAL analysis. It is not capable of predicting what kind of dynamics the perturbation will evolve to. Thus when a fixed point becomes unstable, the non-transient dynamics could be another fixed point or a periodic behavior or chaos or something else. DERIVATION OF LINEAR STABILITY CRITERION FOR ONE-VARIABLE MAPS - Let's derive our criterion for linear stability so that we can see precisely what kinds of assumptions and arguments are used. We will repeat these ideas later for flows and for maps with N variables, the general case. Let's assume that we have a map x[n+1]=g(x[n])) with fixed point x* and let's try to inquire how a starting point x[0] evolves under the map that starts off extremely close to the fixed point. We will write our starting point as x[0] = x* + eps , where eps is some tiny nonzero number, |eps| << |x*|. Here is a key insight: by continuity of the function g(x), we can see that the smaller the perturbation eps, the longer that successive iterates x[1], x[2], etc will stay close to the fixed points. It is then useful to write ALL the members of the orbit growing out of x[0] as perturbations around the fixed point so let's write: x[i] = x* + y[i] , where y[0] = eps and we guess that y[i] will be tiny for many iterations provided that eps is sufficiently small. We can then write the equation: x[n+1] = g( x[n] ) , in terms of the new variable y[i] that we consider small: x* + y[i+1] = g( x* + y[i] ) . Here comes the second crucial step in our argument. If y[i] is small compared to the magnitude of x*, we can use a Taylor series to expand the function g(x* + y[i]) around the value x*: g(x* + y[i]) = g(x*) + g'(x*) y[i] + h.o.t.s where h.o.t.s denotes "Higher Order TermS", which would look like this: (1/2!) g''(x*) y[i]^2 + (1/3!) g'''(x*) y[i]^3 + ... If y[i] is tiny, y[i]^2 is tinier still and y[i]^3 is even smaller and to good approximation we can ignore all terms that are higher order than the first in y[i]. This corresponds to replacing the function g(x* + y[i]) by a function linear in y[i] and is the origin of the term "linear stability". So ignoring these higher-order terms in the Taylor series, we get: x* + y[i+1] ~= g(x*) + g'(x*) y[i] , where the notation ~= means "approximately equal to". But now we know that g(x*)=x* since, by definition, x* is a fixed point of the equation x[n+1]=g(x[n])). We can then cancel the x* and g(x*) from both sides and end up with a new equation just involving y[i]: y[i+1] = g'(x*) y[i] , which holds as long as y[i] is so tiny that we can ignore higher-order powers of y[i]. But we know how to solve this linear map, the solution is: y[i] = ( g'(x*) )^i y[0] = ( g'(x*) )^i eps , and the ith perturbation will decay if the constant multiplicative factor has magnitude less than one or will grow if it has magnitude greater than one, i.e., there is linear stability if: | g'(x*) | < 1 , and linear instability otherwise. For the rare case that g'(x*) has magnitude one exactly (so value g'=1 or g'=-1), it looks like y[i] will not grow or decay but this conclusion is wrong since we then have to include higher-order terms in the Taylor series to resolve this situation. In any case, we have succeeded in deriving our criterion for stability and instability of a fixed point. The crucial assumption is that the perturbation around the fixed fixed point is sufficiently small that we can use a Taylor series expansion to lowest (linear) order in the perturbation. NOTE: If g'(x*) < -1, i.e., the multiplicative factor is negative and larger than one in magnitude, then one gets an oscillatory instability since the perturbation changes sign as it grows. This explains a detail that we saw in at earlier Mathematica calculation: for r > 1, the fixed point x* = 1 - 1/r , CHALLENGE: For dynamics in a bounded domain, e.g., the variable x in the interval [0,1] for the logistic map, how should stability be defined for points like x=0 or x=1 that occur at the border of the domain? WHY THE LINEAR-STABILITY CRITERION APPLIES TO CYCLES - We noted that a member p of a 2-cycle of a map x[n+1]=g(x[n]) will satisfy the equation p = g( g(p) ) , but it is NOT immediately obvious why we can apply the linear stability criterion to g(g(x) since this is NOT the equation that determines the dynamics, which is given by x[n+1] = g(x[n]). So how was I able to get away with this in our discussion above? It should be clear to you that from a map x[n+1]=g(x[n]), I can derive a NEW map that connects every other point of the time series x[n]: x[n+2] = g( g( x[n] ) ) . I can then pretend that this is the only dynamical system that I am interested in and can ask the question what are the fixed points of this map and what determines their linear stability. Hopefully you can see that linear stability or instability for this map implies linear stability or instability of the 2-cycle in the original map. CHALLENGE: Give a careful proof of this fact. CHALLENGE: CHECK OF YOUR UNDERSTANDING OF LINEAR STABILITY THEORY - Assume that you have a one-variable map x[n+1]=g(x[n]) and assume that this map has a fixed point x* that, for given parameter values, is weakly unstable with the specific value: g'(x*) = -1.001 which has magnitude just barely larger than 1. Assume that you choose as an initial condition a value close to the fixed point: x[0] = x* + 0.01 Without using a calculator, what is the value of x[20]? Solution: The perturbation to the fixed point is dx[0]=0.01 and linear stability theory tells us how this perturbation grows under successive iterations of the map: dx[i] = [ g'(x*) ]^i dx[0] We have i=20, g'(x*) = -1.001 = - (1 + eps) where eps=0.001 and so: dx[20] ~= (-1.001)^20 dx[0] ~= (-1)^20 (1 + 0.001)^20 0.01 By Newton's binomial theorem, (1 + eps)^20 ~= 1 + 20 eps + h.o.t.s = 1 + 20 (0.001) = 1.02 and so: x[20] = x* + dx[20] = x* + 0.00102 After 20 iterations, one has moved only a little bit further, 2%, from the fixed point. This confirms also that a linearization is still a good approximation for this many iterations.