Difference Equations and Chaos in Mathematica

Symbolic and numerical mathematics at work

Rob is director of numeric technology for Wolfram Research, and Mark is a kernel developer specializing in numeric and algorithmic development. They can be contacted at [email protected] and [email protected], respectively.

A difference equation (or map) of the form x_n+1=f(x_n,x_n-1,...) which, together with some specified values or initial conditions, defines a sequence {x_n}. Despite the seemingly simple form, difference equations have a variety of applications and can display a range of dynamics. Since maps describe iterative processes, they come up frequently in computer science. Also, many of the approximations in numerical analysis (such as numerical solutions of differential equations) typically approximate continuous dynamical systems using discrete systems of difference equations. Modeling a map using a computer is equivalent to studying the process of functional composition or functional iteration.

An interesting set of discoveries was made for a restricted class of maps of the form x_n+1=f(x_n) with 0<=f(x_n)<=1 when 0<=x_n<=1. In "Simple Mathematical Models with Very Complicated Dynamics" (Nature, No. 26, 1976), R.M. May used the logistic map, where f(x)=rx(1-x). A Mathematica definition is given in Listing One:(a) Notice how LogisticFunction can be defined using currying; for example, f[x,y] can be written as f[x][y]. Some functional languages only support functions of one argument, and this is the only way that multivariate functions can be implemented (see The Mathematica Programmer II, by R.E. Maeder, Academic Press, 1996).

The logistic map models populations with discrete generations (such as grasshoppers) and is normalized to carrying capacity. The parameter r is determined by the base birth rate. When the population becomes large, the net birth rate decreases in a nonlinear way because of lack of resources. When 0<=r<=4, the logistic function satisfies the conditions above. May discovered that as the parameter was increased from 0 to 4, there were different regimes of behavior, divided by so-called bifurcations. For a sufficiently small r, the population would die out. Next, there was a set of values of r where the population would approach a stable equilibrium number. Above that, he found population cycles, where the population would vary in alternate generations, but as the parameter was increased, the length of the cycles doubled, from 2 generations, to 4, to 8, and so on. Beyond this, the population showed seemingly random behavior -- the population would vary from generation to generation in an unpredictable chaotic way. In fact, this provides a classic example of a transition from regular behavior to deterministic chaos, which is characterized by unpredictable behavior with extreme sensitivity to initial values.

Mathematica provides an environment that makes it possible to explore problems such as difference maps. It includes tools for symbolic computation, numerical computation, and graphics. These tools use the same fundamental data structure so the results from one command can be shared with others in a natural way. The glue that unifies the system is the Mathematica programming language.

In this article, we will set up programs in Mathematica and use links to external programs to demonstrate dynamical properties. Using the interactive notebook version (available electronically at http:// www.wolfram.com/special/iteration/), you can investigate them further on your own.

Simple Iterations

The Mathematica NestList command repeatedly composes (or iterates) a function, given an initial condition and returns the values in a list. Using the parametric function LogisticFunction with r=2.5, you can form the sequence of iterates in Listing One:(b) and visualize the results as in Figure 1. You can see that the population approaches an equilibrium value. You can find (an approximation to) the equilibrium value by looking at the last iterate; see Listing One:(c).

More interesting is that with r=3.1, there is a period-two cycle, as Listing One:(d) and Figure 2 illustrate. The cycle is stable -- no matter what initial value you start with, eventually the population settles down to a steady state and cycles between the two values in Figure 2. Of greater interest in the scientific community at large, however, was that as the parameter r was increased, the period of the cycle doubled. Eventually, for large enough r, the period has in effect doubled infinitely many times. Just above this value, there are cycles of every period 2^k, k=0,1,2,... -- but none are stable. In fact, the population varies from year to year in an unpredictable, seemingly random way, known as "deterministic chaos." For example, with r=3.6, Listing One:(e) results in Figure 3.

The population never settles down to a consistent cycle and is said to be "chaotic." One of the key components of chaos is sensitivity to initial conditions. For instance, if you start with a value near Listing One:(e), the distance between iterates eventually becomes large, as Listing One:(f) and Figure 4 illustrate.

The exponential rate of divergence is reflected by the Lyapunov exponent of the system. We will investigate the sensitivity to initial conditions in more detail and show how the Lyapunov exponent of the system can be computed.

The Bifurcation Diagram

As the parameter r varies, the point at which the period of iterates doubles is known as a "bifurcation point." To illustrate, we present programs that let you investigate the sequence of period-doubling bifurcations and the transition to chaos. A bifurcation diagram lets you visualize what values the population (or iterates of a map) take on as a function of the parameter. The idea works like this: You start with some initial value and iterate the map many times. Then iterate the map several more times, and remove duplicate values among those iterates. If the parameter is in a range where there is a stable periodic cycle, then the number of values you have will be the period of the cycle. On the other hand, if you are in a range where there is chaos, then all of the values will be different, and you will see a range of points. To successfully make a bifurcation diagram, you need to iterate the map many times for many different values of the parameter.

Mathematica's programming language lets you set up a general framework for solving a class of problems. Instead of writing a function that works for a particular map (such as the logistic map), you can design your program to accept an arbitrary specification for the map and its variables. This is possible because Mathematica lets you use higher-order functions to create values of functional type at run time.

A common procedural way of building up the points to plot in the bifurcation diagram is to iterate the map separately for each value of r. However, because Mathematica's arithmetic is distributive over lists, we can do the iterations for a set of values of r at once. Listing Two generates the points to make a bifurcation diagram. Notice how the combination of high-level functions makes the program short. For example, Union removes duplicates from a list of values, taking into account that rounding errors require that floating-point numbers should be compared using a tolerance to allow for discrepancy of some of the least significant bits. This lets users focus on conceptual details, not programming details. The use of functional programming constructs makes the program general and extensible and virtually loop free. (For reasons that will become clearer, we construct the function using a separate routine called MakeMapFunction.) Trying the BifurcationDiagram program for a small number of values of r lets you see the form of the output in Listing Three:(a) and also that between r=3.3 and r=3.4, there has been a bifurcation from period 2 to period 4. We can specify a different map, such as the sine map. Listing Three:(b) precomputes a numerical approximation to the constant using N to avoid doing the coercion at each iteration. This map also undergoes at least one bifurcation.

You can see from the output of the logistic map that, for a few values, there is a bifurcation between 3.4 and 3.5. We are ready to generate the data which can be used to make a bifurcation diagram. The BifurcationDiagram program generates a decent diagram, though it may take some time. Since Mathematica is a typeless or polymorphic interpreted language, it has to determine types at run time. Hence, it will not be as fast as compiled code. For simple calculations like this, the overhead is noticeable. For larger problems (such as in matrix algebra), most of the work is done in intrinsic functions and not by the interpreter, so overhead is minor. There are two relatively simple ways to speed up the calculation. First, you can use the Mathematica compiler; second, you can set up an external program that will do the iterations and send the result back to Mathematica for subsequent plotting or analysis.

Most of the work in the BifurcationDiagram program is spent in the evaluation of the function defining the map. By compiling the function definition we can speed up evaluation considerably. Mathematica's numerical compiler compiles Mathematica commands into pseudocode instructions, which are interpreted using machine numbers. The result is a CompiledFunction object, which contains the sequence of virtual-machine instructions describing how the compiled function should execute. Unlike binary machine-code instructions, the virtual-machine description is portable, allowing compiled programs to be transferred between different versions of Mathematica running on different platforms. The compiler supports various types of numbers as well as tensors. The rank of the tensor is specified together with the type at compile time, where a rank 1 tensor is a vector or list. By overwriting the definition of our map constructor function, we can create a compiled map operating on a range of parameters and values; see Listing Three:(c). The Compile command takes type information. In this case, r and x are specified as being rank 1 tensors (vectors) of machine floating-point real numbers.

In production-quality code, you could add an option to BifurcationDiagram to let users select between the compiled and uncompiled versions. Many of Mathematica's numerical functions have the option Compiled to do this; see Listing 3(d).

The compiled version gives the same results as the uncompiled version, but takes a fraction (roughly 1/6) of the time to produce data for a bifurcation diagram, letting you visualize how the period-doubling bifurcations cascade into chaos; see Listing Three:(e) and Figure 5. Using different settings, you can zoom into different portions of the plot to see the fine scale of the graph. If you do lots of zooming in, Mathematica compiled code cannot generally match the speed of code written in a compiled language such as C, particularly for simple iteration examples like this. Why? First, even though the pseudocode produced by the Mathematica compiler is relatively efficient, it does not produce binary code at the machine-instruction level. Second, although it is generally hidden from users, Mathematica is thorough in its machine exception/overflow checking. For example, with many programs, Listing Four:(a) returns an IEEE floating-point infinity, rather than an actual value.

Mathematica reverts to arbitrary-precision software arithmetic, which has a larger exponent range. If a run-time error occurs during the execution of the virtual machine instructions in a CompiledFunction object, the uncompiled Mathematica function continues the evaluation; see Listing Four:(b). This behavior is important because it lets user programs continue to execute even if a compiled step is deep inside a complicated program.

Linking to C

For many programs, it is possible to use C or other compiled code for the time-critical part of the calculation. Mathematica's MathLink protocol for communication with external programs makes this possible. The protocol is general, and MathLink programs have been written that let you access Mathematica from within, for example, Microsoft Word or Excel. MathLink lets you transmit data based on the tree data structure, which is uniform in Mathematica. You can establish communication between processes on one or several computers. (In fact, Mathematica uses MathLink to communicate between its own kernel and front end.)

Listing Five shows how to link Mathematica to a C program that carries out the iterations. The time-critical portion of the calculation is the actual iteration. Also, since C does not do operations for an entire list at a time, the C program will do the computation of a single value of r at a time. To generate the program for a general functional map, the Mathematica Splice is used to set up a code-generation scheme. In Listing Five, the Mathematica MathLinkSpliceAndCompile function generates the C code for a given map, compiles it, sets up a communication link to the compiled executable, and calls it to do the desired iterations. MathLinkBifurcationDiagram uses Install to set up communication between Mathematica and the C program, then call a Mathematica function, which in turn calls the C map_iterate function; see Listing Six. Listing Seven is the map.mc file.

When Splice["map.mc"] is executed, the items inside the delimiters <* and *> are replaced with Mathematica evaluations in a new file, map.c. In Mathematica, a Block is a programming construct where variables keep their names, but are assigned local values. Since Splice is called from within the Block, it has a defined value and it evaluates to the user's function in terms of the variables r and x. It will work as long as simple functions that have standard C library implementations are called. If, for example, you need access to special functions, it is not difficult to design the code to call back to Mathematica to get those values. Also, you can use more sophisticated code-generation capabilities by doing functional transformations on your expressions in Mathematica, such as the Horner or nested form of a polynomial. You can also use common-subexpression optimization to save some numerical work. (For more information, see Mark Sofroniou's "Extending the Built-in Format Rules" in The Mathematica Journal, 3 3, 1993, and "Expression Optimization in Mathematica," in Elements of Mathematica Programming, Troels Petersen ed., TELOS/Springer Verlag, 1998.) Listing Eight is the C file generated by MathLinkBifurcationDiagram[r x (1- x), {r, 3, 4,.0005}, {x,0.1},10000,100].

The next stage is to specify how the C function is accessed in Mathematica. There is a template mechanism provided with MathLink that makes this straightforward. Listing Nine (map.tm) specifies the correspondence between Mathematica and C functions: It defines MLMapIterate and MLLyapunovExponent to be the Mathematica functions corresponding to map_iterate and lyapunov_exponent, respectively. Unlike in Mathematica, which uses reference counting to reclaim memory, we have to focus on programming details like memory allocation and deallocation. The MathLink function (including the time for compilation) runs faster than the function produced by the Mathematica compiler, and gives an essentially identical plot. Consequently, you can consider interactively zooming in on particular ranges of r to investigate details of the bifurcations or other features of the plot.

The Lyapunov Exponent

The Lyapunov exponent measures the exponential rate of divergence (or convergence) of nearby initial conditions. Figure 6 shows how it is calculated for a given sequence of iterates, x_i. Numerically, l can be approximated using a finite number of iterations. The Mathematica MathLinkLyapunovExponent function uses the lyapunov_exponent function in the generated C program to compute numerical values; see Listing Ten.

Now we can make a plot of the Lyapunov exponent as a function of the parameter using Listing Eleven:(a) and Figure 7. Where the Lyapunov exponent is negative, there is convergence to a fixed point or cycle. Where it is positive, nearby initial conditions diverge. You can see this by comparing with the bifurcation diagram. Notice that even beyond the transition to chaos, there are values of r for which the Lyapunov exponent is negative. These correspond to "windows" where cycles of different prime number periods appear and undergo their own period-doubling bifurcations. An example is where a period-three cycle appears and bifurcates near r=3.85. Since the Lyapunov exponent as a function of r has structure on any scale, the resolution of the plot is insufficient to capture the entire structure of the function. However, the Mathematica function makes it trivial to zoom in on any portion of the range. Listing Eleven:(b) and Figure 8 show a part of the bifurcation diagram and the Lyapunov exponent for the values of r right near the onset of chaos at approximately 3.5699.

Loss of Precision

One consequence of a positive Lyapunov exponent is numerical loss of precision. When there is exponential divergence of nearby initial conditions, even small roundoff errors will eventually be magnified. Let's take a closer look at what is happening at the value r=4. At this value, our approximations show that the Lyapunov exponent is as defined in Listing Twelve:(a), which is very near to Log[2]. Actually, as Listing Twelve:(b) shows for r=4, there is a transformation that makes it easier to understand the divergence. With one more step, it becomes Listing Twelve:(c). The iterations are simple. The message from Solve indicates that since the Sin function is not one-to-one, only the principle value has been found (q_n+1=2(q_n + p)) is also a solution for example). If the q_n are thought of as binary numbers, the transformed map is just a shift operator. Since the distance between two nearby points is doubling at each iteration, this shows that the Lyapunov exponent should be exactly ln(2). Another way to look at this is that small errors will be doubled at each iteration: In other words, one bit of precision is lost per iteration. Consequently, if you are doing the calculations with IEEE double-precision numbers, which have 53 bits of mantissa precision, you can expect to have lost all precision by 53 iterations. To verify this, start with the initial value sin² (1) approximated by the machine floating-point number in Listing Twelve:(d). On the other hand, the exact value of the 53rd iterate should be Listing Twelve:(e).

Evaluation using machine numbers is inaccurate, since the C math library function sin() (which Mathematica uses for machine numbers) needs to do argument reduction before it approximates the sine. Argument reduction on a number of this size will result in severe loss of precision (in this case in a single step). Mathematica's adaptive-precision evaluation gets a result that is correct to the requested precision; see Listing Twelve:(f).

N works by adaptively raising the precision to which subexpressions are approximated until a result is found to the requested precision. As you can see in Listing Twelve:(g) and Figure 9, none of the digits of the 53rd iterate computed with machine numbers are correct. Figure9 shows successive loss of bits. One solution is to start with higher precision. For instance, Listing Twelve:(h) obtains a result with 18 digits assured to be correct. You may have noticed that the loss of precision exceeded 53 bits (32 digits~106 bits). The reason for this is that significance arithmetic works by doing a worst-case analysis assuming all errors are independent. In this case, the combination of operations is not as bad as it could possibly be; when you know how much precision is lost per iteration, you can override the automatic precision control. Listing Twelve:(i), for instance, only loses 53 bits.

Since all these plots were calculated using machine numbers (in many cases with thousands of iterations with a positive Lyapunov exponent), you may be wondering if they have any value at all. Though not strictly accurate, the graphs do show the overall (or ergodic) properties of the map as they are intended to. If you needed to have a particular value of some iterate, you would need to use higher precision. However, this could take substantial time if it was thousands of iterates out, since you would need to start with thousands of bits of precision. Nonetheless, the framework of arbitrary-precision arithmetic is there in Mathematica when you need it.

Conclusion

Even though Mathematica has symbolic algorithms for solving difference equations, no general-solution method exists when the equations are nonlinear. Therefore, we often have no option but to resort to numerical methods and experimentation to investigate the dynamics. The various arithmetic models available in Mathematica enable a wide study of phenomena that would otherwise fall outside the scope of conventional programming languages. For example, there has been a recent resurgence of interest in the use of interval arithmetic methods to obtain validated computer-based proofs of certain physical properties.

For More Information

Wolfram Research
100 Trade Center Drive
Champaign, IL 61820
217-398-0700
http://www.wolfram.com/

Listing One

In[1]:= LogisticFunction[r_][x_] := r x (1 - x);

</p>
(b)
<pre>In[2]:= iterates = NestList[LogisticFunction[2.5], 0.1, 30];        ListPlot[iterates, PlotJoined->True,PlotRange->All];


</p>
(c)
<pre>In[3]:= Last[iterates]Out[3]= 0.6


</p>
(d)
<pre>In[4]:= iterates = NestList[LogisticFunction[3.1], 0.1, 30];        ListPlot[iterates, PlotJoined->True];


</p>
(e)
<pre>In[5]:= iterates= NestList[LogisticFunction[3.6], 0.1, 50];        ListPlot[iterates, PlotJoined->True];


</p>
(f)
<pre>In[6]:= iterates1 = NestList[LogisticFunction[3.6], 0.1001, 50];        ListPlot[iterates- iterates1];

Back to Article

Listing Two

In[7]:= BifurcationDiagram[f_,{r_, rmin_, rmax_, rstep_}, {x_, x0_}, start_,     combine_] :=
        Module[{R,temp, MapFunction},           
            R = Table[r,{r,rmin,rmax,rstep}]; 
              (* The range of values of the parameter *)                
            MapFunction = MakeMapFunction[{r,x},f];
              (* Construct the iterating function *)            
            temp = Nest[MapFunction[R, #]&, x0 + 0.R,start+1];
              (* Starting iterates *)           
            temp = NestList[MapFunction[R, #]&,temp,combine-1];
              (* Subsequent iterates *)         
            temp = Map[ Union, Transpose[temp] ]; 
              (* Remove duplicate values from cycles *)         
            Flatten[ MapThread[Thread[{#1,#2}]&,{R,temp}],1]];


</p>
        MakeMapFunction[{r_,x_},f_]:= Function[{r,x},f];

Back to Article

Listing Three

In[8]:= BifurcationDiagram[r x (1 -x), {r, 3.1, 3.5,.1},{x,.1},10000,100]

</p>
Out[8]= {{3.1, 0.558014}, {3.1, 0.764567}, {3.2, 0.513045}, {3.2, 0.799455},
        {3.3, 0.479427}, {3.3, 0.823603}, {3.4, 0.451963}, {3.4, 0.451963}, 
        {3.4, 0.842154}, {3.4, 0.842154}, {3.5, 0.38282}, {3.5, 0.500884}, 
        {3.5, 0.826941}, {3.5, 0.874997}}


</p>
(b)
<pre>In[9]:= BifurcationDiagram[a Sin[N[p]x],{a,.5,.6,.1},{x,.1},10000,  100]


</p>
Out[9]= {{0.5, 0.5}, {0.6, 0.580782}, {0.6, 0.580782}}


</p>
(c)
<pre>In[10]:= (* Construct a compiled version of the iterating function *) MakeMapFunction[r_,x_,f_]:=  Compile[{{r, _Real, 1},{x, _Real, 1}},f];


</p>
(d)
<pre>In[11]:= Options[Plot,Compiled]Out[11]= {Compiled -> True}


</p>
(e)
<pre>In[12]:= ListPlot[BifurcationDiagram[r x (1 - x),{r,3,4,0.001},{x,.1}, 10000,100],  
                    PlotStyle->AbsolutePointSize[0.0001]];

Back to Article

Listing Four

In[13]:= Exp[10000.]

</p>
Out[13]= 8.806818225663 104342


</p>
(b)
<pre>In[14]:= Compile[{{x,_Real}},Exp[x]][10000.]

</p>
CompiledFunction::cfn: Numerical error encountered at instruction 3; proceeding
with uncompiled evaluation.


</p>
Out[14]= 8.806818225663 104342

Back to Article

Listing Five

In[15]:= MathLinkSpliceAndCompile[f_,  rvar_, xvar_] := (* Block is a programming construct which gives local values 
           to the symbols given in the first argument.  This is 
           different than Module, which creates unique local symbols 
           to avoid any possibility of conflict.  Block is needed here
           because x and r need to be the same symbols which the user 
           specified as the variables *)                
        Block[{MapFunction, r,x},                       
            MapFunction[r_][x_] = (f /. {rvar->r,xvar->x});             
              (* The Splice command creates a new file, with 
                 <* MapFunction[r][x] *> replaced with its current 
                 (local) value in the Mathematica session *)            
            Splice["map.mc"];                   
              (* These two commands acess the (Linux) operating system 
                  to remove the old executable and compile the code *)  
            Run["rm -rf map"];                  
            Run["mcc -o map map.c map.tm >& mcc.out"];]

Back to Article

Listing Six

In[16]:= MathLinkBifurcationDiagram[f_, {rvar_, rmin_, rmax_, rstep_}, {xvar_, x0_}, start_, save_] :=   
        Module[{link, BifurcationList},         
            MathLinkSpliceAndCompile[f,rvar,xvar];                      
              (* Install starts the executable (map) and sets up 
                 a communication link between it and the                
                 Mathematica session.  If for some reason--say,
                 the code didn't compile properly--the connnection 
                 cannot be established, Install returns $Failed *)      
            link = Install["map"];                      
            If[link === $Failed,Return[$Failed]];               
              (* The template file map.tm below associates the C 
                  function map_iterate with the Mathematica             
                  function MLMapIterate *)              
            BifurcationList = Table[
                {r, #}& /@  Union[MLMapIterate[r, x0, start, save]],    
                {r,rmin,rmax, rstep}];          
            Uninstall[link];
            Flatten[BifurcationList,1]];

Back to Article

Listing Seven

#include "mathlink.h"#include "mdefs.h"
#include <stdlib.h>
#include <math.h>


</p>
int main(argc, argv)
{
    int argc; char* argv[];
    return MLMain(argc, argv);
}


</p>
void map_iterate(double r, double x, int start, int combine)
{
    int i;
    double *result;


</p>
    for (i = 0; i < start; i++) 
        x = <* MapFunction[r][x] *>;
    
    result = (double *) malloc(start*sizeof(double));
    for (i = 0; i < combine; i++) 
    {
        x = <* MapFunction[r][x] *>;
        result[i] = x;
    }


</p>
    MLPutRealList(stdlink, result, combine);


</p>
    free(result);
}


</p>
double lyapunov_exponent(double r, double x, int iterates)
{
    int i;
    double lexp = 0.;


</p>
    for (i = 0; i < iterates; i++) 
    {
        x= <* MapFunction[r][x] *>;
        lexp += log(fabs(<* Collect[D[MapFunction[r][x],x],r] *>));
        if (!finite(lexp)) return 0.;
    }


</p>
    return lexp/iterates;
}

Back to Article

Listing Eight

#include "mathlink.h"#include "mdefs.h"
#include <stdlib.h>
#include <math.h>


</p>
int main(argc, argv)
{
    int argc; char* argv[];
    return MLMain(argc, argv);
}


</p>
void map_iterate(double r, double x, int start, int combine)
{
    int i;
    double *result;


</p>
    for (i = 0; i < start; i++) 
        x = r*(1 - x)*x;
    
    result = (double *) malloc(start*sizeof(double));
    for (i = 0; i < combine; i++) 
    {
        x = r*(1 - x)*x;
        result[i] = x;
    }


</p>
    MLPutRealList(stdlink, result, combine);


</p>
    free(result);
}


</p>
double lyapunov_exponent(double r, double x, int iterates)
{
    int i;
    double lexp = 0.;


</p>
    for (i = 0; i < iterates; i++) 
    {
        x= r*(1 - x)*x;
        lexp += log(fabs(r*(1-2*x)))
        if (!finite(lexp)) return 0.;
    }


</p>
    return lexp/iterates;
}

Back to Article

Listing Nine

:Begin::Function:      map_iterate
:Pattern:       MLMapIterate[r_, x_, start_Integer, combine_Integer]
:Arguments:     { r, x, start, combine }
:ArgumentTypes: { Real, Real, Integer, Integer }
:ReturnType:    Manual
:End:


</p>
:Begin:
:Function:      lyapunov_exponent
:Pattern:       MLLyapunovExponent[r_, x_, iterates_Integer]
:Arguments:     { r, x, iterates }
:ArgumentTypes: { Real, Real, Integer }
:ReturnType:    Real
:End:

Back to Article

Listing Ten

In[18]:= MathLinkLyapunovExponent[f_,  {rvar_, rmin_, rmax_, rstep_}, {xvar_, x0_}, iterates_] :=         
        Module[{link, LEList},          
            MathLinkSpliceAndCompile[f,rvar,xvar];              
            link = Install["map"];              
            If[link === $Failed,Return[$Failed]];               
              (* The template file map.tm below associates the 
                 C function map_iterate with the Mathematica 
                 function MLMapIterate *)               
            LEList =Table[{r,MLLyapunovExponent[r,x0,iterates]},
                          {r,rmin,rmax,rstep}];
            Uninstall[link];            
            LEList];

Back to Article

Listing Eleven

In[19]:= ListPlot[MathLinkLyapunovExponent[                r x (1 - x),{r,3,4,0.001}, {x, .1},10000],  
            PlotJoined->True];


</p>
(b)
<pre>In[20]:= Show[GraphicsArray[{                  {ListPlot[MathLinkBifurcationDiagram[
                    r x (1 - x),{r,3.5698,3.57,.000001},{x, .1},100000,1000],
                 PlotStyle->AbsolutePointSize[0.0001],
                 PlotRange->{0.8045,.805}]},        
            {ListPlot[MathLinkLyapunovExponent[            
                    r x (1 - x),{r,3.5698,3.57,.000001},{x, .1},100000],
                  PlotJoined->True]}}]];

Back to Article

Listing Twelve

In[21]:= MathLinkLyapunovExponent[r x (1 - x), {r,4,4,1},{x,.1},100000]

</p>
Out[21]= {{4,0.693142}}


</p>
(b)
<pre>In[22]:= Simplify[xn + 1 == 4*xn*(1 - xn)/. xn_ :> Sin[qn]2]

</p>
Out[22]= Sin[qn+1]2  == Sin[qn]2
(c)
<pre>In[23]:= PowerExpand[Solve[%, qn+1]

</p>
Solve::ifun: Inverse functions are being used by Solve, so some solutions may
not be found.


</p>
Out[23]= {{q1+n -> 2*qn}}


</p>
(d)
<pre>In[24]:= Nest[LogisticFunction[4], N[Sin[1]2], 53]

</p>
Out[24]= 0.279937


</p>
(e)
<pre>In[25]:= exact = Sin[253]2;

</p>
(f)
<pre>In[26]:= N[exact,50]

</p>
Out[26]= 0.72067529373649285832360770596944319298599918355010


</p>
(g)
<pre>In[27]:= ListPlot[Log[2, Abs[NestList[LogisticFunction[4], N[Sin[1]2], 53] -      Table[N[Sin[2j]2, 20], {j, 0, 53}]]]];


</p>
(h)
<pre>In[28]:= Nest[LogisticFunction[4], N[Sin[1]2, 50], 53]

</p>
Out[28]= 0.720675293736492858


</p>
(i)
<pre>In[29]:= loss = Log[10., 2];          prec = 50 + loss;   
         Nest[(prec -= loss; LogisticFunction[4][SetPrecision[#1, prec]])&, 
              N[Sin[1]2, prec], 53]


</p>
Out[29]= 0.720675293736492858323607705969443

Back to Article

DDJ