This example is a complete proof of the existence of 116 periodic orbits for the Lorenz system with chaotic parameter values s=10, r=28., q=8/3.
This example shows also how to speed up computations by means of static memory allocation.
Allocating objects on storage (malloc or new) is time consuming. If the dimension of the phase space is fixed and known at compile time you can use static memory allocation like
which is much faster than malloc or new. The CAPD library provides a header file capd/fdcapdlib.h that defines most important types for the dimension specified by the user.
In this example we will use 3-dimensional objects for integration of the Lorenz system and variable size vectors and matrices when computing Interval Newton Operator (this is not time-critical operation in this example).
#include <iostream>
#include <sstream>
#include "capd/capdlib.h"
#define CAPD_DEFAULT_DIMENSION 3
#define CAPD_USER_NAMESPACE capd3
#include "capd/fdcapdlib.h"
#undef CAPD_DEFAULT_DIMENSION
#undef CAPD_USER_NAMESPACE
#include "LorenzPeriodicOrbits.dat"
bool provePeriodicOrbit(capd3::IPoincareMap& pm,
interval X[],
int dim)
{
IVector imCenter(dim), Y(dim,X);
IMatrix D(dim,dim);
for(int i=0;i<dim;i+=2){
int prev = (dim+i-2)%dim;
imCenter[i] = center[i] - x[0];
imCenter[i+1] = center[i+1] - x[1];
x = pm(s2,DP);
DP = pm.computeDP(x,DP);
D[i][i] = 1;
D[i + 1][i + 1] = 1;
D[i][prev] -= DP[0][0];
D[i][prev + 1] -= DP[0][1];
D[i + 1][prev] -= DP[1][0];
D[i + 1][prev + 1] -= DP[1][1];
}
}
{
double tolerance = 1e-7;
try{
capd3::IMap vectorField(
"par:q;var:x,y,z;fun:10*(y-x),x*(28-z)-y,x*y-q*z;");
capd3::IOdeSolver solver(vectorField,
order);
capd3::ICoordinateSection section(3,2,27.);
capd3::IPoincareMap pm(solver,section);
solver.setAbsoluteTolerance(tolerance);
solver.setRelativeTolerance(tolerance);
const int numberOfPeriodicOrbits = 116;
double xMin, xMax, yMin,yMax, temp;
int period, validated=0;
for(int i = 0;i<numberOfPeriodicOrbits;++i){
istringstream in(lorenzPOData[i]);
in >> period;
period*=2;
for(int j=0;j<period;++j){
in >> xMin >> xMax >> yMin >> yMax >> temp >> temp;
}
validated += provePeriodicOrbit(pm,X,2*period);
cout << "already validated: " << validated << " out of " << (i+1) << endl;
}
} catch(exception& e) {
cout << "\n\nException caught: "<< e.what() << endl;
}
return 0;
}
capd::dynset::C0HOSet< CAPD_USER_NAMESPACE::C0Rect2Set > C0HORect2Set
Definition: typedefs.h:48
capd::dynset::C1DoubletonSet< CAPD_USER_NAMESPACE::IMatrix, C1Rect2Policies > C1Rect2Set
Definition: typedefs.h:54
#define main()
Definition: krak-lib.h:388
capd::map::Map< IMatrix > IMap
Definition: mapLib.h:24
int order
Definition: tayltst.cpp:31
bool subsetInterior(const IntervalObject &v1, const IntervalObject &v2)
checks if IntervalObject v1 is contained in interior of IntervalObject v2
Definition: iobject.hpp:143
capd::vectalg::Vector< interval, 3 > IVector
Definition: vecttst.cpp:28
capd::vectalg::Matrix< interval, 3, 3 > IMatrix
Definition: vecttst.cpp:29
MatrixType krawczykInverse(const MatrixType &A)
Definition: Matrix_Interval.hpp:53
GeometricBound midVector(const GeometricBound &x)
Definition: GeometricBound.h:273
This class provides a trait of being set of a given type, i.e.
Definition: DagIndexer.h:22