Example Documentation

Here is an example of how to use the "d2GTO_s" class to compute an attracting invariant manifold of a dynamical system for different parameters. The initial invariant manifold is a stable invariant torus of a vector field (the "new" Lorenz model).


// In this example we continue an attracting invariant torus of a vector field,
// the "new" Lorenz model.


#include "d2GTO_s.h"
#include "cylindricalMap.h"
#include <math.h>
#include "fstream.h"
#include "d1graphP.h"

// parameters of the dynamical system
const REAL a = 0.25;
const REAL b = 4.0;
const REAL F = 1.8;
const REAL G = 1.65;

// discretization parameters, indicating the size of the grid
const int n1=64,n2=64;
const int twon1=2*n1;
const int ambDim=3;                          // ambient dimension
const int order=1;                              // order of approximation
const int numNodes=n1*n2;                // total number of grid points
const int numCoarseElts=2*n1*n2;    // number of elements in initial discretization
const REAL twoPi=2.0*Pi;


// the vector field and its derivative
VEC evalFnp(VEC,REAL);
MAT derivFnp(VEC,REAL);


main()
{
   int i,j;
      
   
// create the dynamical system, a timeDisc:
   int hasInverse=0,iterates=1,numDeriv=0;
   REAL timeStep=0.05,cpinit=0.0;
   REAL accuracy = 1.0e-8;
   REAL derivAcc = 1.0e-4;
   REAL dfrh=0.07;
   timeDisc *timeDPtr = new timeDisc(hasInverse,numDeriv,ambDim,iterates,timeStep,dfrh,
				     accuracy,derivAcc,cpinit,"rk4",evalFnp,derivFnp);
   assert(timeDPtr!=0);
   

// create the d2LTManifold
   d2LTManifold *manPtr;
   int refNumElts=40960,refNumNodes=16384;
   int hoNumNodes=147456,hoOrder=3;
   int numTDs=0;
   manPtr = new d2LTManifold(ambDim,refNumElts-numCoarseElts,numTDs,hoOrder,hoNumNodes,refNumNodes);
   cerr << "reading baseMan data" << endl;
   manPtr->readFromFile("torusOrder3_1.dat");
   
      
// create the d2graphMan
   d2graphMan *gMan1Ptr;
   cerr << "creating graphman 1" << endl;
   gMan1Ptr = new d2graphMan(manPtr);
   assert(gMan1Ptr!=0);
   
   
// create the d2GTO_s
   d2GTO_s *gtPtr;
   int numBigPatchElts=96;
   int numPatchElts=70;
   int MAXITS=1000;
   int rInitSrchData=1,useInitEltPtr=1;
   REAL TOLF=1.0e-8,TOLMIN=1.0e-10,TOLX=1.0e-11,STPMX=1.0;
   
   cerr << "creating graph transform" << endl;
   gtPtr = new d2GTO_s(timeDPtr,gMan1Ptr,rInitSrchData,useInitEltPtr,"srch128NLMrk4_2.dat",
		       numPatchElts,numBigPatchElts,MAXITS,TOLF,TOLMIN,TOLX,STPMX);
   assert(gtPtr!=0);
   
   
// parameters for d2GTO_s::iterate
   REAL projectNeigh=0.006;
   REAL gtAcc = 1.0e-4;
   int fullOutput=0,gtMaxit=99;

// data structure to pass to d2GTO_s::iterate
   cerr << "creating basePtPtr" << endl;
   baseVector<VEC *> *basePtPtr = new baseVector<VEC *>(hoNumNodes);
   assert(basePtPtr!=0);
   for (i=0; i<hoNumNodes; i++) {
      (*basePtPtr)[i] = new VEC(3);
      assert((*basePtPtr)[i]!=0);
   }

   
   cerr << "beginning iteration" << endl;
   gtPtr->iterate(fullOutput,gtMaxit,gtAcc,projectNeigh,basePtPtr,"nlmHO","dummy.dat");
   
   return 0;
}




// new Lorenz model, arg = (x,y,z).
VEC evalFnp(VEC arg,REAL eps)
{
   VEC val(3);
   
   val[0] = -arg[1]*arg[1] - arg[2]*arg[2] - a*arg[0] + a*F;
   val[1] = arg[0]*arg[1] - b*arg[0]*arg[2] - arg[1] + G;
   val[2] = b*arg[0]*arg[1] + arg[0]*arg[2] - arg[2];
   
   val = -1.0*val;        // to make the torus stable
   
   return val;
}

MAT derivFnp(VEC arg,REAL eps)
{
   MAT val(3,3);
   
   // first equation
   val(0,0) = -a;
   val(0,1) = -2.0*arg[1];
   val(0,2) = -2.0*arg[2];
   
   // second equation
   val(1,0) = arg[1] - b*arg[2];
   val(1,1) = arg[0] - 1.0;
   val(1,2) = -b*arg[0];
   
   // third equation
   val(2,0) = b*arg[1] + arg[2];
   val(2,1) = b*arg[0];
   val(2,2) = arg[0] - 1.0;
   
   val = ((REAL) -1.0)*val;
   
   return val;
}