Chaos

Functions

We recall the midpoint RK method from Chapter 4:

RKODE[G_, x0_,f0_, dx_, Ns_] := Module[{k1,k2,retvals,index,xn,fn},
retvals = Table[{0.0,0.0},{i,1,Ns}];
xn = x0;
fn = f0;
retvals[[1]] = {x0,f0};
For[index = 2, index ≤ Ns, index=index+1,
k1 = dx G[xn,fn];
k2 = dx G[xn + dx/2, fn + k1/2];
fn = fn + k2;
xn = xn + dx;
retvals[[index]] = {xn,fn};
];
Return[retvals];
]

Getvals generates the sequence (Subscript[u, n+1])= Gup((Subscript[u, n])) given Gup, u0, and a parameter relevant to the map (in the case of the logistic map, G(u) = a u (1-u), the param would be “a”) - it takes Ns steps.

Getvals[Gup_, u0_, param_, Ns_] := Module[{uvals,index},
uvals = Table[0.0,{j,1,Ns}];
uvals[[1]] = u0;
For[index = 2, index ≤ Ns, index=index+1,
uvals[[index]] = Gup[uvals[[index-1]], param];
];
Return[uvals];
]

Sine and Henon Maps - Figures 13.1 and 13.2

We can run Getvals with the sine map:

sineu[un_,param_] := param Sin[un]

u1 = Getvals[sineu, .1,2.2,50];
u2 = Getvals[sineu,.1,2.3,50];

and plot both period one and period two stable solutions

G0 = ListPlot[u1,PlotRange->All,PlotStyle->GrayLevel[0.5]];
G1 = ListPlot[u2,PlotRange->All,PlotStyle->GrayLevel[0.5]];
GraphicsGrid[{{G0,G1}}]

The Henon map is a two-dimensional map, but Getvals can be used to generate its sequence:

Henonu[un_,param_] := {param[[1]] - param[[2]] un[[2]] - un[[1]]^2,un[[1]]}

uH = Getvals[Henonu,{.1,.1},{0,1},1000];
ListPlot[uH,PlotStyle->GrayLevel[0.5]]

Periodicity - Figures 13.3 and 13.4

The Logistic map is defined below:

Gup[u_,a_] := a u(1-u)

The period-one solutions can be found from the intersection of Gup with the line (Subscript[u, n+1])= (Subscript[u, n]) - the magnitude of the slope of Gup at the point of intersection tells us whether the solution is stable or not.

Plot[{u,Gup[u,1.0],Gup[u,2.0],Gup[u,3.0]},{u,0,1},PlotStyle->{GrayLevel[0],GrayLevel[0.25], GrayLevel[0.5], GrayLevel[0.75]}]

The period-two solutions can be found from the intersection of Gup applied twice, and the (Subscript[u, n+2])= (Subscript[u, n]) line.

Plot[{u,Gup[Gup[u,4.0],4.0]},{u,0,1},PlotStyle->{GrayLevel[0],GrayLevel[0.5]}]

Bifurcation Diagrams - Figure 13.5 and 13.6

This next function takes a set of values (inu, in vector form) from some map, and returns the unique values (up to tolerance eps) found in the last m entries.  It uses the built-in command “Union” to generate unique values using the “SameTest” defined below (i.e. two numbers are treated as the same if they are within eps of each other).

FindPeriods[inu_, m_, eps_] := Module[{Ns,udata,retvals},
Ns = Length[inu];
udata = Table[inu[[j]],{j,Ns-m+1,Ns}];
retvals = Union[udata, SameTest-> (Abs[#1-#2] ≤ eps&)];
Return[retvals];
]

Define G for the logistic map:

Gl[u_,a_] := a u (1-u)

Get one hundred values at a = 2 starting from u0 = 0.75:

ut = Getvals[Gl,0.75,2.0,200];

There is a stable fixed point at 1-1/a = 1/2, and if we look at the last ten entries of ut, we find a single value:

FindPeriods[ut, 100,10^(-8)]

A plot of the solution shows a small “transient” region (since we only look for periods in the last hundred data points, we miss this transient behavior in FindPeriods), where the value of u changes, before being attracted to the stable fixed point at .5.

G1 = ListPlot[ut,PlotStyle->GrayLevel[0.5]]

Moving up in a, we have a pair of stable, period-two solutions:

ut = Getvals[Gl, .1, 3.2,200];
FindPeriods[ut,100,10^(-8)]

G2 = ListPlot[ut,PlotStyle->GrayLevel[0.5]]

At a = 3.5, we have four period-four solutions:

ut = Getvals[Gl,0.8,3.5,200];
FindPeriods[ut, 100,10^(-8)]

and again, plotting the solution, we see a transient region, but ultimately, the four stable periodic points dominate.

G4 = ListPlot[ut,PlotStyle->GrayLevel[0.5]]

For a = 3.8, the periodicity of the solution exceeds our window of points (so the solution is, from a numerical point of view, either chaotic, or has period greater than 100)

ut = Getvals[Gl,0.8,3.8,200];
FindPeriods[ut, 100,10^(-8)]

Gc = ListPlot[ut,PlotStyle->GrayLevel[0.5]]

Plotting these four different types of solution together gives Figure 13.6.

GraphicsGrid[{{G1,G2},{G4,Gc}}]

Now we’ll make a bifurcation diagram - first generate a table with entries of the form {a, {unique limiting values} } for the logistic map with parameter a, then we’ll plot it.

avaltab = Table[{a, FindPeriods[Getvals[Gl, Pi/4., a, 1000], 100, 10^(-4)]},{a,1.1,4,.1}];

The next function takes the list of the values of a and all limit points, and returns a plot of these.

PlotFurcate[avals_] := Module[{ret,pts,index,entry,jndex},
pts = {PointSize[0.001]};
For[index = 1, index ≤ Length[avals], index=index+1,
entry = avals[[index]];
pts = Join[pts,Table[Point[{entry[[1]], entry[[2,j]]}],{j,1,Length[entry[[2]]]}]];
];
ret = Show[Graphics[{PointSize[.01/2],pts}],Axes->True,AxesOrigin->{1,0}];
Return[ret];
]

Graphics[PlotFurcate[avaltab]]

Now let’s increase the resolution:

avaltab = Table[{a, FindPeriods[Getvals[Gl, Pi/4., a, 1000], 100, 10^(-4)]},{a,1.1,4,.01}];

PlotFurcate[avaltab]

And again!  This is Figure 13.5.

avaltab = Table[{a, FindPeriods[Getvals[Gl, Pi/4., a, 1000], 100, 10^(-4)]},{a,1.1,4,.001}];

PlotFurcate[avaltab]

You can see, below, the values right around the bifurcation at a = 3 - to estimate the Feigenbaum number, we need to be able to reliably identify the bifurcation point - that is difficult numerically - for 2.993 and 3.002, we have a clear period-one and pair of period-two solutions, respectively.  But closer to the bifurcation point itself, both the size of our window (m, in FindPeriods, defining where we look so as to avoid counting the transient part of the solution) and the tolerance with which we identify unique values (eps, in FindPeriods) matter in determining exactly where to pinpoint the bifurcation.

avaltab[[1890;;1905]]

Initial Conditions - Figures 13.7, 13.8,13.9 and 13.10

The logistic map, again:

Gl[u_,a_] := a u (1-u)

In the non-chaotic regime, “all” starting points end at the same fixed point, and we’ll test this by starting at x0 = .1 -> .9 in steps of .1, this gives Figure 13.7:

Us = Table[Getvals[Gl,x0, 2.5, 10],{x0,.1,.9,.1}];
GG = Partition[Table[ListPlot[Us[[jj]],PlotRange->{0,1},PlotStyle->GrayLevel[0.5]],{jj,1,Length[Us]}],3];
GraphicsGrid[GG]

When the parameter a is in the chaotic regime, solutions do not converge to a fixed point (or set of fixed points with different periods), and the plot demonstrating this is Figure 13.8:

Us = Table[Getvals[Gl,x0,3.8, 10],{x0,.1,.9,.1}];
GG = Partition[Table[ListPlot[Us[[jj]],PlotRange->{0,1},PlotStyle->GrayLevel[0.5]],{jj,1,Length[Us]}],3];
GraphicsGrid[GG]

We can look at how solutions that are initially close con/di-verge as a function of step number, below we start off two solutions with slightly different initial values, x0, and plot the difference between the two solutions.  The a here is associated with a stable fixed point, so we expect the two solutions to both converge to the same fixed point, and their difference should go to zero:

diffseq = Getvals[Gl,.101,2.5,15] - Getvals[Gl,.1,2.5,15];
ListPlot[Abs[diffseq],PlotRange->All,PlotStyle->GrayLevel[0.5]]

In this next plot, we choose an a that is in the chaotic region, so the two solutions should diverge:

diffseq = Getvals[Gl,.101,3.8,15] - Getvals[Gl,.1,3.8,15];
ListPlot[Abs[diffseq],PlotRange->All,PlotStyle->GrayLevel[0.5]]

The pair of plots above form Figure 13.9.

The Lyapunov exponent can be estimated, from (13.38) using the function below - you provide the update map Gup, its derivative Gprime, the parameters (in param) and number of steps to take (Ns).  The function generates a solution starting from x0 = .1 (arbitrary, but it shouldn’t matter to the exponent calculation):

GetLyapunov[Gup_,Gprime_, param_,Ns_ ] := Module[{uvals,retval},
uvals = Getvals[Gup, .1, param, Ns];
retval = Sum[Log[Abs[Gprime[uvals[[j]], param]]],{j,1,Length[uvals]}]/Length[uvals];
Return[retval];
]

Using the logistic map, we have derivative:

DGl[u_,a_] := a-2 a u

and the Lyapunov exponent for a = 2.5 is:

ly = GetLyapunov[Gl,DGl, 2.5,10000]

The exponent will be negative when there are stable periodic solutions, and positive when the solutions are chaotic.  Below, we plot the log of the difference of two trajectories that start from nearby initial values (.1 versus .101) using a = 2.5, and we put the line with slope ly (the Lyapunov exponent we just calculated, for a = 2.5) on top, this is Figure 13.10.

diffseq = Log[Abs[Getvals[Gl,.101,2.5,15] - Getvals[Gl,.1,2.5,15]]];
G0 = ListPlot[diffseq,PlotRange->All,PlotStyle->GrayLevel[0.5]];
G1 = Plot[ly n + (diffseq[[5]] - ly 5),{n,1,15},PlotStyle->GrayLevel[0.5]];
Show[G0,G1]

For a = 3.8, the Lyapunov exponent is positive, indicating that this value for a gives chaotic solutions:

ly = GetLyapunov[Gl,DGl, 3.8,10000]

Lorenz Equations - Figures 13.11, 13.12, and 13.13

We’ll use the midpoint Runge-Kutta routine to solve the ODEs in (13.41) with parameters:

sigma = 10;
b = 8/3;

Here’s the right-hand side:

G[x_,f_] := {sigma (f[[2]] - f[[1]]), f[[1]](r-f[[3]]) - f[[2]], f[[1]] f[[2]] - b f[[3]]}

Run the solution out and plot the points in three dimensions - there is an attractive fixed point for r = 20:

r = 20;
fvals = RKODE[G,0.,{2.,5.,5.},.005,10000];
datapts = Table[Point[{fvals[[j,2,1]],fvals[[j,2,2]],fvals[[j,2,3]]}],{j,1,Length[fvals]}];
datapts = Join[{PointSize[0.01/4], datapts}];
Ga = Graphics3D[datapts,Axes->True]

We can start from two nearby initial values for x, say x0 = 2 and x0 = 2 + (10^(-7)), and calculate the difference in the resulting trajectories (using the norm, in three dimensions)

eps = 10^(-7);
fvals = RKODE[G,0.,{2.,5.,5.},.005,10000];
fvalse = RKODE[G,0.,{2.+eps,5.,5.},.005,10000];

res = Table[{fvals[[j,1]],Norm[fvalse[[j,2]] - fvals[[j,2]]]},{j,1,Length[fvals]}];
ListLinePlot[res,PlotRange->All,PlotStyle->GrayLevel[0.5]]

The trajectory and convergence plots form Figure 13.11.

Now the same experiment, but for r = 28, in the chaotic region for this problem - there is no stable fixed point, and the solution must continue to cycle around the two lobes:

r = 28;
fvals = RKODE[G,0.,{2.,5.,5.},.005,10000];
datapts = Table[Point[{fvals[[j,2,1]],fvals[[j,2,2]],fvals[[j,2,3]]}],{j,1,Length[fvals]}];
datapts = Join[{PointSize[0.01/4], datapts}];
Graphics3D[datapts,Axes->True]

We can again look at the divergence (this time) of initially close points:

eps = 10^(-7);
fvals = RKODE[G,0.,{2.,5.,5.},.005,10000];
fvalse = RKODE[G,0.,{2.+eps,5.,5.},.005,10000];

res = Table[{fvals[[j,1]],Norm[fvalse[[j,2]] - fvals[[j,2]]]},{j,1,Length[fvals]}];
ListLinePlot[res,PlotRange->All,PlotStyle->GrayLevel[0.5]]

the difference between the two trajectories is bounded, since the solutions each oscillate between the lobes.  The trajectory and divergence plot make up Figure 13.12.

Finally, we can project the chaotic solution into the two dimensional x-z plane with y = 0, that gives Figure 13.13:

r = 28;
fvals = RKODE[G,0.,{2.,5.,5.},.005,12000];

dataflat= Table[Point[{fvals[[j,2,1]],fvals[[j,2,3]]}],{j,1,Length[fvals]}];
GA = {{Graphics[{AbsolutePointSize[.1],dataflat[[1;;2000]]},Axes->True, AxesOrigin->{0,0}],Graphics[{AbsolutePointSize[.1],dataflat[[1;;4000]]},Axes->True, AxesOrigin->{0,0}],
Graphics[{AbsolutePointSize[.1],dataflat[[1;;6000]]},Axes->True, AxesOrigin->{0,0}]},{Graphics[{AbsolutePointSize[.1],dataflat[[1;;8000]]},Axes->True, AxesOrigin->{0,0}],
Graphics[{AbsolutePointSize[.1],dataflat[[1;;10000]]},Axes->True, AxesOrigin->{0,0}],Graphics[{AbsolutePointSize[.1],dataflat[[1;;12000]]},Axes->True, AxesOrigin->{0,0}]}};
GraphicsGrid[GA]

Mandelbrot set - Figures 13.14 and 13.15

To define the Mandelbrot set, we need to determine if a sequence has finite norm as n goes to infinity.  We can estimate the norm for “large” n by generating the sequence and checking the absolute value of the last element - the function below returns that absolute value, given n (the number of steps to take) and c, the starting location for the sequence (in the complex plane):

MandelIter[n_,c_] := Module[{val,index},
val = c;
For[index =1, index ≤ n, index=index+1,
val = val^2 + c;
If[Abs[val] ≥ INFINITY, index = n;];
];
Return[val];
]

INFINITY is defined to be (10^10) for our purposes in this problem (although you could just use 2 here).

INFINITY = 10^10;

This function generates a grid of points in the complex plane, and runs the sequence that defines the Mandelbrot set for n = 30 steps - if the norm of the 30th element is less than INFINITY, the point is included in the set, otherwise, the point is assigned the value of the log of the 30th element’s absolute value (for picture purposes):

MandelGrid[crange_]:=Module[{points,abs,MAX},
points = Table[{x,y,If[( abs = Abs[MandelIter[30,N[x+ I y]]]) ≤ INFINITY,-1.0,Log[abs]]},{x,crange[[1,1]], crange[[1,2]], (crange[[1,2]] - crange[[1,1]])/crange[[3]]},{y,crange[[2,1]], crange[[2,2]], (crange[[2,2]]-crange[[2,1]])/crange[[3]]}];
points = Partition[Flatten[points],3];
MAX = Max[Abs[Table[points[[jj,3]],{jj,1,Length[points]}]]];
Print[MAX];
points = Table[{points[[jj,1]],points[[jj,2]], If[points[[jj,3]] == -1.0, -1.0, points[[jj,3]]/MAX]},{jj,1,Length[points]}];
Return[points];
]

The range we’ll use is

crange = {{-2.5,1.0},{-1.5,1.5},500};

pts = MandelGrid[crange];

f3[x_] :=Abs[Log[x]]

ptsout = Table[If[pts[[jj,3]] == -1.0, {GrayLevel[0.0],Point[{pts[[jj,1]],pts[[jj,2]]}]}, {GrayLevel[f3[pts[[jj,3]]]],Point[{pts[[jj,1]],pts[[jj,2]]}]}],{jj,1,Length[pts]}];

To generate the grayscale, we take the Log again (via f3), this gives Figure 13.14

Graphics[ptsout,Axes->True]

Using a different range, we can probe different regions of the set (Figure 13.15):

crange = {{-.5,0.1},{0.5,1.2},500};

pts = MandelGrid[crange];

ptsout = Table[If[pts[[jj,3]] == -1.0, {GrayLevel[0.0],Point[{pts[[jj,1]],pts[[jj,2]]}]}, {GrayLevel[f3[pts[[jj,3]]]],Point[{pts[[jj,1]],pts[[jj,2]]}]}],{jj,1,Length[pts]}];

Graphics[ptsout,Axes->True]