Here is an example of how to use the "d2TCGTO_s" class to continue an attracting invariant manifold of a dynamical system over a parameter range. The initial invariant manifold is a stable invariant torus of a map (the "fattened sink map").
// In this example we continue an attracting invariant torus of the "fattened sink map",
// a map on R x S^{1} x S^{1}.
#include "d2TCGTO_s.h"
#include "map.h"
#include <math.h>
#include "fstream.h"
// parameters of the dynamical system
const REAL alpha = 0.5;
const REAL beta = 0.25;
// discretization parameters indicating size of grid
const int n1=50,n2=50;
const int twon1=2*n1;
int ambDim=3; // the ambient dimension
int order=1; // the order of discretization
int numNodes=(n1+1)*(n2+1); // the total number of grid points
int numCoarseElts=2*n1*n2; // the initial number of elements
const REAL twoPi=2.0*Pi;
// the dynamical system and its derivative
VEC evalFnp(VEC,REAL);
MAT derivFnp(VEC,REAL);
main()
{
int i,j,k;
// Create the dynamical system
int hasInverse=0;
REAL cpinit=0.0;
map *mapPtr = new map(hasInverse,ambDim,cpinit,evalFnp,0,derivFnp);
assert(mapPtr!=0);
// set fitData for the d2LTManifold
int nodesPerElt = (order+1)*(order+2)/2;
Matrix<int> fitData(numCoarseElts,5+nodesPerElt);
for (j=0; j<n1; j++)
for (i=0; i<n2; i++) {
fitData(i*twon1+2*j,0) = i*twon1+(2*j-1+twon1)%twon1;
fitData(i*twon1+2*j,1) = i*twon1+2*j+1;
fitData(i*twon1+2*j,2) = (i*twon1+2*j+twon1+1)%numCoarseElts;
fitData(i*twon1+2*j+1,0) = i*twon1+(2*j+2)%twon1;
fitData(i*twon1+2*j+1,1) = i*twon1+2*j;
fitData(i*twon1+2*j+1,2) = (i*twon1+2*j-twon1+numCoarseElts)%numCoarseElts;
}
for (j=0; j<n1; j++)
for (i=0; i<n2; i++) {
fitData(i*twon1+2*j,3) = -1;
fitData(i*twon1+2*j+1,3) = -1;
}
for (j=0; j<n1; j++)
for (i=0; i<n2; i++) {
fitData(i*twon1+2*j,4) = 0;
fitData(i*twon1+2*j+1,4) = 0;
}
for (j=0; j<n1; j++)
for (i=0; i<n2; i++) {
fitData(i*twon1+2*j,5) = ((n1+1)*(i+1)+j)%numNodes;
fitData(i*twon1+2*j,6) = ((n1+1)*i+j)%numNodes;
fitData(i*twon1+2*j,7) = ((n1+1)*(i+1)+(j+1)%(n2+1))%numNodes;
fitData(i*twon1+2*j+1,5) = (n1+1)*i+(j+1)%(n2+1);
fitData(i*twon1+2*j+1,6) = ((n1+1)*(i+1)+(j+1)%(n2+1))%numNodes;
fitData(i*twon1+2*j+1,7) = ((n1+1)*i+j)%numNodes;
}
// set nodeList for the d2LTManifold, and set which nodes are on which boundaries.
// (includes specification of initial tangent and normal data)
baseVector<d2GridPt *> nodeList(numNodes);
VEC point(ambDim),trans(3);
MAT normal(ambDim,ambDim-2),tangent(ambDim,2);
REAL tau,theta;
tau = 0.0;
theta = 0.0;
int numTruncDims=2;
Matrix<int> *bndryNodePtr;
bndryNodePtr = new Matrix<int>(numNodes,numTruncDims);
assert(bndryNodePtr!=0);
for (j=0; j<n1+1; j++) {
for (i=0; i<n2+1; i++) {
point[0] = 0.0;
point[1] = theta;
point[2] = tau;
trans[0] = 1.0;
trans[1] = 0.0;
trans[2] = 0.0;
normal.setColumn(0,trans);
trans[0] = 0.0;
trans[1] = 1.0;
trans[2] = 0.0;
tangent.setColumn(0,trans);
trans[0] = 0.0;
trans[1] = 0.0;
trans[2] = 1.0;
tangent.setColumn(1,trans);
nodeList[(n1+1)*i+j] = new d2GridPt(1,1,2,ambDim,point,tangent,normal);
assert(nodeList[(n1+1)*i+j] != 0);
if (j==0) (*bndryNodePtr)((n1+1)*i+j,0)=-1;
if (j==1) (*bndryNodePtr)((n1+1)*i+j,0)=-1;
if (j==n1) (*bndryNodePtr)((n1+1)*i+j,0)=1;
if (j==n1-1) (*bndryNodePtr)((n1+1)*i+j,0)=1;
if (i==0) (*bndryNodePtr)((n1+1)*i+j,1)=-1;
if (i==1) (*bndryNodePtr)((n1+1)*i+j,1)=-1;
if (i==n2) (*bndryNodePtr)((n1+1)*i+j,1)=1;
if (i==n2-1) (*bndryNodePtr)((n1+1)*i+j,1)=1;
tau+=twoPi/((REAL) n2);
}
tau=0.0;
theta+=twoPi/((REAL) n1);
}
// create the d2LTManifold
d2LTManifold *manPtr;
int irrGrid=1,wAElts=1;
manPtr = new d2LTManifold(ambDim,numCoarseElts,numTruncDims,order,numNodes,numNodes,
nodeList,fitData,irrGrid,wAElts);
// create the d2graphMan
d2graphMan *gMan1Ptr;
gMan1Ptr = new d2graphMan(manPtr);
// data structure for d2TCGTO_s::iterate method
baseVector<VEC *> *basePtPtr = new baseVector<VEC *>(numNodes);
assert(basePtPtr!=0);
for (i=0; i<numNodes; i++) {
(*basePtPtr)[i] = new VEC(3);
assert((*basePtPtr)[i]!=0);
}
// create the d2TCGTO_s
d2TCGTO_s *gtPtr;
int MAXITS=1000;
REAL TOLF=1.0e-6,TOLMIN=1.0e-8,TOLX=1.0e-9,STPMX=100.0;
REAL bndryTol=1.0e-2;
int rInitSrchData = 1,wInitSrchData = 0;
int numPatchElts = 50;
gtPtr = new d2TCGTO_s(mapPtr,gMan1Ptr,bndryNodePtr,bndryTol,numTruncDims,rInitSrchData,
"srchDatMS.dat",wInitSrchData,numPatchElts,MAXITS,TOLF,TOLMIN,TOLX,STPMX);
assert(gtPtr!=0);
// parameters for d2TCGTO_s::iterate/d2TCGTO_s::continuation
REAL projectNeigh=1.0;
REAL gtAcc=1.0e-4;
int fullOutput=0,gtMaxit=30;
// gtPtr->iterate(fullOutput,gtMaxit,gtAcc,projectNeigh,basePtPtr,"MS_A");
// parameters for d2TCGTO_s::continuation
int contMesh=5;
REAL start=0.0,end=0.5;
REAL cosTol=0.0; // for adaptive refinement, maximum slope.
// cerr << "angle tolerence, degrees: " << 180.0*acos(cosTol)/Pi << endl;
cout << gtPtr->continuation(contMesh,start,end,projectNeigh,cosTol,gtAcc,"MS_AA");
delete gtPtr;
for (i=0; i<numNodes; i++) delete (*basePtPtr)[i];
delete basePtPtr;
delete gMan1Ptr;
delete manPtr;
delete mapPtr;
return 0;
}
VEC evalFnp(VEC arg,REAL eps)
{
VEC val(3);
val[0] = beta*arg[0] + eps*sin(arg[2]);
val[1] = arg[1] + alpha*sin(arg[1]) + eps*arg[0];
val[2] = arg[2] - alpha*cos(arg[1])*sin(arg[2]) + eps*arg[0];
return val;
}
MAT derivFnp(VEC arg,REAL eps)
{
MAT val(3,3);
// first equation
val(0,0) = beta;
val(0,1) = 0.0;
val(0,2) = eps*cos(arg[2]);
// second equation
val(1,0) = eps;
val(1,1) = 1.0 + alpha*cos(arg[1]);
val(1,2) = 0.0;
// third equation
val(2,0) = eps;
val(2,1) = alpha*sin(arg[1])*sin(arg[2]);
val(2,2) = 1.0 - alpha*cos(arg[2])*cos(arg[1]);
return val;
}
written by Dimitri van Heesch, © 1997-1998