Here is an example of how to use the "d1GTO_s" class to compute an attracting invariant manifold of a dynamical system for different parameters. The initial invariant manifold is a stable invariant closed curve of an ODE (the Rossler system).
// This example computes a stable periodic orbit of the Rossler equation
// and an invariant transversal to the orbit.
#include "d1GTO_s.h"
#include "d1LocInvTrans_s.h"
#include "timeDisc.h"
const int itn = 3; // itn = the ambient dimension
REAL lambda = 2.2; // parameter in the Rossler equation
VEC evalFnp(VEC); // the Rossler equation
MAT derivFnp(VEC); // the space derivative of the Rossler equation
main()
{
int i,j,k;
// parameters for the time-discretization class:
const int iterates = 1;
const int hasInverse = 1;
int dn1MAXITS=1000;
REAL TOLF=1.0e-6,TOLMIN=1.0e-6,TOLX=1.0e-7,STPMX=100.0;
REAL timeStep = 0.5;
// explicit constructor for the time-discretization class:
timeDisc *timeDPtr = new timeDisc(hasInverse,itn,iterates,timeStep,"rk4",evalFnp,
derivFnp,dn1MAXITS,TOLF,TOLMIN,TOLX,STPMX);
// also possible:
// new timeDisc(hasInverse,itn,iterates,timeStep,"rk4",evalFnp,derivFnp)
/*
Below we construct an initial "d1smallGrid", which is a grid of points (nodes)
in ambient space, and at each point a vector in ambient space indicating
the tangent direction. This grid represents a 1-dimensional manifold.
The nodes are partitioned into a number of "elements", which all contain
the same number of nodes. Later we will fit a polynomial to the points
on each element.
*/
const int itm = 3; // the number of nodes per element
const int itN = 32; // the number of elements
const int itNodes = (itm-1)*itN + 1; // the total number of nodes
// itpnt is a vector of vectors, representing all the points of the grid:
baseVector <VEC> itpnt(itNodes);
// itTE is a vector of matrices, representing all the tangent data for the grid:
baseVector <MAT> itTE(itNodes);
// We will combine the data of the previous two lines and put it in this structure,
// which is a matrix of d1baseGridPts (Each d1baseGridPt is a point with an
// attached vector in ambient space.) Each row of the matrix represents an element.
baseMatrix<d1baseGridPt> gridPnts(itN+1, itm-1);
REAL r = 2.0*pow(2.0,0.5); // parameter for the initial data
// A circle in 3-dimensions with radius r:
for(i = 0; i < itNodes; i++) {
itpnt[i] = *(new VEC(itn));
itpnt[i][0] = r*sin(2*Pi*i/(itNodes-1));
itpnt[i][1] = r*cos(2*Pi*i/(itNodes-1));
itpnt[i][2] = 0.0;
}
// The tangent data is specified merely by taking the difference
// of two consecutive points, and normalizing the result:
REAL norm; // for normalizing the tangent data
VEC diff(itn); // will hold the difference of two consecutive points
for(i = 0; i < itNodes-1; i++) {
itTE[i] = *(new MAT(itn,1));
diff = itpnt[i+1] - itpnt[i];
norm = diff.norm();
itTE[i].setColumn(0,(1.0/norm)*diff);
}
// Below we use the data specified above to fill in the matrix
// of d1baseGridPts called "gridPnts". The address of this matrix
// will be passed to the explicit constructor of the d1smallGrid class.
for (i=0; i < itN; i++)
for (j=0; j<itm-1; j++) {
k = i*(itm-1) + j;
gridPnts(i,j) = *(new d1baseGridPt(itn, itpnt[k], itTE[k]));
}
gridPnts(itN,0) = gridPnts(0,0); // first node equals last node (a closed curve).
// The explicit constructor for the "d1smallGrid":
d1smallGrid *grid2Ptr = new d1smallGrid(itn,itNodes,itN,itm,&gridPnts);
/*
Here we create a "d1Grid" corresponding to our "d1smallGrid". The
difference between a d1Grid and a d1smallGrid is that a d1Grid has
at each point an orthonormal set spanning the hyperplane normal to
the tangent vector. The constructor below automatically fills in
this data. normTol determines how close to "normal" this set is.
*/
REAL normTol = 1.0E-4;
d1Grid *gridPtr = new d1Grid(*grid2Ptr,normTol);
/*
Now we create a d1baseMan object, which represents our initial manifold,
over which we will iterate the algorithm. polyType is the type of
polynomial to use to represent the manifold, for example "Lagrange"
or "Hermite".
*/
static char polyType[9] = "Lagrange";
d1baseMan *baseMnPtr = new d1baseMan(polyType,itm,itN,gridPtr);
/*
Here we create a d1graphMan object, which represents a manifold obtained
at a step of our algorithm which is graphed over the d1baseMan.
*/
d1graphMan *gmanPtr = new d1graphMan(polyType,baseMnPtr,grid2Ptr);
/*
The "d1GTO_s" class represents the interface to the graph-transform
algorithm, applied to 1-dimensional manifolds, using orthogonal
coordinates, in the case of a stable manifold.
The explicit constructor's arguments are a pointer to a map representing
the dynamical system, and a pointer to a "d1graphMan".
*/
d1GTO_s *graphTransPtr = new d1GTO_s(timeDPtr, gmanPtr);
// parameters for the d1GTO_s class "iterate" method:
VEC gtAcc(2);
gtAcc[0] = 1.0E-4;
gtAcc[1] = 10.0;
int whichAlg=2,fullOutput=1,gtMaxit = 11;
int MAXIT = 100,mesh = 1,maxSearch = 100;
REAL maxDist = 100.0,xacc = 1.0E-6,projectNeigh = 1.0;
// The d1GTO_s class "iterate" method needs to be passed a d1graphMan,
// in addition to the one used to construct the class,
// so it can use one as the "current" manifold of the step, and one as
// the "next" manifold of the step:
d1graphMan *nextGManPtr = new d1graphMan(polyType,itn,itm,itN);
cout << graphTransPtr->iterate(whichAlg,fullOutput,gtMaxit,MAXIT,mesh,maxSearch,
gtAcc,maxDist,xacc,projectNeigh,nextGManPtr,"weg7s") << endl;
// Set the tangent data of the d1graphMan *gmanPtr, by taking
// derivatives of it's defining polynomials:
gmanPtr->setTangentData();
// Write the data of the d1graphMan *gmanPtr to a file called "smallGrid1.dat",
// so that the work done in the above iteration doesn't need to be re-done:
gmanPtr->GridPtr()->writeToFile("smallGrid1.dat");
// Read data determining a "d1smallGrid" from the file "smallGrid1.dat",
// to be used to create a new d1baseMan with which to restart the algorithm
grid2Ptr->readFromFile("smallGrid1.dat");
gridPtr = new d1Grid(*grid2Ptr,normTol);
baseMnPtr = new d1baseMan(polyType,itm,itN,gridPtr);
gmanPtr = new d1graphMan(polyType,baseMnPtr,grid2Ptr);
// At this point the d1baseMan *baseMnPtr and the d1graphMan *gmanPtr
// represent the same manifold.
graphTransPtr = new d1GTO_s(timeDPtr, gmanPtr);
nextGManPtr = new d1graphMan(polyType,itn,itm,itN);
gtMaxit = 100;
cout << graphTransPtr->iterate(whichAlg,fullOutput,gtMaxit,MAXIT,mesh,maxSearch,
gtAcc,maxDist,xacc,projectNeigh,nextGManPtr,"weg7t") << endl;
gmanPtr->setTangentData();
gmanPtr->GridPtr()->writeToFile("smallGrid1.dat");
// Now we will take the invariant manifold obtained above, and try to
// find the invariant stable bundle.
grid2Ptr->readFromFile("smallGrid1.dat");
gridPtr = new d1Grid(*grid2Ptr,normTol);
baseMnPtr = new d1baseMan(polyType,itm,itN,gridPtr);
// The "vecBundle" class is a matrix of matrices, representing,
// at each grid point, a matrix whose columns form an orthonormal
// set spanning a subspace of the ambient space. *sBndlPtr below will
// represent the stable bundle:
vecBundle *sBndlPtr = new vecBundle(itN,itm,1,itn-1);
// The "vecBundle" class provides for interpolation of the subspaces.
// Here we set the polynomials used for interpolation:
sBndlPtr->setPolynomials(baseMnPtr);
int haveCCS = 1;
d1LocInvTrans_s *sBndlLocPtr = new d1LocInvTrans_s(haveCCS,timeDPtr,baseMnPtr,
sBndlPtr);
// parameters for the "d1LocInvTrans_s" iterate_alg1 method (explained below):
REAL gtacc = 1.0E-6;
REAL contrTol = 50.0;
REAL lipTol = 2.0;
cout << sBndlLocPtr->iterate_alg1(mesh,maxSearch,MAXIT,xacc,gtMaxit,
gtacc,contrTol,lipTol) << endl;
return 0;
}
// Rossler equation.
VEC evalFnp(VEC arg)
{
assert(arg.getSize() == itn);
VEC val(itn);
val[0] = - arg[1] - arg[2];
val[1] = arg[0] + 0.2*arg[1];
val[2] = 0.2 + arg[2]*(arg[0] - lambda);
return val;
}
MAT derivFnp(VEC arg)
{
assert(arg.getSize() == itn);
MAT val(itn,itn);
val(0,0) = 0.0;
val(1,0) = 1.0;
val(2,0) = arg[2];
val(0,1) = -1.0;
val(1,1) = 0.2;
val(2,1) = 0.0;
val(0,2) = -1.0;
val(1,2) = 0.0;
val(2,2) = (arg[0] - lambda);
return val;
}
written by Dimitri van Heesch, © 1997-1998