Legal Notice, Disclaimer, and Book Ordering Information:

(* Legal Notice, Disclaimer, and Book Ordering Information:

Last Modified: FIO June 8, 2002

Mathematica For Physics: 2nd Edition

Robert L. Zimmerman and Fredrick I. Olness

WebSite:  
  www.physics.smu.edu/~olness
  darkwing.uoregon.edu/~phys600/

MathSource Number: 0206-862

ISBN 0-201-53796-6
For ordering information, call 1-800-282-0693
Addison-Wesley Publishing Company

Communication with the authors:
Fredrick I. Olness:  olness@mail.physics.smu.edu
Robert L. Zimmerman: bob@zim.uoregon.edu

Copyright 2002, Addison-Wesley Publishing Company, Inc. The material in
this file may be distributed freely  so long as the content remains
unchanged, and the copyright and  reference notices are included. The
publisher grants permission for the noncommercial use of these programs
and program segments.  All other uses require the prior written consent
of the publisher.    

All rights are reserved with regard to the  material in the
"Mathematica For Physics" text, and may not be reproduced, stored in a
retrieval system, or transmitted,  in any form or by any means,
electronic, mechanical, photocopying,  recording, or otherwise, without
the prior written permission  of the publisher. 

The programs and applications presented in this book have been included
for  their instructional value. They have been tested with care but are
not  guaranteed for any particular purpose. The publisher does not
offer any  warranties or representations, nor does it accept any
liabilities with  respect to the programs or applications. 

 *)

5. Discrete Dynamical Systems

5.1 Introduction

Off[General::spell];  
Off[General::spell1];

5.2 Logistic Map

Overview of the Logistic Map

Problem 1: Logistic Map

Solution

Clear["Global`*"]

Part a

map[ μ_]  := μ # (1 - #) &;

{x0,  map[ μ][
        x0], map[μ][map[ μ][x0]]}

orbit[μ_, x0_, n_ ] := 
      NestList[ μ  # (1 - #) &, x0 , n ]

orbit[μ , x0 , 2 ]

orbit[2.9  , 0.1, 20 ]

Take[orbit[2.9, 0.1, 300 ], -8]

orbitplot[μ_, start_, n_  , 
        Opts___] := ListPlot[
        NestList[ μ  # (1 - #) &, start , 
          n ], 
Opts , 
PlotStyle -> {AbsolutePointSize[4],  RGBColor[1, 0, 0]} , 
        GridLines -> 
          Automatic,  
Epilog -> ListPlot[
              NestList[μ  # (1 - #) &, start, 
                n], 
PlotJoined -> 
                True , 
DisplayFunction -> 
                Identity
][[1]] 
]

orbitplot[2.9 , 0.1, 50, 
        AxesLabel -> {n,  (Subscript[x, i])} ];

Part b

list = 
        orbitplot[3.4, 0.1, 
          50,  
PlotRange -> {Automatic, {0, 
                1}},    
AxesLabel -> {n,  
              (Subscript[x, i])}];

eq1= Take[
        orbit[3.4,  .1, 200 ], -8]

Union[eq1, 
      SameTest -> (Abs[( #2  -  #1 )] < 
              (10^(-12)) &)]

list = 
        orbitplot[3.54,  .1, 
          50, 
PlotRange -> {Automatic, {0, 
                1}} ,  
AxesLabel -> {n,  
              (Subscript[x, i])}];

eq2= Take[
        orbit[3.54, 0.1, 200 ], -8]

Union[eq2, 
      SameTest -> (Abs[(#2 - #1)] < (10^(-4)) &)]

list = 
        orbitplot[3.55,  .1, 
          50, 
PlotRange -> {Automatic, {0, 
                1}}, 
AxesLabel -> {n,  
              (Subscript[x, i])}];

eq3= Take[
        orbit[3.55, 0.1, 1000 ], -16]

Union[eq3, 
      SameTest -> (Abs[(#2 - #1)] < (10^(-12)) &)]

list = 
        orbitplot[3.7,  .1, 200, 
          PlotRange -> {Automatic, {0, 1}}];

Problem 2: Logistic Fixed Points

Solution

Clear["Global`*"]

Part a

period[
        1] = μ x0 (1 - x0) == x0 // Solve[#, x0] &   // 
        Flatten

Nest[ μ  # (1 - #) &, x0, 2 ] == x0

eq1 =  Nest[ μ  # (1 - #) &, x0, 2 ] == x0 // 
            Solve[#, x0] &  // FullSimplify  // Flatten

period[2] = Complement[eq1, period[1]]

fixedpoint[μ_, n_, 
        opts___] :=
Module[{temp1, 
          x0}, 
temp1 = 
          Nest[ μ  # (1 - #) &, x0, n ] == x0 // 
            NSolve[#, x0, opts] & ;    Select[x0 /. temp1, 
          Abs[Im[#]] < 
              (10^(-10))  &]  
]

{pt1, pt2, pt3, pt4} =   fixedpoint[7/2, # ] & /@ Range[4] // 
          Rationalize;
{pt1, pt2, pt3, pt4} // ColumnForm

Complement[ pt2, pt1 ,  
      SameTest -> (Abs[( #2  -  #1 )] < 
              (10^(-3)) &)]

Complement[ pt4,  pt2 ,  
      SameTest -> (Abs[( #2  -  #1 )] < 
              (10^(-4)) &)]

Part b

Plot[{7/2  # (1 - #) &[x] , x} , {x, 0, 
          1}, 
PlotStyle -> {Thickness[0.01], 
            RGBColor[1, 0, 0]}, 
GridLines -> 
          Automatic, 
        Epilog -> {PointSize[0.04], Hue[0.3], 
            Map[Point, 
              Transpose[{ fixedpoint[7/2, 1], 
                  fixedpoint[7/2, 1]}]]}
];

fixedpointG[μ_, n_,  Opts___]  := 
      Plot[{x, Nest[ μ # (1 - #) &, x, n ]}  // Evaluate, {x, 0, 
          1}, 
Opts, 
PlotStyle -> {{Thickness[0.01],  Hue[Random[]]}, {Hue[Random[]]}},  
        DisplayFunction -> 
          Identity , 
Epilog -> {PointSize[0.04], 
            Hue[0.3], 
            Map[Point, 
              Transpose[{ fixedpoint[μ, n], 
                  fixedpoint[μ, n]} ] ]}
]

list = fixedpointG[7/2, # , 
              PlotLabel -> #] &   /@   Range[4, 
            1, -1] ;

Show[GraphicsArray[Partition[list, 2]]];

Show[ list , 
        DisplayFunction -> $DisplayFunction];

Part c

list2 = fixedpointG[#, 
              2,  
Ticks -> {{0, 1/2, 1} , {0, 1/2, 
                    1}},  
PlotRange -> {{0, 1}, {0, 
                    1}}, 
PlotLabel -> # ]  &    /@    { 2.5, 3,  3.5}  ;

Show[ 
        GraphicsArray[{{list2[[1]]}, {list2[[2]] ,  
              list2[[3]]}} ] , 
        PlotRange -> All ];

Problem 3: Logistic Cobwebs.

Solution

Clear["Global`*"]

Part a

cobweb[f_ , x0_, n_,  
        opts___]  := 
Module[ {temp1, temp2, points, s , 
          x}, 
temp1[s_]  :=  
          Nest[f[#]  &, x0 , 
            s ];  
temp2  =   {temp1[# ], 
                temp1[# ], temp1[# ], temp1[# + 1]} & /@ 
            Range[1, 
              n] ; 
points = {x0, 0, x0, f[x0], 
                temp2 } // Flatten  // 
            Partition[#, 2] &  ; 
Plot[{x, f[x]},  {x, 
            0, 1},  Epilog -> {Hue[ .9], Line[points]}, 
          opts]
]

cobweb[7/2 # (1 - #) &,  5 /7, 30 ] ;

list =  cobweb[7/2 # (1 - #) &,  # , 10 ,  
              DisplayFunction -> Identity ] & /@ { 3/7,  
            1/2 };

Show[GraphicsArray[ list  ]];

Part b

cobweb[   3.99 # (1 - #) &,  1/2 , 150  ];

Problem 4: Logistic Bifurcations

Clear["Global`*"]

Solution

Lvalues = 
        Compile[{μ}, ({μ, #} &)  /@  
            Drop[NestList[μ # (1 - #) &,  .62, 100],  
              75]];

bifurcation[μ0_,μ1_,n_]:= Flatten[Table[Lvalues[μ],{μ,μ0,μ1,(μ1-μ0)/n}],1]

Spoints = bifurcation[3,3.99,400];

pt1 =  
        ListPlot[Spoints , 
          GridLines -> 
            Automatic, 
PlotStyle -> {Hue[ .9], 
              AbsolutePointSize[
                1.4]}, 
Epilog -> {
Thickness[ .007],  
Line[{{# , 
                        0}, { # , 1}}] &   /@   {3.45, 3.54, 
                  3.564 }}
] ;

Problem 5: Logistic Lyapunov Exponent and Entropy

Solution

Clear["Global`*"]

Part a

f[μ_]  := μ*# (1 - #) &;

LGder[μ_] := Log[
          Abs[μ* (1 - 2 #)]] &

data = 
        Table[
{μ, 
            1/1000   LGder[μ ][ 
                  NestList[μ*# (1 - #) &, 0.1, 999]] /. 
              List -> Plus },    {μ, 3.5, 4.0, 
            0.001}];

pt1 = 
        ListPlot[
          data, 
PlotRange -> {-1 , 
              1}, 
PlotJoined -> 
            True, 
PlotStyle -> 
            Thickness[ .001],  
AxesLabel -> {μ, λ}];

Part b

Entropy[
        p_List] := - p . Log[p]

fdata[μ_] := 
      Drop[NestList[μ # (1 - #) & ,  .1, 500], 400]

Needs["Statistics`DataManipulation`"]

?BinCounts

prob[μ_]:=Select[
        1/500  BinCounts[fdata[μ] , {0.0, 1, 1/500 }] , 
        Positive]

entropylist = 
        Table[{μ, Entropy[prob[μ]] }, {μ, 3.5, 4, 
            0.005}];

\(\(\(pt2 = 
        ListPlot[
          entropylist, 
PlotJoined -> 
            True, 
GridLines -> 
            Automatic, 
PlotStyle -> 
            Thickness[ .005], \ 
AxesLabel -> {μ, "Entropy"}];\)\(\ \ \)\)\)

Show[pt1, 
        pt2, 
PlotRange -> {{3.5, 3.9}, {-1.5, 
              1.5}}, 
Ticks -> {Automatic, 
            None}];

5.3 Other Maps

Overview of Other Maps

Problem 1: Salmon Map

Solution

Clear["Global`*"]

Part a

map[μ_] := #  Exp[μ (1 - #)] &

orbitG[μ_, start_, n_ , Opts___]:=ListPlot[
        NestList[map[μ ], start , n ], Opts ,  
        DisplayFunction -> Identity]

values = { 1.8 , 2.2 , 2.6, 2.99};

list = orbitG[#,  .4, 
              50, 
PlotStyle -> 
                AbsolutePointSize[2] , 
Ticks -> 
                False , 
PlotLabel -> # 
] &   /@   values ;

Show[
        GraphicsArray [  {list[[{1, 2}]],  
            list[[{3, 4}]]}  ]];

Take[NestList[map[#], 0.5 , 200 ], -4] & /@ { 1.8 , 
          2.2  , 2.6}//ColumnForm

Part b

Svalues = Compile[{μ},({μ,#}&) /@ Drop[NestList[ # Exp[μ(1 - #)]&, 
0.5, 200],175]];

bifurcation[μ0_, μ1_, n_] := 
  Flatten[Table[Svalues[μ],{μ, μ0, μ1, (μ1 - μ0)/n}], 1]

Spoints = bifurcation[1.9,3.5,1000];

pt1 = 
        ListPlot[
          Spoints, 
GridLines -> 
            Automatic, 
PlotStyle -> {Hue[  .9] , 
              AbsolutePointSize[ .5]} ] ;

Part c

f[μ_] := # Exp[μ (1 - #)] &;

    
fder[μ_] :=  
      Log[Abs[(1 - μ #) Exp[μ (1 - #)]]] &

points[n0_]  := 
        Table[{μ, 
            Apply[Plus, fder[μ ][NestList[f[μ ], 0.1, n0 - 1]]/
                n0]}, {μ, 2 , 3.5, 1/n0}];

pt2 = 
        ListPlot[
          points[300] ,  
PlotJoined -> 
            True, 
PlotStyle -> 
            Thickness[ .006]];

pt3 = Show[pt1, pt2];

Part d

Needs["Statistics`NormalDistribution`"]

?NormalDistribution

WNdata=Table[Random[NormalDistribution[1,0.5]], {500}];

chaos=NestList[map[2.9],0.95,500];

Show[
GraphicsArray[
ListPlot[#, 
PlotStyle -> 
                  AbsolutePointSize[2], 
 
                Ticks -> {{100, 300, 500}, {0,  1, 
                      2}},  
DisplayFunction -> 
                  Identity
] &  /@  {WNdata, 
              chaos} 
]
];

Show[
        GraphicsArray[ListPlot[
                Transpose[{Drop[#1, -1], 
                    Drop[#1, 1]}], 
PlotRange -> 
                  All, 
AxesOrigin -> {0, 
                    0}, 
PlotStyle -> 
                  AbsolutePointSize[
                    2], 
DisplayFunction -> 
                  Identity
] &   /@  {WNdata, 
              chaos}
]
];

Problem 2: Sine-Circle Map

Clear["Global`*"]

Solution

Part a

map[ Ω_, K_]  := 
      Mod[( # + Ω - K/(2 π) Sin[2 π #]), 
          1 ] &

cobweb[f_ ,  x0_, n_,  
        opts___]  :=  
Module[ {temp1, temp2, points, 
          s , x}, 
temp1[s_]  :=  
          Nest[f[#]  &, x0 , 
            s ]; 
temp2  =   {temp1[# ], temp1[# ], 
                temp1[# ], temp1[# + 1]} & /@ 
            Range[1, 
              n] ; 
points = {x0, 0, x0, f[x0], 
                temp2 } // Flatten  // 
            Partition[#, 2] &  ; 
Plot[{x, f[x]},  {x, 
            0, 1},  Epilog -> {Hue[ .9], Line[points]}, 
          opts]
]

cobweb[map[1/7, 0], 1/10 , 50 ] ;

Take[
        NestList[map[1/7, 0], 1/10 , 
          400]  ]//Union

cobweb[map[Sqrt[2], 0], 1/10 , 50 ] ;

list =  
cobweb[map[#, 0.9] , 0.45, 
              50,  
Ticks -> {{0,  .5, 1}, 
                  Automatic}, 
PlotLabel -> #, 
              DisplayFunction -> 
                Identity 
] &   /@  {0.04, 0.61, 
            0.65}  ;

Show[
        GraphicsArray[{{list[[1]]}, list[[{2, 3}]]}]];

Take[
      NestList[map[0.04, 0.9 ], 0.45 , 200], -9 ]

Take[
        NestList[map[0.65, 0.9 ], 0.45, 
          200], -9 ]//Union

Part b

orbit  =  
Compile[{K}, 
Map[{K, #} &,  
Drop[ 
                NestList[ 
                  Mod[( # +  .65 - 
                          K/(2 π) Sin[2 π #] ) , 
                      1 ] & , 0.1, 100],  
                75]
]
] ;

bifurcation[K0_, K1_, n_]  :=  
        Flatten[Table[orbit[K], {K, K0, K1, (K1 - K0)/n}], 1]

ListPlot[ bifurcation[0,6,800],PlotStyle -> AbsolutePointSize[1.1]];

Part c

wind[ Ω_, 
        K_ , θ_] := 
Module[{temp1 }, 
temp1  =  
          NestList[( # + Ω - 
                  K/(2 π) Sin[2 π #] ) &, θ , 
            699 ];  
(1/600 ) ( 
            Last[ temp1] - First[temp1])
]

pt1 =  
        ListPlot[{#, wind[#, 1.2  ,  .3]} & /@ 
            Range[ .25 ,  .75,  .008], 
PlotJoined -> 
            True, 
PlotStyle -> {Thickness[ .006], 
              Hue[ .6]}, 
GridLines -> 
            Automatic , 
AspectRatio -> 
            1 ];

Problem 3: Taylor-Greene-Chirikov map

Solution

Clear["Global`*"]

Part a

map[k_][{θ_, i_}] = { 
          Mod[θ + i + k  Sin[θ], 2 π], 
          Mod[i + k Sin[θ], 2 π]};

orbit[θ0_, i0_] := 
      NestList[map[0.97  ], {θ0, i0}, 150]

\(orbitG[θ0_, i0_, Opts___]\  := \ 
ListPlot[\ 
        orbit[\ θ0, i0\ ]\ , 
Opts\ , 
        Frame -> True, Axes -> None, 
        FrameLabel -> {"θ", "i"}\ , 
PlotStyle -> {Hue[Random[]], 
            PointSize[ .009\ ]}\ , 
DisplayFunction -> 
          Identity
]\)

Show[ orbitG[Random[Real, {0, 2 π}], 
              Random[Real, {0, 2 π}]] & /@ 
          Range[60 ], 
DisplayFunction -> $DisplayFunction, 
GridLines -> 
          Automatic, 
AspectRatio -> 
          Automatic];

Part b

Chθ[{θ0_, i0_}] := 
      Map[First, orbit[θ0, i0]]

pt1 = ListPlot[
              Take[Chθ[# ], -40] , 
PlotJoined -> True, 
AxesLabel -> { 
                  n, θ }, 
DisplayFunction -> 
                Identity
] &   /@   {{ 3,  
              0.1 }, {0.65, 3}} ;

Show[GraphicsArray[{pt1 }]];

Part c

wind[  K_ , θ_, 
        i_] :=  
Module[{map1, orbit2, wind3,  θa, 
          pb}, 
map1[{θa_, pb_}] := { θa + 
              pb + K  Sin[ θa], 
            pb + K  Sin[ θa]} ; 
orbit2 = 
          NestList[map1, {θ , i }, 
            400  ] ; 
wind3 =  
          Map[First, orbit2]  ;  
(1/400 ) ( 
            Last[ wind3] - First[wind3])
]

ListPlot[{#, wind[ 0.45  , 0.2 , #]} & /@ 
          Range[0.15, 0.24, 0.005], 
PlotJoined -> 
          True, 
PlotStyle -> {Thickness[0.005], 
            Hue[0.6]}, 
PlotLabel -> {"K=0.47,\
θ=0.2"}, 
AxesLabel  -> {i, 
            None }];

Problem 4: Henon Map

Solution

Clear["Global`*"]

Part a

map[a_, b_][{x_, y_}]:={1 - a (x^2) + y , 
        b x  }

pt1 = 
ListPlot[ 
NestList[
            map[1.4,  0.3], { 0 , 0.1}, 
            20000 ], 
PlotStyle -> 
            AbsolutePointSize[2],  
GridLines -> 
            Automatic ];

Show[pt1, PlotRange -> {{0.3, 0.47}, {0.2, 0.217}}, 
        GridLines -> Automatic];

Part b

Hx[a_] :=   Map[First, 
        Drop[ NestList[  map[a, 0.3 ] , {0.3, 0.6}, 200], 
          150] ]

bifurcation[a0_, a1_, n_]  :=  
        Flatten[Table[Hx[a], {a, a0, a1, (a1 - a0)/n}],  1]

pts = bifurcation[1.05 , 1.1, 300 ] ;

ListPlot[
        pts, 
PlotStyle -> 
          AbsolutePointSize[1.4], 
AspectRatio -> 
          1, 
Ticks ->  {{}, 
            Automatic }];

5.4 Fractals

Overview of Fractals

Problem 1:  Mandelbrot Set

Solution

Clear["Global`*"]

Part a

MandO = 
        Compile[{{c, _Complex}},  
          FixedPointList[#^2 + c &, 0, 69, 
            SameTest -> (Abs[#2] > 2.0 &)]] ;

Take[MandO[1/10], -1] // Chop

MandO[ 1 ] // Chop

Part b

Mandelbrot[c_] := Length[MandO[c]]

Plot3D[
        Mandelbrot[x + I y], 
{x, -2.0, 
          0.5}, {y, -1.1, 1.1},  
AspectRatio -> 
          Automatic,  
ColorFunction -> 
          Hue,   
PlotRange -> 
          All , 
PlotPoints -> 100];

DensityPlot[ 
        Mandelbrot[
          x + I y], 
{x, -2,  .5} , {y, -1, 
          1},  
Mesh -> 
          False,   
AspectRatio -> 
          Automatic, 
Frame -> 
          True, 
ColorFunction -> 
          Hue, 
PlotPoints -> 100 ];

DensityPlot[
        Mandelbrot[
          x + I y], 
{x, -1.45, -1.25}, {y, -0.12, 0.12} ,  
Mesh -> 
          False,   
AspectRatio -> 
          Automatic, 
Frame -> 
          True, 
ColorFunction -> 
          Hue, 
PlotPoints -> 
          100];

Problem 2:  Julia Set

Solution

Clear["Global`*"]

Part a

Juliao = 
        Compile[{{z0, _Complex}},  
          FixedPointList[#^2 + c &, z0, 50, 
            SameTest -> (Abs[#2] > 2.0 &)]] ;

c = -1; 
Juliao[0.1 ] // Chop

Juliao[1.9 ] // Chop

Part b

Julia = 
        Compile[{{z0, _Complex}}, 
          Length[FixedPointList[#^2 + c &, z0, 50, 
              SameTest -> (Abs[#2] > 2.0 &)]]];

c =  -1;

Plot3D[
        Julia[x + I y], 
{x, -2.0, 2}, {y, -1 , 
          1 }, 
AspectRatio -> 
          Automatic,  
ColorFunction -> 
          Hue , 
PlotPoints -> 
          100,  
PlotRange -> All ];

DensityPlot[
        Julia[x + I y], 
{x, -2, 2}, {y, -1 , 
          1 },   
Mesh -> 
          False,   
AspectRatio -> 
          Automatic, 
Frame -> 
          True, 
PlotPoints -> 
          100, 
ColorFunction -> Hue ];

5.5 Exercises

Exercise 1: Gaussian Map

Problem 2: Return Maps for the Henon Map

Exercise 3: Cubic Maps

Exercise 4: Two Dimensional Map

Exercise  5: Tent map