Physics 213 Friday 8-29-03 IMPORTANT POINTS FROM PREVIOUS LECTURE: - Definition of nonlinear dynamics: properties of systems that evolve in time. - Definition of a map: a dynamical system with discrete time values, e.g., the logistic map. Definition of a flow, evolution equations that involve a continuous time. - A nonlinear dynamical system is nonlinear in the sense that the rule used to update the state of the system is nonlinear. Thus in a map of the form: x[n+1] = f( x[n] ) , the function f(x) is nonlinear if it does not satisfy the superposition property: f( c1 x1 + c2 x2 ) = c1 f(x1) + c2 f(x2) , where c1, c2 are arbitrary real numbers and x1, x2 are arbitrary real numbers. More generally, the updating rule may take some set of objects in a vector space V1 (numbers, functions, vectors, matrices, matrices of functions, etc) and change them into some other set of objects that form some other vector space V2. Then the above condition is then the same condition for linearity atd again c1, c2 are arbitrary real numbers, but now x1,x2 are elements of some vector space V1, and f(x1), f(x2), and f(c1 x1 + c2 x2) are elements of a second vector space V2. Challenge: is the rule that associates a 2x2 matrix M with its trace tr(M) a nonlinear rule? This maps the vector space of 2x2 matrices to the vector space of real numbers. What about the rule that associates the determinant det(M) of a matrix M with M? THE LOGISTIC MAP - Based on our discussions in class, I hope you are convinced that one of the simplest possible choices for the function g(N) would be a quadratic map, which is a _nonlinear_ function of the variable N: g(N) = [ r (1 - N/K) ] N , = a N + b N^2 , for constants r, K or constants a=r and b=-r/K. I have written the logistic map in this strange form to emphasize how a biologist would think of this model. Here nonlinear has a special meaning in the context of dynamical systems: the rule for evaluating g is not linear in that it does not obey a superposition principle. This in turn means that the evolution equation itself does not obey a superposition principal, the sum of two solutions is itself not a solution. CHALLENGE: Is the map N[t+1] = r N[t] + b for b != 0 a linear or nonlinear dynamical system? The multiplicative factor r is a reproduction rate since the simple linear model N[t+1] = r N[t] , can be analytically solved and has geometric growth (r > 1) or decay (r < 1): N[t] = r^t N[0] . Note: I use Mathematica notation for exponentiation here: r^t means the variable r raised to the t power. The quantity in brackets, [r (1 - N/K)], can then be thought of as a population-dependent reproduction rate since it acts as a multiplicative factor on number of insects N[t]. For small populations, N << K, this factor looks like the constant r by itself, but for larger populations approaching the CARRYING CAPACITY K of the population, this factor goes to zero and the population heads towards extinction. Most simple models of insect populations have this property: finite reproduction rate at low values of N[t], decreasing reproduction rate at high values of N[t]. The number K is an important property of the insect population and indicates the maximum number of individuals that can occur or the characteristic magnitude of insects expected. - The quadratic term has the important property of causing SATURATION of the exponential growth. Saturation is a hallmark of nonlinear systems, important but not the most interesting. - Note: Another way to derive the logistic equation is to assume (more realistically, we could hope) that a Taylor series to lowest order terms around the known value g(0)=0 might yield a useful model. I.e., we approximate the unknown function g(N) by g(N) = g(0) + g'(0) N + [g''(0)/2!] N^2 + h.o.t.s , where "h.o.t.s" means "higher-order terms", i.e., terms involving higher powers of the small quantity N. Since g(0)=0, we recover the logistic equation with quadratic term, with a=g'(0), b=g''(0)/2!. PLEASE REVIEW TAYLOR SERIES in your favorite introductory calculus book. They will be one of the fundmental tools we will be using over and over again in this course and that you will use over and over again after this course. HOW MANY PARAMETERS? NOT TWO BUT ONE - The dynamical system represented by the logistic model, N[t+1] = r (1 - N/K) N , appears to have TWO parameters r and K. Using an argument that we will use many times later in the course, I would like to argue that there is really only ONE effective parameter in this equation. We can see this by observing that the units with which we measure the magnitude of the population N[t] is arbitrary. A different choice simply scales (multiplies by a constant value) the present value by some constant C. So let's choose a new variable: x[t] = N[t] / N0 , which is dimensionless (no units) since we have chosen to normalize N[t] by some not-yet-specified unit of size N0 ("N zero"). We will then try to choose the constant N0 to eliminate a parameter (if possible). Substituting into the logistic equation gives: (N0 x[t+1]) = r (N0 x[t]) - (r / K) (N0 x[t])^2 . Dividing by N0 gives: x[t+1] = r x[t] - (r N0 / K) x[t]^2 , and we see that we can effectively remove one of the parameters, here K, by rescaling the unit of size to be the carrying capacity K, i.e., we choose N0 to have the value K N0 = K . The logistic equation then becomes: x[t+1] = r x[t] (1 - x[t]) . with a single parameter r and with x[t] a real number in the interval [0,1]. Any logistic equation with parameters r and K will have a solution that corresponds to a solution of this equation and vice versa. The big advantage of this rescaling is that we have a one-, rather than two-, dimensional parameter space to think about. CHALLENGES: - A particle of constant mass m is attached to a spring of spring constant k so Newton's equations of motion F=ma ("force equals mass times acceleration") give the equation: m (d^2 x / dt^2) = - k x . How many effective parameters are there in this equation? - How many effective parameters are there in the complex Ginzburg-Landau equation, a partial differential equation of importance in many areas of physics and engineering: tau u_t = eps u + D u_{xx} - E u^3 ? Here tau, eps, D, and E are constants, u_t denotes the partial derivative of u(t,x) w.r.t. the variable t, and u_{xx} denotes the second-order partial derivative of u(t,x) w.r.t. the variable x twice. EXAMPLES OF LOGISTIC MAP DYNAMICS - The Java demo http://www.expm.t.u-tokyo.ac.jp/~kanamaru/Chaos/e/Logits/logits.html provides a quick simple way to play with the dynamics of the logistic map. - The MMa function http://www.phy.duke.edu/~hsg/213/misc/plotlogistic.nb provides a more powerful way to explore the logistic map since it stores the generated numbers in a list that can be compared with other lists or anaylzed later. - Using Mathematica, showed briefly some example of dynamics that comes out of the logistic map. One sees decay toward 0 (extinction), toward a finite constant value, toward an oscillatory value, and highly irregular behavior as the reproduction rate r is varied. What was not well described was how the results may depend on the initial value x[0]. Some playing around will convince you that MOST initial conditions in the interval [0,1] lead to the same asymptotic behavior. TRANSIENTS - What happens if we start with an initial population x[0] that is close to but not exactly equal to the fixed point x=0 or the fixed point x=1-1/r? What happens depends on the value of the parameter r and is generally NOT easy to analyze except by numerical means. - What you see clearly in these plots are TRANSIENTS, there is some short-lived atypical behavior followed by statistically stationary behavior (fixed point, 2-cycle). It is an extremely important point that, for most nonlinear dynamical systems, it is the non-transient behavior that is most important. Questions: - How long does one have to wait for transients to decay? - What determines the transient time? - Do transients depend on initial population and on parameter r? NUMERICAL EXPLORATIONS OF THE LOGISTIC MAP - Let's use a computer code and Mathematica to take a quick and qualitative look at some of the dynamics of the logistic map, i.e., to see what some of its solutions look like as a function of time. A computer solution for the logistic map will consist of a finite set of N numbers x[0], x[1], x[2], ..., x[N] where each number is related to the previous one by the logistic map, i.e., given x[i] at time i, we form the product r*x[i]*(1-x[i]) and this gives us x[i+1]. So to simulate the logistic map on a computer, we need to specify THREE things: - an initial condition x[0] ; - a value for the parameter r ; - some details of how we do the computation, e.g., do we use single precision or double precision numbers (which have approximately 7 and 14 significant digits respectively), what computer we use, and what compiler we use. This last requirement is essential if you want someone else to be able to reproduce your results or to understand why they can't reproduce your results. - A little thinking suggests that we can't just plug in arbitrary values for the initial condition x[0] and for the parameter r. E.g., if we consider the logistic equation as a biological model, we need to have x[t] non-negative at all times since it is not meaningful to talk about a negative number of insects. But then each value x[i], including the starting value x[0] but satisfy 0 <= x[i] <= 1 , else x[i+1] will be negative (assuming r is positive). - Let's say we do choose x[0] in the interval [0,1] to start a sequence. Is there any guarantee that x[i] will remain in this interval many iterations down the road? Not necessarily but we can, in fact, choose the reproduction rate r to guarantee this property. No matter what value x[0] takes on in the interval [0,1], one can never get a next value x[1] that is greater than the maximum of the curve: f(x) = r x (1 - x) , on this interval. If we choose r sufficiently small that the maximum is never larger than 1, we are home free. Calculus (or graphics) tells us that the maximum occurs at x=1/2 with a value of f(1/2)=r/4. So if we restrict the value of r to lie in the following interval 0 <= r <= 4, AND if the initial condition x[0] is chosen in the interval [0,1], then the sequence x[i] generated by the logistic map will never leave the unit interval. This is nice. - WARNING: Our conclusion above, that r should lie in the interval [0,4] in order that x[i] stays within [0,1] forever, is correct mathematically but ignores a fundamental fact of digital computers, that all calculations are carried out with a finite precision and so every calculation can incur a small but nonzero error. As an example, for the value r=4, it is, in fact, possible for the function 4 x (1 - x) to evaluate to a number just slightly larger than 1 because of round-off errors, e.g., if x is close to, but not equal to, the value x=1/2, e.g., if x=0.500000001 or x=0.499999999. The next value x[i+1] will be negative and then the sequence is in trouble. CHALLENGE: What is the largest value of 4x(1-x) you can achieve for any floating point number x in the interval [0,1], i.e., by how much can you exceed the mathematical value of 1 at the maximum? NUMERICAL PROPERTIES OF THE LOGISTIC EQUATION x[i+1] = r x[i] (1 - x[i]) - Now that we have some idea of what range of values we can try for the initial value x0 and for the reproduction rate r, let's look at some actual time series. The following is a short Mathematica function called "plotlogistic" that creates and plots a sequence of numbers whose dynamics is generated by the logistic equation: (* MMa function for iterating and plotting the logistic equation *) plotlogistic[ r_ , (* reproductive rate *) x0_ , (* initial value *) n_ (* number of elements in sequence *) ] := Module[ { (* local variables *) i, xlist } , xlist = Table[ 0., {n} ] ; (* allocate space *) xlist[[1]] = x0 ; For[ i = 1 , i < n, i++ , xlist[[i+1]] = r xlist[[i]] (1. - xlist[[i]] ) ] ; ListPlot[ xlist, PlotJoined -> True, (* connect points with lines *) PlotRange -> {0.,1.} (* vertical range [0,1] *) ] ; Return[ xlist ] ] (* end of plotlogistic *) There are numerous features of the Mathematica language that you will need to learn on your own during this course. I will discuss a few of them in lecture and give a few pointers here. As usual for learning any new computer language, it is often best to start with a working code and to modify it, rather than to write a program from scratch. Some comments on the MMa function plotlogistic: - The "blank" notation "_" that indicates a variable of a function, e.g., "fill in the blanks..." - The use of comments, i.e., text between the symbols (* and *). Anything between these symbols is ignored by Mathematica. - The use of delayed evaluation ":=" when defining a function. It is QUITE important that you use delayed evaluation, not the equals sign "=", when defining a MMa function for reasons that are fundamental to any computer algebra program. - The use of "Module" to group statements together. This is analogous to the body of a function, e.g., the braces {...} in a C code. - The collection of LOCAL variables in a list at the beginning of a Module. Any variable that you use that is not in this list is visible outside the function and retains its value after you call this function. When writing MMa code, you often want to put most of your Module variables in a local list so that you don't accidently collide with a previously defined variable of the same name. - The concept of "lists" which is the MMa way of grouping numbers together as a vector. Lists of lists are used to represent a matrix. - The use of the Table function as a widely used trick to generate a list of a given length. This is roughly similar to "malloc" in C. - The double-bracket notation [[1]] to indicate the index of a list. Thus xlist[[1]] is the first element of a list name xlist. A single bracket, e.g., f[x], denotes function evaluation of a function f at value x. And parentheses are used purely for grouping, e.g, "a(2)" is not the function a evaluated with argument 2 but the symbol "a" multiplied by the quantity 2, i.e., the symbol '2a'. This notation can be confusing at first. - The For statement in MMa works just like the "for" statement in C and C++ except that you have to be careful about using commas versus semicolons. - In a Module statement, all successive statements are separated by a semicolon. A common and frustrating mistake in Mathematica is to leave out a semicolon by accident. Mathematica will then try to multiply one expression by the next, usually giving a cryptic error message. (Try leaving out a semicolon intentionally in this program to see what happens.) - The use of ListPlot to plot a list of numbers. Like many MMA functions, one can add "options" which modify the plot. In the above, I use two options, PlotJoined -> True to connect successive points, and PlotRange -> {0.,1.} to make the vertical axis of the plot span the range 0 to 1. You should leave out each of these statements in turn and try using plotlogistic to see the effect of these options. By the way, the arrow notation "->" means "substitute" or "replace". The ListPlot function has internal variables called PlotJoined and PlotRange and the options are a simple way of replacing these variables with the desired values. - Instead of looping with a For statement, one could use the NestList function. I find the For statement more readable. - You can make the points in ListPlot bigger by using the option: Prolog -> AbsolutePointSize[2] - WARNING: The function plotlogistic.m is dangerous to use because I haven't added any tests to make sure that the arguments r, x0, and n have sane values, e.g., that they evaluate to numbers, that r is positive and between 0 and 4, that x0 lies in [0,1], and that n is a positive integer. Good software practice is to make sure that all input parameters are validated thoroughly BEFORE you begin any calculations, with invalid data being reported and the code terminated. - We can finally try out our code with a few examples, e.g., try executing the following: plotlogistic[ 0.9, 0.1, 50 ] (* try it out *) plotlogistic[ 0.9, 0.9, 50 ] We see decay toward extinction, x[i] -> 0, for most initial conditions. plotlogistic[ 2.8, 0.1, 50 ] (* try it out *) plotlogistic[ 2.8, 0.9, 50 ] We see growth, a few oscillations, and asymptotically a steady state, x[i] -> constant > 0. plotlogistic[ 3.2, 0.1, 50 ] (* try it out *) We see oscillatory growth but now the oscillations don't die out. The oscillations have period 2, i.e., the values repeat every other year. plotlogistic[ 3.9, 0.1, 50 ] Whoa, crazy behavior here! What is going on? Can insect populations really bounce around like this? Please explore some other choices of reproduction rate r and initial value x[0] so that you can build up your own intuition about the dynamics of the logistic map. CHALLENGE: What happens if you type something crazy on purpose, e.g., plotlogistic[ 0.9, 50 ] (* one argument left out *) plotlogistic[ 0.9, 0.1, 50., 4 ] (* too many arguments *) plotlogistit[ a0, 0.3, 50 ] (* symbolic argument a0 *) Can you understand what MMa does and why? BIFURCATION DIAGRAMS - Make sure you understand how Figure 10.2.7 on page 357 of Strogatz was created and how to interpret this plot. - You should look through the MMa code bifurcationdiagram.nb (on the course web page via the Miscellaneous Files link http://www.phy.duke.edu/~hsg/213/misc/ which gives one example of how you would create a diagram of this kind. Notice the subtlety that one has to specify a time for transients to decay an acceptable amount. Further, parts of the diagram have hysteresis so some details of the figure depend on which initial conditions were used to start the iterations for a fixed value of the parameter r. Can you find example of coexisting solutions? CHAOS - Bounded nonperiodic behavior. This can result when all fixed points and all periodic orbits are unstable (a sufficient but not necessary condition). Any initial condition not starting exactly on a fixed point or starting exactly on a periodic orbit will be unable to converge towards one of these unstable objects and the resulting dynamics is irregular and random looking. We will develop tools later in the course to quantify this irregularity. A numerical comment: it is often the case that the numerical value corresponding to a fixed point or to a 2-cycle corresponds to a number that can not be represented exactly in the set of machines numbers of a computer (quite roughly, all numbers of the form: 1.d1 d2 d3 ... d52 X 2^m , where the mantissa "1.d1..." is a number with binary digits (each di is either 0 or 1) and where m is a positive, zero, or negative exponent determined by the precision being used.) Thus any irrational number such as the Sqrt[2] or cube root of 3 requires an infinity of binary or decimal digits and can only be approximated on a digital computer with a finite number of bits of information. This means that, in many cases, it is IMPOSSIBLE to represent a fixed exactly in the computer and ALL numerical initial values will move away from the fixed point. SONOFICATION OF THE LOGISTIC MAP - If your computer has a sound card, you might enjoy LISTENING instead of seeing the intricate bifurcation diagram of the logistic map. Type in and execute the following code and see if you can understand it: logistic[ n_Integer ] := Module[ {f, t, x} , f = Compile[ {x,t}, Evaluate[ (3 + t/n) x (1 - x) ] ] ; FoldList[ f, 0.233, Range[n] ] ] ListPlay[ logistic[8000], SampleRate -> 2000 ] This is again only approximate because of finite observation times of infinitely-long lived transients and a finite sampling rate.