CAPD DynSys Library  5.2.0
Enclosure for a trajectory between time steps

The sources of the following examples can be found in the examples/encloseTrajectory directory of the CAPD package.

Enclosure for a trajectory between time steps

#include "capd/capdlib.h"
using namespace capd;
using namespace std;
int main()
{
try{
cout.precision(12);
// This is the vector field for the Planar Restricted Circular 3 Body Problem
IMap vectorField("par:mu,mj;var:x,y,dx,dy;fun:dx,dy,x-mj*(x+mu)*sqrt((x+mu)^2+y^2)^(-3)-mu*(x-mj)*sqrt((x-mj)^2+y^2)^(-3)+2*dy,y*(1-mj*sqrt((x+mu)^2+y^2)^(-3)-mu*sqrt((x-mj)^2+y^2)^(-3))-2*dx;");
// mass ratio
interval mu=interval(123.0)/interval(10000.0);
interval mj=1.0-mu;
// set parameter values
vectorField.setParameter("mu",mu);
vectorField.setParameter("mj",mj);
// The solver uses high order enclosure method to verify the existence of the solution.
// The order will be set to 20.
IOdeSolver solver(vectorField,20);
ITimeMap timeMap(solver);
// This is our initial condition
IVector x(4);
x[0]=0.9928634178;
x[1]=0.0;
x[2]=0.0;
x[3]=2.129213043;
// define a doubleton representation of the interval vector x
C0Rect2Set s(x);
// Here we start to integrate. The time of integration is set to T=5.
double T=5;
timeMap.stopAfterStep(true);
interval prevTime(0.);
do
{
timeMap(T,s);
interval stepMade = solver.getStep();
cout << "\nstep made: " << stepMade;
// This is how we can extract an information
// about the trajectory between time steps.
// The type CurveType is a function defined
// on the interval [0,stepMade].
// It can be evaluated at a point (or interval).
// The curve can be also differentiated wrt to time.
// We can also extract from it the 1-st order derivatives wrt.
const IOdeSolver::SolutionCurve& curve = solver.getCurve();
interval domain = interval(0,1)*stepMade;
// Here we use a uniform grid of last time step made
// to enclose the trajectory between time steps.
// You can use your own favorite subdivision, perhaps nonuniform,
// depending on the problem you want to solve.
int grid=2;
for(int i=0;i<grid;++i)
{
interval subsetOfDomain = interval(i,i+1)*stepMade/grid;
// The above interval does not need to be a subset of domain.
// This is due to rounding to floating point numbers.
// We take the intersection with the domain.
intersection(domain,subsetOfDomain,subsetOfDomain);
// Here we evaluated curve at the interval subsetOfDomain.
// v will contain rigorous bound for the trajectory for this time interval.
IVector v = curve(subsetOfDomain);
std::cout << "\nenclosure for t=" << prevTime + subsetOfDomain << ": " << v;
std::cout << "\ndiam(enclosure): " << diam(v);
}
prevTime = timeMap.getCurrentTime();
cout << "\ncurrent time: " << prevTime << endl << endl;
}while(!timeMap.completed());
}catch(exception& e)
{
cout << "\n\nException caught!\n" << e.what() << endl << endl;
}
} // END

Enclosure for derivatives (monodromy matrix) between time steps

#include "capd/capdlib.h"
using namespace capd;
using namespace std;
int main()
{
try{
cout.precision(12);
// This is the vector field for the Planar Restricted Circular 3 Body Problem
IMap vectorField("par:mu,mj;var:x,y,dx,dy;fun:dx,dy,x-mj*(x+mu)*sqrt((x+mu)^2+y^2)^(-3)-mu*(x-mj)*sqrt((x-mj)^2+y^2)^(-3)+2*dy,y*(1-mj*sqrt((x+mu)^2+y^2)^(-3)-mu*sqrt((x-mj)^2+y^2)^(-3))-2*dx;");
// mass ratio
interval mu=interval(123.0)/interval(10000.0);
interval mj=1.0-mu;
// set parameter values
vectorField.setParameter("mu",mu);
vectorField.setParameter("mj",mj);
// The solver uses high order enclosure method to verify the existence of the solution.
// The order will be set to 20.
IOdeSolver solver(vectorField,20);
ITimeMap timeMap(solver);
// This is our initial condition
IVector x(4);
x[0]=0.9928634178;
x[1]=0.0;
x[2]=0.0;
x[3]=2.129213043;
// Define a doubleton representation of the interval vector x
// and the derivative wrt initial condition.
C1Rect2Set s(x);
// Here we start to integrate. The time of integration is set to T=1.
double T=1;
timeMap.stopAfterStep(true);
solver.setAbsoluteTolerance(1e-12);
solver.setRelativeTolerance(1e-12);
interval prevTime(0.);
IMatrix prevMatrix = IMatrix::Identity(4);
do
{
timeMap(T,s);
interval stepMade = solver.getStep();
cout << "\nstep made: " << stepMade;
// This is how we can extract an information
// about the trajectory between time steps.
// The type CurveType is a function defined
// on the interval [0,stepMade].
// It can be evaluated at a point (or interval).
// The curve can be also differentiated wrt to time.
// We can also extract from it the 1-st order derivatives wrt.
const IOdeSolver::SolutionCurve& curve = solver.getCurve();
interval domain = interval(0,1)*stepMade;
// Here we use a uniform grid of last time step made
// to enclose the trajectory between time steps.
// You can use your own favourite subdivision, perhaps nonuniform,
// depending on the problem you want to solve.
int grid=10;
for(int i=0;i<grid;++i)
{
interval subsetOfDomain = interval(i,i+1)*stepMade/grid;
// The above interval does not need to be a subset of domain.
// This is due to rounding to floating point numbers.
// We take the intersection with the domain.
intersection(domain,subsetOfDomain,subsetOfDomain);
// Here we evaluated curve at the interval subsetOfDomain.
// v will contain rigorous bound for the trajectory for this time interval.
IVector v = curve(subsetOfDomain);
std::cout << "\nenclosure for t=" << prevTime + subsetOfDomain << ": " << v;
std::cout << "\ndiam(enclosure): " << diam(v);
// Here we evaluated derivative wrt x of the curve at the interval subsetOfDomain.
// M will contain rigorous bound for the monodromy matrix when integrate with with initial conditions
// x = x(prevTime) and d/dx x(prevTime) = Id.
// To obtain d/dx x(0) with multiply M by d/dx(prevTime).
IMatrix M = curve[subsetOfDomain];//*prevMatrix;
std::cout << "\nenclosure of monodromy matrix for t=" << prevTime + subsetOfDomain << ": " << M;
std::cout << "\ndiam(enclosure of monodromy matrix): " << diam(M);
}
prevTime = timeMap.getCurrentTime();
cout << "\ncurrent time: " << prevTime << endl << endl;
// here we update prev matrix
prevMatrix = IMatrix(s);
}while(!timeMap.completed());
}catch(exception& e)
{
cout << "\n\nException caught!\n" << e.what() << endl << endl;
}
} // END
capd::diffAlgebra::Curve< capd::diffAlgebra::BasicCurve< typename MapT::MatrixType > >
capd::poincare::TimeMap
TimeMap class provides methods for transport of sets (or points) by a given flow over some time inter...
Definition: TimeMap.h:33
capd::vectalg::Vector
Definition: ColumnVector.h:177
capd::vectalg::Matrix::Identity
static Matrix Identity(size_type dim)
Definition: Matrix.hpp:109
capd::IMatrix
capd::vectalg::Matrix< capd::DInterval, CAPD_DEFAULT_DIMENSION, CAPD_DEFAULT_DIMENSION > IMatrix
Definition: typedefs.h:34
capd::intervals::diam
Interval< T_Bound, T_Rnd > diam(const Interval< T_Bound, T_Rnd > &)
upper bound for a diameter of an interval
Definition: Interval_Fun.hpp:132
capd
Definition: atom.h:31
capd::dynset::C1DoubletonSet
The C1 set is represented as doubleton: x + C*r0 + B*r;.
Definition: C1DoubletonSet.h:39
capd::dynsys::OdeSolver
Definition: OdeSolver.h:39
capd::map::Map
This class is used to represent a map .
Definition: Map.h:124
capd::dynset::C0DoubletonSet
The set is represented as doubleton: x + C*r0 + B*r; and is moved by the following method.
Definition: C0DoubletonSet.h:36
main
#define main()
Definition: krak-lib.h:382
capd::intervals::intersection
bool intersection(Interval< T_Bound, T_Rnd > A_iv1, Interval< T_Bound, T_Rnd > A_iv2, Interval< T_Bound, T_Rnd > &A_rIntersection)
Intersection of two intervals.
Definition: Interval_Fun.hpp:163
capd::intervals::Interval
Definition of template class Interval.
Definition: Interval.h:41
capd::vectalg::Matrix
Definition: ColumnVector.h:174
capd::interval
intervals::DoubleInterval interval
Definition: DoubleInterval.h:36