Class hierarchy   Compound list   File list   Compound Members   File Members   Related Pages   Examples  

---

Here is an example of how to use the "GTO_hyp" class to continue a saddle-type invariant manifold of a dynamical system over a parameter range. The initial invariant manifold is a saddle-type invariant closed curve of an ODE (the Lorenz system).

// In this example we continue a saddle type periodic orbit of the Lorenz system, 
// integrated using the Taylor series method.

#include "GTO_hyp.h"
#include "routines.h"
#include "svdcmp.h"
#include "ludcmp.h"
#include <fstream.h>

// parameters of the Lorenz system:
const REAL paraSigma = 10.0;
const REAL paraBeta = 8.0/3.0;
const REAL paraR = 20.0;

// Evaluate the vector field (the Lorenz system)
VEC evalOde(VEC,REAL);

// For reading in the initial data:
void readEigenVectors(char [],baseVector<MAT> &,int);
void readInitialPoints(char [],baseVector<VEC> &,VEC &,int);

// For the Taylor series method and its derivative:
const int tsmOrder=20;  // order of the Taylor series method
baseVector<VEC *> *phiPtr;
baseVector<VEC *> *xPtr;
VEC T1(7),T2(10),tsmTimestepSeq(tsmOrder);
// timestep for the Taylor series method:
const int tsmSteps=10;
const REAL tsmTimestep=0.9/((REAL) tsmSteps);
VEC tsmMapVal(VEC,REAL);
void getTaylorCoeff1(int,REAL);
MAT tsmMapDeriv(VEC,REAL);
MAT tsmMapDerivStep(MAT,REAL);
void getTaylorCoeff2(int,REAL);

main()
{
   int i,j,k;
   
// Fill in data and allocate memory for data structures for the Taylor series method
   for (i=0; i<tsmOrder; i++) tsmTimestepSeq[i] = pow(tsmTimestep,i);
   phiPtr = new baseVector<VEC *>(tsmOrder);
   assert(phiPtr!=0);
   for (i=0; i<tsmOrder; i++) {
      (*phiPtr)[i] = new VEC(3);
      assert((*phiPtr)[i]!=0);
   }
   xPtr = new baseVector<VEC *>(tsmOrder);
   assert(xPtr!=0);
   for (i=0; i<tsmOrder; i++) {
      (*xPtr)[i] = new VEC(4);
      assert((*xPtr)[i]!=0);
   }
   
// Create the dynamical system, the Lorenz model integrated by Taylor series method.
   REAL cpinit = 20.0;
   map *mapPtr = new map(0,3,cpinit,tsmMapVal,0,tsmMapDeriv);
   assert(mapPtr!=0);

// Read in points on the initial invariant curve. To obtain a value on the initial 
// invariant curve, we have iterated Newton's method on a Poincare map. (Poincare
// map method)
   const int numSteps=2000;
   VEC timeStepSequence(numSteps);
   baseVector<VEC> valueSequence(numSteps);
   readInitialPoints("points.dat",valueSequence,timeStepSequence,numSteps);
   
// Read in the initial stable/unstable splitting. This is determined by
// the eigenvectors of the Poincare map derivative, for the 
// Poincare maps defined on the orthogonal transversal to the invariant curve
// at each point of the discrete approximation. The derivative of the Poincare
// map was approximated by integrating the variational equation. 
   VEC testVec1(3),testVec2(3);
   MAT dummy(3,2);
   baseVector<MAT> invBundle(numSteps);
   for (i=0; i<numSteps; i++) invBundle[i] = dummy;
   cout << "reading eigenvectors" << endl;
   readEigenVectors("eVecs.dat",invBundle,numSteps);
   
// ensure that eigenvectors on adjacent elements are oriented the same way
   for (i=0; i<numSteps-1; i++) {
      for (j=0; j<2; j++) {
	 testVec1 = (invBundle[i]).getColumn(j);
	 testVec2 = (invBundle[i+1]).getColumn(j);
	 if (testVec1*testVec2 < 0.0) 
	    (invBundle[i+1]).setColumn(j,-1.0*testVec2);
      }
   }
      
      
      
   // ************************************************************************** //
   // ****************** Create Objects for Graph Transform Iteration ******************** //
   // ************************************************************************** //
   
   const int ambDim=3;
   const int stableDim=1;
   const int unstableDim=1;
   const int numTruncDims=0;
   const int factor=10;
   const int numElts=numSteps/factor;
   const int nodesPerElt=2;
   
   // ********************************************************************** //
   // ***************** allocate and assign gridPtr for manifold_hyp ***************** //
   
   cout << "allocating and assigning gridPtr" << endl;
   
   baseMatrix<VEC *> *gridPtr = new baseMatrix<VEC *>(numElts,nodesPerElt);
   assert(gridPtr!=0);
   
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++) {
	 (*gridPtr)(i,j) = new VEC(ambDim);
	 assert((*gridPtr)(i,j)!=0);
      }
   
   
   for (i=0; i<numElts; i++) {
      for (j=0; j<nodesPerElt-1; j++) {
	 k = i*(nodesPerElt-1) + j;
	 *((*gridPtr)(i,j)) = valueSequence[factor*k];
      }
   }
   
   for (i=0; i<numElts-1; i++)
      *((*gridPtr)(i,nodesPerElt-1)) = *((*gridPtr)(i+1,0));
   
   *((*gridPtr)(numElts-1,nodesPerElt-1)) = *((*gridPtr)(0,0));
   
   
   // for second manifold_hyp:
   
   baseMatrix<VEC *> *secondgridPtr = new baseMatrix<VEC *>(numElts,nodesPerElt);
   assert(secondgridPtr!=0);
   
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++) {
	 (*secondgridPtr)(i,j) = new VEC(ambDim);
	 assert((*secondgridPtr)(i,j)!=0);
      }
   
   
   for (i=0; i<numElts; i++) {
      for (j=0; j<nodesPerElt-1; j++) {
	 *((*secondgridPtr)(i,j)) = *((*gridPtr)(i,j));
      }
   }
   
   for (i=0; i<numElts-1; i++)
      *((*secondgridPtr)(i,nodesPerElt-1)) = *((*gridPtr)(i,nodesPerElt-1));
   
   *((*secondgridPtr)(numElts-1,nodesPerElt-1)) = *((*gridPtr)(numElts-1,nodesPerElt-1));
   
   
   
   // ****************************************************************************** //
   // ***************** allocate and assign tangentBasisPtr for manifold_hyp ***************** //
   
   cout << "allocating and assigning tangentBasisPtr" << endl;
   
   baseMatrix<MAT *> *tangentBasisPtr = new baseMatrix<MAT *>(numElts,2);
   assert(tangentBasisPtr!=0);
   
   for (i=0; i<numElts; i++)
      for (j=0; j<2; j++) {
	 (*tangentBasisPtr)(i,j) = new MAT(ambDim,1);
	 assert((*tangentBasisPtr)(i,j)!=0);
      }
   
   REAL norm;
   
   for (i=0; i<numElts; i++) {
      genUse = evalOde(valueSequence[factor*i],paraR);
      
      norm = genUse.norm();
      assert(norm>zerotol);
      genUse = (1.0/norm)*genUse;
      
      *((*tangentBasisPtr)(i,0)) = genUse;
   }
   
   
   for (i=0; i<numElts-1; i++) 
      *((*tangentBasisPtr)(i,1)) = *((*tangentBasisPtr)(i+1,0));
   
   *((*tangentBasisPtr)(numElts-1,1)) = *((*tangentBasisPtr)(0,0));
   
   
   // for second manifold_hyp:
   
   baseMatrix<MAT *> *secondtanBasisPtr = new baseMatrix<MAT *>(numElts,2);
   assert(secondtanBasisPtr!=0);
   
   for (i=0; i<numElts; i++)
      for (j=0; j<2; j++) {
	 (*secondtanBasisPtr)(i,j) = new MAT(ambDim,1);
	 assert((*secondtanBasisPtr)(i,j)!=0);
      }
   
   
   for (i=0; i<numElts; i++)
      *((*secondtanBasisPtr)(i,0)) = *((*tangentBasisPtr)(i,0));
   
   for (i=0; i<numElts-1; i++) 
      *((*secondtanBasisPtr)(i,1)) = *((*tangentBasisPtr)(i,1));
   
   *((*secondtanBasisPtr)(numElts-1,1)) = *((*tangentBasisPtr)(numElts-1,1));
   
   
   // ************************************************************************* //
   // ***************** allocate and assign stBGridPtr for manifold_hyp ***************** //
   
   cout << "allocating and assigning stBGridPtr" << endl;
   
   baseMatrix<VEC *> *stBGridPtr = new baseMatrix<VEC *>(numElts,nodesPerElt);
   assert(stBGridPtr!=0);
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++) {
	 (*stBGridPtr)(i,j) = new VEC(ambDim);
	 assert((*stBGridPtr)(i,j)!=0);
      }
   
   
   for (i=0; i<numElts; i++) {
      for (j=0; j<nodesPerElt-1; j++) {
	 k = i*(nodesPerElt-1) + j;
	 *((*stBGridPtr)(i,j)) = (invBundle[factor*k]).getColumn(1) + *((*gridPtr)(i,j));
      }
   }
   
   for (i=0; i<numElts-1; i++)
      *((*stBGridPtr)(i,nodesPerElt-1)) = *((*stBGridPtr)(i+1,0));
   
   *((*stBGridPtr)(numElts-1,nodesPerElt-1)) = *((*stBGridPtr)(0,0));
   
   
   // for second manifold_hyp:
   
   baseMatrix<VEC *> *secondstBGridPtr = new baseMatrix<VEC *>(numElts,nodesPerElt);
   assert(secondstBGridPtr!=0);
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++) {
	 (*secondstBGridPtr)(i,j) = new VEC(ambDim);
	 assert((*secondstBGridPtr)(i,j)!=0);
      }
   
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++)
	 *((*secondstBGridPtr)(i,j)) = *((*stBGridPtr)(i,j));
   
   
   // ************************************************************************* //
   // ***************** allocate and assign unBGridPtr for manifold_hyp ***************** //
   
   cout << "allocating and assigning unBGridPtr" << endl;
   
   baseMatrix<VEC *> *unBGridPtr = new baseMatrix<VEC *>(numElts,nodesPerElt);
   assert(unBGridPtr!=0);
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++) {
	 (*unBGridPtr)(i,j) = new VEC(ambDim);
	 assert((*unBGridPtr)(i,j)!=0);
      }
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt-1; j++) {
	 k = i*(nodesPerElt-1) + j;
	 *((*unBGridPtr)(i,j)) = (invBundle[factor*k]).getColumn(0) + *((*gridPtr)(i,j));
      }
   
   for (i=0; i<numElts-1; i++)
      *((*unBGridPtr)(i,nodesPerElt-1)) = *((*unBGridPtr)(i+1,0));
   
   *((*unBGridPtr)(numElts-1,nodesPerElt-1)) = *((*unBGridPtr)(0,0));
   
   
   // for second manifold_hyp:
   
   baseMatrix<VEC *> *secondunBGridPtr = new baseMatrix<VEC *>(numElts,nodesPerElt);
   assert(secondunBGridPtr!=0);
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++) {
	 (*secondunBGridPtr)(i,j) = new VEC(ambDim);
	 assert((*secondunBGridPtr)(i,j)!=0);
      }
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++)
	 *((*secondunBGridPtr)(i,j)) = *((*unBGridPtr)(i,j));
   
   
   // ************************************************************************** //
   // *************************** create manifold_hyp ****************************** //
   
   cout << "allocating and assigning manifold_hyp" << endl;
   
   manifold_hyp *saddleManPtr = new manifold_hyp(ambDim,stableDim,unstableDim,numTruncDims,numElts,
				   nodesPerElt,gridPtr,tangentBasisPtr,
				   stBGridPtr,unBGridPtr);
   assert(saddleManPtr!=0);
   
   saddleManPtr->checkGridPts(2.0E-15);
   saddleManPtr->checkTangentPts(2.0E-15);
   saddleManPtr->checkSplitPts(4.0E-15);
   saddleManPtr->checkSplit(4.0E-15);
      
   
   // second manifold_hyp for obtaining the new base manifold:
   
   manifold_hyp *wNBManPtr = new manifold_hyp(ambDim,stableDim,unstableDim,numTruncDims,numElts,
				   nodesPerElt,secondgridPtr,secondtanBasisPtr,
				   secondstBGridPtr,secondunBGridPtr);
   assert(wNBManPtr!=0);
   
   
   
   // ***************** allocate and assign grid2Ptr for graphMan_hyp ***************** //
   
   cout << "allocating and assigning grid2Ptr" << endl;
   
   baseMatrix<VEC *> *grid2Ptr = new baseMatrix<VEC *>(numElts,nodesPerElt);
   assert(grid2Ptr!=0);
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++) {
	 (*grid2Ptr)(i,j) = new VEC(ambDim);
	 assert((*grid2Ptr)(i,j)!=0);
      }
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++)
	 *((*grid2Ptr)(i,j)) = *((*gridPtr)(i,j));
   
   
   // ************************* create graphMan_hyp ************************* //
   
   cout << "allocating and assigning graphMan_hyp" << endl;
   
   graphMan_hyp *gManPtr = new graphMan_hyp(saddleManPtr,grid2Ptr);
   assert(gManPtr!=0);
   
   gManPtr->zero();
   
   // ***************** allocate and assign grid3Ptr for graphMan_hyp 2 ***************** //
   
   cout << "allocating and assigning grid3Ptr" << endl;
   
   baseMatrix<VEC *> *grid3Ptr = new baseMatrix<VEC *>(numElts,nodesPerElt);
   assert(grid3Ptr!=0);
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++) {
	 (*grid3Ptr)(i,j) = new VEC(ambDim);
	 assert((*grid3Ptr)(i,j)!=0);
      }
   
   for (i=0; i<numElts; i++)
      for (j=0; j<nodesPerElt; j++)
	 *((*grid3Ptr)(i,j)) = *((*gridPtr)(i,j));
   
   
   // ********************** create graphMan_hyp 2 ********************** //
   
   cout << "allocating and assigning graphMan_hyp 2" << endl;
   
   graphMan_hyp *gMan2Ptr = new graphMan_hyp(saddleManPtr,grid3Ptr);
   assert(gManPtr!=0);
   
   gMan2Ptr->zero();
   
   // ********************** create GTO_hyp ********************** //

   cout << "allocating and assigning GTO_hyp, tsm integrator" << endl;
   cout <<"numElts = " << numElts << endl;
      
   int bracketMesh=10;
   REAL xacc=1.0E-5,pNeigh=0.5,bndryTol=0.1;
   Matrix<int> *bndryEltPtr;

   GTO_hyp *gTPtr = new GTO_hyp(bracketMesh,xacc,pNeigh,mapPtr,wNBManPtr,
		   gManPtr,gMan2Ptr,bndryTol,bndryEltPtr);
   assert(gTPtr!=0);
      
// parameters for GTO_hyp::continuation
   REAL Start=20.0,End=10.0,gtAcc=1.0E-4,splitgtAcc=1.0E-12;
   int contMesh=500;      
      
   gTPtr->continuation(contMesh,Start,End,gtAcc,splitgtAcc,"hLrEE");
   
   return 0;
}




// *****************************************************************************************************


// Defines the dynamical system (Taylor series method applied to the Lorenz system).
VEC tsmMapVal(VEC arg,REAL eps)
{
   VEC value(3);
   REAL store;
   int i,k;
   
   for (i=0; i<tsmSteps; i++) {
      (*((*phiPtr)[0]))[0] = arg[0];
      (*((*phiPtr)[0]))[1] = arg[1];
      (*((*phiPtr)[0]))[2] = arg[2];
      
      for (k=0; k<tsmOrder-1; k++) getTaylorCoeff1(k,eps);
      
      store=0.0;
      for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*phiPtr)[k]))[0]);
      value[0] = store;
      
      store=0.0;
      for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*phiPtr)[k]))[1]);
      value[1] = store;
      
      store=0.0;
      for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*phiPtr)[k]))[2]);
      value[2] = store;
      
      arg=value;
   }
   
   return value;
}

// Defines the dynamical system (Taylor series method).
// this method assumes the 0th through (k)th taylor coefficients
// have been computed (for phi and x), k>=0
void getTaylorCoeff1(int k,REAL eps)
{
   REAL store;
   int j;
   
   T1[0] = paraSigma*((*((*phiPtr)[k]))[1]);
   T1[1] = -paraSigma*((*((*phiPtr)[k]))[0]);
   
   (*((*phiPtr)[k+1]))[0] = (T1[0] + T1[1])/((REAL) (k+1));
   
   T1[2] = eps*((*((*phiPtr)[k]))[0]);
   T1[3] = -((*((*phiPtr)[k]))[1]);
   
   store=0.0;
   for (j=0; j<=k; j++) store += ((*((*phiPtr)[j]))[0])*((*((*phiPtr)[k-j]))[2]);
   T1[4] = -store;
   
   (*((*phiPtr)[k+1]))[1] = (T1[2] + T1[3] + T1[4])/((REAL) (k+1));
   
   store=0.0;
   for (j=0; j<=k; j++) store += ((*((*phiPtr)[j]))[0])*((*((*phiPtr)[k-j]))[1]);
   T1[5] = store;
   
   T1[6] = -paraBeta*((*((*phiPtr)[k]))[2]);
   
   (*((*phiPtr)[k+1]))[2] = (T1[5] + T1[6])/((REAL) (k+1));
   
}


// Defines the derivative of the dynamical system (Taylor series method).
MAT tsmMapDeriv(VEC arg,REAL eps)
{
   MAT value(3,3);
   REAL store;
   int i,k;
   
   value(0,0)=1.0;
   value(1,0)=0.0;
   value(2,0)=0.0;
   
   value(0,1)=0.0;
   value(1,1)=1.0;
   value(2,1)=0.0;
   
   value(0,2)=0.0;
   value(1,2)=0.0;
   value(2,2)=1.0;
   
   for (i=0; i<tsmSteps; i++) {
      (*((*phiPtr)[0]))[0] = arg[0];
      (*((*phiPtr)[0]))[1] = arg[1];
      (*((*phiPtr)[0]))[2] = arg[2];
      
      for (k=0; k<tsmOrder-1; k++) getTaylorCoeff1(k,eps);
      
      value = tsmMapDerivStep(value,eps);
      
      store=0.0;
      for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*phiPtr)[k]))[0]);
      arg[0] = store;
      
      store=0.0;
      for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*phiPtr)[k]))[1]);
      arg[1] = store;
      
      store=0.0;
      for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*phiPtr)[k]))[2]);
      arg[2] = store;
   }
   
   return value;
}


// Defines the derivative of the dynamical system (Taylor series method).
MAT tsmMapDerivStep(MAT arg2,REAL eps)
{
   MAT value(3,3);
   REAL store;
   int k;
   
   // fill in (*((*xPtr)[0]))[i], i=0,1,2,3.
   
   (*((*xPtr)[0]))[0] = arg2(0,0);
   (*((*xPtr)[0]))[1] = arg2(1,0);
   (*((*xPtr)[0]))[2] = arg2(2,0);
   (*((*xPtr)[0]))[3] = 0.0;
   
   for (k=0; k<tsmOrder-1; k++) getTaylorCoeff2(k,eps);
   
   store=0.0;
   for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*xPtr)[k]))[0]);
   value(0,0) = store;
   
   store=0.0;
   for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*xPtr)[k]))[1]);
   value(1,0) = store;
   
   store=0.0;
   for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*xPtr)[k]))[2]);
   value(2,0) = store;
   
   (*((*xPtr)[0]))[0] = arg2(0,1);
   (*((*xPtr)[0]))[1] = arg2(1,1);
   (*((*xPtr)[0]))[2] = arg2(2,1);
   (*((*xPtr)[0]))[3] = 0.0;
   
   for (k=0; k<tsmOrder-1; k++) getTaylorCoeff2(k,eps);
   
   store=0.0;
   for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*xPtr)[k]))[0]);
   value(0,1) = store;
   
   store=0.0;
   for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*xPtr)[k]))[1]);
   value(1,1) = store;
   
   store=0.0;
   for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*xPtr)[k]))[2]);
   value(2,1) = store;
   
   (*((*xPtr)[0]))[0] = arg2(0,2);
   (*((*xPtr)[0]))[1] = arg2(1,2);
   (*((*xPtr)[0]))[2] = arg2(2,2);
   (*((*xPtr)[0]))[3] = 0.0;
   
   for (k=0; k<tsmOrder-1; k++) getTaylorCoeff2(k,eps);
   
   store=0.0;
   for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*xPtr)[k]))[0]);
   value(0,2) = store;
   
   store=0.0;
   for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*xPtr)[k]))[1]);
   value(1,2) = store;
   
   store=0.0;
   for (k=0; k<tsmOrder; k++) store += tsmTimestepSeq[k]*((*((*xPtr)[k]))[2]);
   value(2,2) = store;
   
   
   return value;
}


// Defines the dynamical system (Taylor series method).
// this method assumes getTaylorCoeff1(k-1) has been called.
void getTaylorCoeff2(int k,REAL eps)
{
   REAL store;
   int j;
   
   T2[0] = -paraSigma*((*((*xPtr)[k]))[0]);
   T2[1] = paraSigma*((*((*xPtr)[k]))[1]);
   
   (*((*xPtr)[k+1]))[0] = (T2[0] + T2[1])/((REAL) (k+1));
   
   T2[2] = eps*((*((*xPtr)[k]))[0]);
   
   store=0.0;
   for (j=0; j<=k; j++) store += ((*((*phiPtr)[j]))[2])*((*((*xPtr)[k-j]))[0]);
   T2[3] = -store;
   
   T2[4] = -((*((*xPtr)[k]))[1]);
   
   store=0.0;
   for (j=0; j<=k; j++) store += ((*((*phiPtr)[j]))[0])*((*((*xPtr)[k-j]))[2]);
   T2[5] = -store;
   
   (*((*xPtr)[k+1]))[1] = (T2[2] + T2[3] + T2[4] + T2[5])/((REAL) (k+1));
   
   store=0.0;
   for (j=0; j<=k; j++) store += ((*((*phiPtr)[j]))[1])*((*((*xPtr)[k-j]))[0]);
   T2[6] = store;
   
   store=0.0;
   for (j=0; j<=k; j++) store += ((*((*phiPtr)[j]))[0])*((*((*xPtr)[k-j]))[1]);
   T2[7] = store;
   
   T2[8] = -paraBeta*((*((*xPtr)[k]))[2]);
   
   (*((*xPtr)[k+1]))[2] = (T2[6] + T2[7] + T2[8])/((REAL) (k+1));
   
   (*((*xPtr)[k+1]))[3] = k>=1 ? 0 : 1;
}


// The Lorenz system vector field.
VEC evalOde(VEC arg,REAL eps)
{
   VEC value(3);
   
   value[0] = paraSigma*(arg[1] - arg[0]);
   value[1] = eps*arg[0] - arg[1] - arg[0]*arg[2];
   value[2] = arg[0]*arg[1] - paraBeta*arg[2];
   
   return value;
}


void readInitialPoints(char filename[],baseVector<VEC> &valueSeq,VEC &timeStepSeq,int size)
{
   long old;
   int i;
   ifstream fin(filename);
   
   old = fin.setf(ios::fixed,ios::floatfield);
   fin.precision(outPrecision);
   
   for (i=0; i<size; i++) {
      fin >> timeStepSeq[i];
      fin >> valueSeq[i];
   }
   
   old = fin.setf(old,ios::floatfield);
   fin.close();
}

void readEigenVectors(char filename[],baseVector<MAT> &invBundle,int size)
{
   long old;
   int i;
   ifstream fin(filename);
   
   old = fin.setf(ios::fixed,ios::floatfield);
   fin.precision(outPrecision);
   
   for (i=0; i<size; i++) {
      fin >> invBundle[i];
   }
   
   old = fin.setf(old,ios::floatfield);
   fin.close();
}

---