Difference Equations and Chaos in Mathematica
by Robert Knapp and mark Sofroniou

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

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

(c)
In[3]:= Last[iterates]
Out[3]= 0.6

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

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

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


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]];

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


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

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}}

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

Out[9]= {{0.5, 0.5}, {0.6, 0.580782}, {0.6, 0.580782}}

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

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

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

Listing Four
(a)
In[13]:= Exp[10000.]

Out[13]= 8.806818225663 104342

(b)
In[14]:= Compile[{{x,_Real}},Exp[x]][10000.]

CompiledFunction::cfn: Numerical error encountered at instruction 3; proceeding
with uncompiled evaluation.

Out[14]= 8.806818225663 104342


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"];]


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]];


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

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

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

    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;
    }

    MLPutRealList(stdlink, result, combine);

    free(result);
}

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

    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.;
    }

    return lexp/iterates;
}


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

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

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

    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;
    }

    MLPutRealList(stdlink, result, combine);

    free(result);
}

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

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

    return lexp/iterates;
}


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:

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


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];


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

(b)
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]}}]];


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

Out[21]= {{4,0.693142}}

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

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

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

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

(d)
In[24]:= Nest[LogisticFunction[4], N[Sin[1]2], 53]

Out[24]= 0.279937

(e)
In[25]:= exact = Sin[253]2;

(f)
In[26]:= N[exact,50]

Out[26]= 0.72067529373649285832360770596944319298599918355010

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

(h)
In[28]:= Nest[LogisticFunction[4], N[Sin[1]2, 50], 53]

Out[28]= 0.720675293736492858

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

Out[29]= 0.720675293736492858323607705969443







DDJ
