Methods / ti-cas.org
	The TI-92 at the edge of chaos

O The MODEL OF VERSHUSLT

About 1845, the Belgian mathematician Pierre Verhulst proposed the following model to describe the population evolution in a given environment, by supposing that this population admits an ideal manpower E.

By noting P₀ the initial manpower, and P_n the observed manpower at date n, one supposes that the quantity (P_n+1-P_n) /P_n (the rate of variation between year n and the year n+1) remains proportional the rate of variation between year n and the year n+1) remains proportional to E-P_n (i-e what's separate current manpower from ideal manpower). While noting K this proportionality factor (positive), we arrive at the formula: P_n+1=P_n+K*P_n*(E-P_n).

This model isn't at all eccentric. It is noted for example that if P_n=E, then for very m>n, P_m=E. Ideal manpower is thus a «stable» situation.

In addition, if P_n<E, P_n+1>P_n, and if P_n>E, P_n+1<P_n: evolution is always made «in direction» of the ideal manpower E. The increase rate (P_n+1-P_n) /P_n, equal to K*(E-P_n), is all the more high since one are «far» of ideal manpower , which seems still rather natural there.

The problem comes from K, which measures the sensitivity system. Larger K is ... more the population reacts in a violent «way » to a variation with the ideal population. In particular, it is thus possible that one passes from a Pn<E manpower to Pn+1>E, and so on...

The problem is simple: Does the population evolve wisely worms ideal manpower , and if not what's occur ?

O REWRITING OF THE PROBLEM

It is simpler to consider the quantities U_n=P_n /E (measuring the relationship between the current population and ideal manpower). We obtain the relation by posing k=K*E : u_n+1=u_n+k*u_n(1-u_n).

Here, we can consider that ideal manpower is 1. It will be advisable to give itself an initial manpower u0 positive (u0<1 if one leaves of a lower manpower than the ideal manpower). With this problem translation, we are in a known ground : it is a question of studying a sequence obeying at a recurrence un+1=f(un), where f is defined by: f(x)=x+k*x*(1-x). The SEQUENCE mode of the TI92 will be very useful because it will enable us to visualize the behavior of the sequence (un).

NOTE : it is shown that this type of problem can still be simplified. One could thus bring back oneself, by a change of variable, with the recurrence un+1=c*un*(1-un) or with the recurrence un+1=un2+c. By convenience, we will remain about it however with the initial formulation.

O PRELIMINAIRES

We thus note f(x)=x+k*x*(1-x), where k is a strictly positive constant, and where x is a positive variable. It is checked that F is increasing afterwards of x=0 to x=a_k=(k+1)/(2k) then decreasing. Its maximum is f(a_k)=(k+1)²/(4k).

It is cancelled into 0, is initially positive, then is cancelled at b_k=(k+1)/k=2a_k.

It is important that the interval I=[0,bk ] (that on which f(x)³0) is stable by f (i.e. f(I)ÍI). That ensures that if one chooses u0 in I, then for any n, un is in I.

To say that I is stable by f, it is to say that f(ak) (the maximum of F on I) is lower than bk. That is written (k+1)²/(4k) £ (k+1)/k, and that is equivalent saying that: k£3. Morality: subsequently, the coefficient k will check 0<k£3, and the initial term u0 will be such as 0£u₀£1+1/k.

Under these conditions, one is certain that, for any n, 0£u_n£1+1/k.

In the study of the sequence defined by u0 and by u_n+1=f(u_n, the graph of f brings many indications, in particular by its placement compared to the line y=x. Two meet with points (0,0) and (1,1), which provides two possible limits for the sequence (un), 0 and 1.

Let us note that 0 cannot be a limi of the sequence (un) (except some exceptional cases where all the un are null starting from a certain row). Intuitively, that comes owing to the fact that the modifications of «population » are always made «in direction » of the ideal population (here 1). One can also notice that if un is close to 0, then: u_n+1=u_n+k*u_n(1-u_n)@(1+k)u_n>u_n.

Another more traditional explanation: a fixed point a of a recurrence vn+1=g(vn) (i.e. a solution of g(x)=x, or a possible limit of the sequence (vn)) is "gravitational" if -1<g ' (a)<1 (all the more gravitational that g'(A) is close to 0) and «repulsive » if |g ' (a)|>1. In our case, f '(0)=1+k>1: the fixed point 0 is thus of repulsive nature.

Of the same f '(1)=1-k<1. This fixed point will be thus «gravitational »if 1-k>-1 , so if k<2 (i.e. if the sensitivity of the system is not too large). The most interesting case will be thus 2<k£3 (very significant system).

O THNX TO THE TI-92 ...

Enough theory. It is time that the TI-92 passes to the action! In the HOME screen , we start by defining the function f, and we evaluate for that the instruction: Define f(x)=x+k*x*(1-x).

Initially, we plot the curves y=f(x) which correspond to several values of K (of 0.5 in 0.5, from 0.5 to 3).

Here procedure: In MODE, select FUNCTION mode. In the «Y=» screen, create the function y1 by giving him the value y1=f(x). In HOME, evaluate: { 0.5,1,1.5,2,2.5,3 } -> k.

O SEQ. MODE

To study the sequence (u_n), we should choose the SEQUENCE mode. In «Y=» screen, one defines the sequence (u_n), and one gives the value of u₀ : u1=f(u1(n-1)), ui1=0.3 (value 0.3, selected here for u₀, is an example).

Let us recall that «nmin » is the index of the first term of the continuation (here u0); «nmax » is the index of the last studied term (choose it rather large); « plotstrt » is the sequence number of the first term to be traced (1 means that one will trace starting from u0). « plotstep » is the variation of index between two traced terms (1 is logical, because one will trace consecutive terms).

The other values define the window (traditional).

The choice of k will be done in HOME. Let us start with value 0.5: 5  ->  k

One enters on the GRAPH screen. The TI-92 traces the first bisectrix y=x, and the curve y=f(x). One actuates F3 to pass in TRACE mode. Supports on the cursor make it possible to " see " the behavior of the sequence.

Here, what one obtains with the data indicated previously :

Gradually, we increased the sensitivity of the system, which starts to react in a rather hard way. The first sign of disorder appeared with k >1: one sees that the system answers to a situation un<1 (lower population than ideal manpower) by un+1>1 (following population higher than this manpower) then by un+2<1, etc.

The tension appears clearly with k =2, for these « jolts » have sorrow to make converge the system towards its ideal manpower.

O INTEGRAL CHAOS

The best (2<K£3) remains to come ... We are indeed ready to study the behavior of the system for the values 2<K£3 (very significant system).

One could think that calm will return. To have the clear heart of it, one decides to set WEB/AUTO (F7 in the screen «Y=») and to pose (in screen WINDOW) «plotstrt=100 »(and always «nmax=200»). Thus the TI92 will trace only the behavior of the sequence, to u100 to u200.

Let us remember that the system always answers (at the time of the passage of un to un+1) by a displacement «in direction » of ideal manpower 1 (because the difference un+1-un is sign of 1-un). And this answer is all the more abrupt since k is large. If k>2, the fixed point 1 (ideal manpower) became « repulsive ». The sequence (un) live a conflict situation with the value 1, which at the same time attracts and pushes back it.

The maximum value k=3 shows that a complete disorder can result from this relation which one could describe as « sadomasochistic. »

O THE ORDER APPEARS

One can wonder if situations of balance do not appear all the same for less extreme values of K (i.e. if 2<k<3).

After one period of relative disorder, the system tends to oscillate between 4 distinct manpower (alternatively lower and higher than ideal manpower 1).

We now will try out other values (by keeping the same data as previously).

One could choose k to obtain cycles length 16, 32, 64, 128, ... ... ! But for other k, no cycle appears (as with k=3)...

Up to now, one varies k, and one saw that the behavior of the sequencen (with u0=0.3) was considerably modified. One can legitimately wonder if the value of u0 has the same influence.

If we take again the preceding values k=2.5, k=2.25, k=2.83, k=2.56, and if we test u0=1.1, u0=0.7, or u0=0.1, we obtain (after a certain time) the same cycles as those already observed.

To be honest, certain particular values of u0 can make exception (let us not speak about u0=0 or u0=1). Thus, with k=2.5, u0=1.2 produces a cycle length 2, on values 1.2 and 0.6: in this case the sequence (un) is quite simply 2-periodical ...

The cycles observed seem «necessary. » They are indeed because they are related to the fixed points of the fof, fofof, applications ... ..., but all that would involve us a little far!

We satisfy with a simple example: if k=2.25, the equalities f(a)=b and f(b)=a imply fof(a)=a and fof(b)=b. The reals a and b are thus two fixed points of the fof application.

O STUDY OF k= 3

In what precedes, it was possible to have a certain stability (cycles length 2,4,8, ... .) in an apparent disorder

The value k=3 seems much more desperate. No ordered evolution could be established. It is however there that us a surprise of size waits. There exists indeed, only in this case (it is a roof), a simple formula to calculate each un !

Let us suppose k=3. The sequence (un) is defined by the data of u0 and the recurrence relation u_n+1=u_n+3u_n(1-u_n)=4u_n-3u_n². One always supposed 0£u₀£(1+1/k)=4/3,which involves, for any n, 0£u_n£4/3. The double inequality 0£u₀£4/3 enables us to pose u₀=4/3*sin²(q), with 0£q£p/2. A small miracle is observed.

u₁ = 4u₀-3u₀² = 16/3*sin²(q)-16/3*sin⁴(q) = 16/3*sin²(q)cos²(q) = 4/3*sin²(2q). An obvious recurrence gives then, for all entire n: un = 4/3*sin²(2ⁿ*q), where q = Arcsin(Ö(3u₀/4)).

We did not finish with k=3, because this formula makes it possible to find values of the initial term u0 for which the sequence (un) is periodic!

Either p a positive integer (p>1). One knows that in the sequence (un), each term is defined only starting from the precedent. The sequence (un) is thus p-periodical, starting from the enteger m, as soon as u_m+p=u_m.

This equality is written: 4/3*sin²(2^m+p*q)=4/3*sin²(2^m*q), & for example : 2^m+p*q=2^m*q+qp (where q is an integer), i-e : q=qp/2^m/(2^p-1).

Let us summarize: with the initial value u₀=4/3*sin²(qp/2^m/(2^p-1)), (m³0, p³2, q³0) the sequence (un) is p-periodical starting from its um term. Here 2 numerical applications (in both cases the layouts are made starting from u0):

Where is the order, where is the disorder? Both seem so overlapping ... And it is not finished! because we will see that in this apparent tumble hide strange figures («the attractile ones). It's FRACTALS (that must say something to you...). Who will win : the order or chaos?

O THE RESEARCH OF THE CYCLES

Let us point out the definition of f : Define f(x)=x+k*x*(1-x)

k is a global variable here, being able to take any value in ]0,3]. One knows that if the initial term u0 of the sequence (un) is in ]0,1+1/k[, then it is the same for all the terms of the sequence (un). It was noticed that, if 0<k£2, the sequence converges towards 1 (in a monotonous way if 0<k£1, nonmonotonous if 1<k<2; quickly if k@1 and slowly if k@0 or k@2).

If 2<k<3, there are cycle-limits. For example, if k=2.25, the sequence (un) «tends »towards the cycle a®b®a, where a=0.715383 and b=1.173506 (f(a)=b and f(b)=a).

We propose to program the search for such cycle-limits, which one saw that they did not depend, except exception, of the initial value u0. We want to write a named function «cycles » taking in input the value to be given to the parameter k and returning the list of the values of the cycle obtained (empty list if no cycle is found).

In other words, it is necessary to find from when the sequence becomes periodic, and between which values it «turns ». The problem is that one does not know which will be the length of this cycle.

The solution is astute. The merit returns apparently to Martin Gardner in an article of Scientific American. The idea is to calculate in parallel x=un and y=u2n (the sequence (y) «advances » twice more quickly than the sequence (x)). If the sequence (un) becomes periodic as from a certain time, then one will necessarily have x=y.

It's then enough to compare for x=un the terms, u_n+1, u_n+2,..., until an equality un+p=un appears (necessarily p£n). We found the cycle, of lenght p : {u_n, u_n+1, u_n+2,..., u_n+p-1}

NOTE : it's only a cycle-limit. The equalities of the u_n+p=u_n type and the periodicity are true only on the computer (because of the precision limited to 12 digits).

The following listing is the translation of the preceding ideas. It will be noted that u0 is selected randomly between 0 and 1.

O The DIAGRAM OF FEIGENBAUM

It is seen that the size of the obtained cycles is not a simple function of K.

For certain values of K (k=2.569, for example), «cycles » does not find a period : either that this period does not exist, or that the value 500 (used in the listing) is too low to highlight it ...

The theoretical idea which comes to mind is to trace the the points (k,y), where, for each value of k, y is one of the points of the cycle which corresponding to K. For each value of K, in X-coordinate, one would be thus brought to trace 1,2,3,4,8... points according to the length of the cycle.

One thinks of using the «cycles» function , but the layout would take time.

One thus proceeds in a way a little «more artisanal. »For a given K, (and for selected value of u0), one starts by calculating the first 100 terms of the sequence, without trace. In this manner, one hopes to avoid the most disordered part (that which precedes the appearance of a cycle).

Then one traces the 100 following points (K in X-coordinate, y=un in ordinate).

It's possible that some points are confused (in particular if one entered during the cyclic period and if the period is very lower than 100).

For the values of k which would lead to very long cycles (128,256...) or with the absence of cycle, the display from these 100 points will give an idea of the disorder which still reigns at this stage of the sequence.

One will carry out layouts over all the width of the screen, therefore for 239 values of K (TI-92/92 Plus). For each k, one calculates 100 values of un, then one traces 100 others of them. Thus, one calculates almost 50000 values. The consequence is immediate: it is long ! You can accelerate the layout (by decreasing value 100 in the listing), but that will be never very precise.

Here, the listing of the program Feigenbm().

Here, what one obtains, with the window: xmin=1.8; xmax=3; ymin=0, ymax=1.4 (the calculating time is 45 minutes! with a TI-92v1.13).

The obtained figure is known under the name of diagram of Feigenbaum (American physicist whose work dates from the Seventies).

Its structure is attractive: one sees many points of junction (when the length of the cycle doubles). Sometimes one passes from an apparent chaos to a calmer situation. As any fractal figure, it has autosimilarity properties (one can zoom and recognize the original figure). We will finish this paragraph by some of these zooms.

The resemblance to the original is striking. The accuracy error towards the first junction is due o the program: it happens that the computation of the first 100 terms (without tracing them), is not enough to enter to the cyclic behavior.

The Feigenbm() program starts to show its limits. Don't forget that the current window represents approximately 2/10000 of the initial window, and that the majority of the "lit" points fall apart from the new window. The inaccuracies (first junction) can be corrected by increasing value 100 in Feigenbm(), therefore by increasing the calculating time ...

It is shown that the values of k which produce junctions (doubly cycles) are approximately: k1=2; k2=2.44949; k3=2.54409; k4=2.56441; k5=2.56876; k6=2.56969; k7= 2.56989...

These values define intervals: d1 = k2-k1, d2 = k3-k2; d3 = k4-k3, ... ...

Finally, the ratios of these lengths: d1/d2, d2/d3, d3/d4, d4/d5, ... ... , tend towards D @ 4.66920, called constant of Feigenbaum (discovery in 1975).

This mysterious constant is universal. It very frequently appears in the fractals, and one also meets it in many physical phenomena (dynamics of the fluids, physics of the lasers, acoustics, electronics...). In fact the constant of Feigenbaum could have the same importance in nature as constant p ...

O THE «BUTTERFLY» EFFECT

This title is derived from an article of meteorologist (E.N.Lorenz), who, in the Sixties, studied how much the forecasting models weather were sensitive to a negligible variation of the initial conditions. It then wondered whether a tornado at the top of Texas could not be caused by the flight of one butterfly to the other end of the sphere...

We will study the butterfly effect on the behavior of the sequences (un) defined by u_n+1=u_n+k*u_n*(1-u_n), with k=3 (thus for the maximum sensitivity). The initial term u0 is in #[ 0 , 1+1/k ] (and it is then the same of all the un).

In the MODE, we choose SEQUENCE. In the «Y= »we pose: u1=u1(n-1)+3*u1(n-1)*(1-u1(n-1)).

We saw previously, that the initial value u₀=4/3*sin²(2p/ 7) led to a sequence (u_n) periodic of period 3. It is checked that: u₁=4/3*sin²(4p/7), u₂=4/3*sin²(8p/7), and u₃=4/3*sin²(16p/7)=u₀. And more generally u_n+3=u_n.

In «Y=», one poses ui1=4/3*sin²(2p/7) (let us recall that it is about the initial term of the first sequence. For us it is u0 ).

One must thus expect to see «rolling »the sequence (un) on the three successive values u0 » 0.815014, u1 »1.26731, u2 »0.251007.

On GRAPH, one passes in mode TRACE, to visualize the behavior of the sequence.

What did it occur? The initial value that we believed to give to u0, ( u₀=4/3*sin²(2p/7) ), was replaced by an approximate value by using the internal precision of the calculator.

This initial error, negligible, did not cease swelling until the behavior of the sequence was completely changed by it (the error between truth un and the TI-92' un have the size of un itself).

The initial butterfly was thus transformed into an enormous ventilator. The TI-92 is not in question. Any numerical system coding reals with a finished number of decimals would undergo this phenomenon.

To describe what we observed, i.e. an extreme dependence compared the initial conditions, in a process which however seems nothing to leave randomly (each term is calculated according to the precedent, by a concrete formula ), the Anglo-Saxons speak about «Deterministic chaos» .

O OTHER EXPERIMENTS

Here another solution of checking the exacerbated sensitivity of the system which we are studying.

In small F6, we choose Square for the two sequences. In F7, we choose CUSTOM for the axes, with u1 (first sequence) in X-coordinate and u2 (second sequence) in ordinate. The TI-92 will build the points of co-ordinates (u1(N),u2(N)). Theoretically these points should practically be on the first bisectrix...

A last experiment should complete to convince you (or to worry you). Let us note that the formula of recurrence u_n+1=u_n+3*u_n(1-u_n) can be written, more simply: u_n+1=4u_n-3u_n².

The same disorder appeared! You will not undoubtedly see your calculator of the same way now... Once again, no system of numerical calculation can resist the butterfly «effect». A computer which would work with 100 decimals would fall into the same trap, but more slowly... The phenomenon which has been just highlighted produces however only on quite particular occasions.

One saw on several occasions how much the value of the coefficient k could have importance in the behavior of the sequences u_n+1=u_n+ku_n*(1-u_n). The low values of K (for ex. 0<k<2) lead to stable «systems».Variations on the initial term do not prevent the sequence to converge towards 1.

The following example offers a striking illustration of it. One chose k=0.05 (for a slow convergence) and two sequencens defined by the usual relation, the initial value being 0.1 for a sequence and 1.4 for the other. One chose xmin=0, xmax=238, ymin=0, ymax=1.4 and mode TIME for the axes.

O LET US OBSERVE THE BUTTERFLIES

Until now we made only observations. We now propose to study the role of the coefficient K in a little more precise way.

The studied sequences check the relation: u_n+1=f(u_n) avec f(x)=x+k*x*(1-x). Let us recall that all the terms of the sequence are in the interval [ 0 , 1+1/k ]. What interests us, it is of knowing how an initial error e0 on the u0 term can be propagated.

Suppose that "en" is the absolute error on the term un. By using a traditional result of analysis: e_n+1 = |f(u_n+e_n)-f(u_n)| » e_n*|f '(u_n)|.

However F ' (X) is equal to: 1+k-2k*x. Thus e_n+1/e_n » |1+k-2k*u_n|. There is an estimation of the factor by which the error is multiplied, in the passage of row N to the row n+1. In the passage of row 0 to row N, the initial error is multiplied by the factor: (e₁ /e₀)(e₂ /e₁)....(e_n /e_n-1)=(e_n/e₀).

If the error underwent a rate of increase constant l, this rate should check: lⁿ=(e₁ /e₀)(e₂ /e₁)...(e_n /e_n-1) or : ln(l)=1/n[ln(e₁ /e₀)+...+ln(e_n/e_n-1)]

And finally, with the notations of the TI-92:

ln(l)=1/n*S(ln(e_m+1 /e_m),m,0,n-1)=1/n*S(ln(|1+k-2k*u_m|),m,0,n-1).

For a given value of u0, one can calculate the sum of the second member, and take exponential result. One thus obtains the average coefficient of amplification of the error on N first rows of the sequence.

Here a program which calculates the average rate of error amplification, of u0 to un, for each value of N.

This average rate is displayed with the value of N which corresponds to him. Syntax is lambda(u₀,k), where u0 is the initial term and K the coefficient.

If the error is initially 1E-12, then at the end of 40 iterations, it is about (240)*(10-12)»1: the error covers completely the signal which one must observe... One will thus not be surprised, in one of the preceding examples, that chaos clearly «appeared » during the calculation of u41.

NB: preceding calculations are correct if the total error remains small.

O CONCLUSION

We have just seen how, with a minimum of programming, and especially thanks to its SEQUENCE environment the TI-92 enabled us to tackle this impassioning subject of the discrete dynamic systems.

One does not know who won ... the order or the chaos (and so finally it was the same thing?).

Let us recognize how much the TI-92 was invaluable to us to undertake this study.

JMF