This document describes computer programs for computing periodic orbits of the Lorenz equations as well as data about the computed periodic orbits.
A good place to start browsing the code is the makefile.
The makefile used icl, the Intel compiler and linker. Switching to gcc will take some work but I hope not too difficult.
The code is linked against Intel MKL for LAPACK and BLAS. Those libraries are used in LinAlg.cpp. If you link against some other version of LAPACK/BLAS, you will have to change the headers in LinAlg.cpp.
The code is also linked against GSL, the gnu scientific library. The library is used by Random.hh.
Finally, it is linked against FFTW2. This library is also used by LinAlg.cpp.
It appears as if this version was built and run on Windows (embarrassingly), which is why the object files are assumed to have the obj extension in the makefile. Moving it to Linux should be easy. There are no system dependencies.
I will now describe the source files in a bottom up fashion.
dgpadm.f is the fortran source for computing the matrix exponential that I downloaded from the Internet.
Vector.hh, Matrix.hh, and SparseMatrix.hh are implemented in LinAlg.cpp. Confusingly, SparseMatrix.hh does not define a class for sparse matrices. It is really just a slightly different interface for dense matrices.
Solve.hh defines a class for RK solvers implemented in Solve.cpp.
The core method for finding periodic orbits is PeriodicOrbit.hh implemented in PeriodicOrbit.cpp and POUtil.cpp. It takes in an initial guess, a vector field, and a jacobian as inputs and refines the initial guess into a more accurate periodic orbit.
Functions to manipulate symbol sequences and construct initial guesses are declared in lorenz1.hh and implemented in lorenz1.cpp.The saveOrbit() function saves a PO in a directory using an index (simply an int). The function loadOrbit() can load a PO of given index into the library (maintained in computer memory). The function loadOrbits() can load all POs within an index range. Given a symbol string, the function returnApprox() generates an initial guess for a PO of that symbol string using all the POs loaded into memory.
Functions to generate all POs with symbol sequences of a certain length or shorter are declared in lorenz2.hh and defined in lorenz2.cpp.
Functions to compute Lyapunov exponents are declared in lorenz3.hh and defined in lorenz3.cpp. The same sources also include functions to find the distance between two POs. Finding the distance between two POs involves careful shifting of phase.
Given any point on the Lorenz attractor, one can find POs arbitrarily close it. lorenz4.hh and lorenz4.cpp demonstrate that.
lorenz5.hh and lorenz5.cpp print some angles related to \(A^{n}B^{n}\) orbits.
The Lorenz vector field, Jacobian, and a function to compute a PO given a symbol string are declared in lorenz.hh and defined in lorenz.cpp
hdim.hh and hdim.cpp have code used to compute the Hausdorff dimension using periodic orbit theory. fredholm.math is a mathematica file used for the same purpose.
driver1.cpp has functions to demonstrate various properties of the Lorenz equations. Above each function, there are comment lines describing what it does. After several attempts, I was able to compute \(A^{25}B^{25}\). Theory predicts that \(A^{24}B\) should also exist. However, the method failed to find that orbit. I tried again several years later using a new idea for generating guess (see the function Lorenz4ex()) and this time the method worked.
Functions used to generate tables are found in tables.cpp.
Matlab scripts used to generate figures are found under FIGS. To run these scripts you first have to place the DATA folders in the locations the scripts expect them to be in.
This download is 69 MB. It includes code and data from multiple precision computations. The computations require the GNU mp library. Readme2.txt in that folder may be of help in browsing this code. This directory also includes ouput*.txt that show traces of certain runs, for example, of initially unsuccessful attempts to compute \(A^{24}B\).
This download is 73MB. Description of subdirectories:
These matlab scripts will be useful for loading and visualizing the data. Therefore we will describe some of them here.
scr1.m shows the basic way to load data:
load out.dat;
orbit = out;
clear in;
x = orbit(1,:);
y = orbit(2,:);
z = orbit(3,:);To load orbit17, replace out.dat by orbit17.dat. You can obtain its period from T17.dat. To plot x, y, z against time, find the length of x, y, z and divide T by that length to get h. \(h\) is the interval between successive points. The Lorenz vector field is encoded in the following vfield.m function:
function vfield = lorenz(orbit)
s = 10; b = 8/3; r = 28;
vfield = zeros(size(orbit));
x = orbit(1,:);
y = orbit(2,:);
z = orbit(3,:);
vfield(1,:) = s*(y-x);
vfield(2,:) = r*x - y - x.*z;
vfield(3,:) = x.*y - b*z; You can find the relative residual error of an orbit after loading it (along with its period) using the following function in relreserr.m:
function rerr = relreserr(orbit,T)
[m,n] = size(orbit);
orbitf = fft(orbit.').';
orbitf = orbitf.* (ones(3,1)*sqrt(-1)*[0,1:n/2-1,0,-n/2+1:-1]);
dorbit = real(ifft(orbitf.').');
w = 2*pi/T;
err = w*dorbit - vfield(orbit);
rerr = norm(err, 'fro')/ (w*norm(dorbit, 'fro'));To find the symbol sequence of an orbit after loading it, use the following symdyn.m function:
function s = symdyn(orbit)
cnt = 0;
[d,n] = size(orbit);
s = [];
for i = 1:n-1
if((orbit(3,i) <= 27) & (orbit(3,i+1)> 27))
cnt = cnt+1;
if(orbit(2,i) < 0)
s = [s,'A'];
else
s = [s,'B'];
end
end
end
if((orbit(3,n) <= 27) & (orbit(3,1)> 27))
cnt = cnt+1;
if(orbit(2,n) < 0)
s = [s,'A'];
else
s = [s,'B'];
end
endDATA1 contains data about 1375 periodic orbits corresponding to symbol sequences of length 13 of less. This is used as the library for all other computation of periodic orbits. A po numbered x is stored as orbitx.dat, its period as Tx.dat, the number of fourier points as nx.dat, and the error as errx.dat. If a lyapunov exponent is computed it is stored as lyapx.dat.
DATA2 contains data about periodic orbits corresponding to symbol sequences of length 20 or less. WinnerInfo has symbol sequence, period, error, and lyap exponent in that order. It also has data about periodic orbits with long and short periods and large and small Lyapunov exponents, for a given symbol length. err.dat has MAX and AVG error.
DATA3 contains data for 25 periodic orbits \(A^nB^n\). It has data for orbits of the form \(A^mB^n\) (index=100*m+n), with m = 25 and n = 24, 23, 22, 21, …7. If also has data to show \(A^{26}B^{26}\) pushed into a local minimum.
DATA4 contains data for 24 periodic orbits \(A^nB\)
DATA5 contains randomly generated pos with symbol strings of length 5 to 47 and their DISTANCE from a po with the same symbol string w/o the last symbol. DATA5/orbit347 = random po with sym string of length 347 DATA5/orbit110 = po \(AB...A^nB^n\) with n = 10 DATA5/orbit240 = po \(AB...A^nB^n\) with n = 20
DATA6 contains data about randomly generated points approximated by periodic orbits. trace*.dat has data from one type of algorithm for approximating the points. TRACE*.dat had data from another type of algorithm. The data in orbit*.dat are for the best-fit orbits computed by the second type of algorithm.
DATA7 contains data for plotting the Cantor book of the lorenz attractor. A.*, B.*, BA.* and AA.* have the data for the different sheets. scr1.m and scr2.m plot that data.
DATA8 contains the file PreciseInitCond.dat with \(100\) digits of the initial conditions for a few periodic orbits. DATA8/matlab contains routines related to data in DATA8/zoom(1/2/3). FIGS contains figures and tables.
DATA9 contains Hausdorff dimension data.
DATA10 had data about angles between stable and unstable directions of AnBn orbits.
The mathematica file used for periodic orbit theory computations is here: fredholm.math.