Long-time integration

Computing solution after a fixed time

For long-time integration we recommend to use ITimeMap. It has almost the same interface as its nonrigorous version DTimeMap - see Long-time integration.

In order to compute trajectory after a given time we have to

• define an instance of ODE solver

  IMap vectorField(...);
  int order = ...;
  IOdeSolver solver(vectorField,order);

• define an instance of ITimeMap

  ITimeMap timeMap(solver);

• specify initial condition

  interval initTime = ...; // needed for nonautonomous systems only
  IVector initialCondition(...);
  C0Rect2Set set(initialCondition,initTime);
  // or C0HORect2Set set(initialCondition,initTime);
  // or C0TripletonSet set(initialCondition,initTime);
  // or C0HOTripletonSet set(initialCondition,initTime);

• and integrate

  interval finalTime = ...;
  IVector y = timeMap(finalTime,set);

Below is an example of integration of the Planar Restricted Circular Three Body Problem.

See also an example of integration of the Rossler system: Integration of an ordinary differential equation

Complete example (from examples/odesrig/ITimeMapExample.cpp):

#include <iostream>
#include "capd/capdlib.h"
using namespace std;
using namespace capd;
using capd::autodiff::Node;
 
/*
* This is vector field of the PCR3BP
* @param in is an array of independent variables
* @param out is an array of dependent variables, i.e. out = f(in)
* @param params - parameters of the map. Here we assume that mu=params[0] is a relative mass of Jupiter
*/
void pcr3bpVectorField(Node /*t*/, Node in[], int /*dimIn*/, Node out[], int /*dimOut*/, Node params[], int /*noParams*/)
{
  // Try to factorize expression as much as possible.
  // Usually we have to define some intermediate quantities.
  Node mj = 1-params[0]; // relative mass of the Sun
  Node xMu = in[0] + params[0];
  Node xMj = in[0] - mj;
  Node xMuSquare = xMu^2; // square
  Node xMjSquare = xMj^2;
  Node ySquare = in[1]^2;
 
  // power -1.5, for rigorous computation use ONLY REPRESENTABLE CONSTANTS.
  // If exponent is not representable or it is an interval then it should be a parameter of the map.
  Node factor1 = mj*((xMuSquare+ySquare)^-1.5);
  Node factor2 = params[0]*((xMjSquare+ySquare)^-1.5);
  out[0] = in[2];
  out[1] = in[3];
  out[2] = in[0] - xMu*factor1 - xMj*factor2 + 2*in[3];
  out[3] = in[1]*(1 - factor1 - factor2) - 2*in[2];
}
 
 
int main(){
  cout.precision(17);
  int dim=4, noParam=1;
  IMap vf(pcr3bpVectorField,dim,dim,noParam);
  // set value of parameter mu, which is relative mass of Jupiter
  // 0 is index of parameter, 0.0009537 is its new value
  interval mu = 0.0009537;
  interval mj = 1.-mu;
  vf.setParameter(0, mu);
 
  // define solver and TimeMap
  IOdeSolver solver(vf,30);
  ITimeMap timeMap(solver);
 
  // initial condition - very close to L1 Lyapunov orbit for Oterma's energey value C=3.03
  interval JacobiConstant = 3.03;
  IVector L1(4);
  L1[0] = 0.9208034913207400196;
  L1[3] = sqrt(sqr(L1[0]) + 2.*(mj/abs(L1[0]+mu) + mu/abs(L1[0]-mj)) + mu*mj - JacobiConstant);
 
  // just to show two different methods: Taylor and Hermite-Obreshkov
  C0TripletonSet set1(L1);
  C0HOTripletonSet set2(L1);
 
  // integrate and print output
  cout << "T=0\n" << L1 << endl;
  // actually this is rigorous bound for the return time of this orbit to Poincare section y=0 as we can see from the output
  interval finalTime1 = 3.082119126376508;
  interval finalTime2 = 3.082119126392904;
 
  cout << "T=" << finalTime1 << "\n" << timeMap(finalTime1,set1) << endl;
  cout << "T=" << finalTime2 << "\n" << timeMap(finalTime2,set2) << endl;
 
  return 0;
}
/* output:
T=0
{[0.92080349132074002, 0.92080349132074002],[0, 0],[0, 0],[0.10444767270690263, 0.10444767270691328]}
T=[3.082119126376508, 3.082119126376508]
{[0.92080349131513572, 0.92080349132163652],[-1.8971283160953771e-12, -9.315581150390686e-14],[-1.4537252665967243e-11, 3.1594168753874883e-13],[0.10444767270599591, 0.10444767271229422]}
T=[3.0821191263929042, 3.0821191263929042]
{[0.92080349131626626, 0.92080349132050598],[1.2496678848824555e-13, 1.3098424504476883e-12],[-1.0157566056247822e-11, -4.6867189830280325e-13],[0.10444767270712088, 0.10444767271116925]}
*/

Enclosure of trajectory between time steps.

Long time integration can be stopped after each step made and then resumed. It is enough to set flag

timeMap.stopAfterStep(true);

and integrate an ODE in a while loop

do{
  timeMap(finalTime,set);
  // do something with intermediate result 
}while(!timeMap.completed());  

After each step one can extract from the solver what time step has been made and get enclosure of trajectory for this time range.

Details are shown in the example bellow (see also another example Enclosure for a trajectory between time steps).

Complete example (from examples/odesrig/ITimeMapEncloseTrajectoryExample.cpp):

#include <iostream>
#include "capd/capdlib.h"
using namespace std;
using namespace capd;
 
int main(){
  cout.precision(17);
 
  // define vector field
  IMap pendulum("time:t;par:omega;var:x,dx;fun:dx,sin(omega*t)-sin(x);");
  pendulum.setParameter("omega",1.);
 
  // define solver and TimeMap
  IOdeSolver solver(pendulum,20);
  ITimeMap timeMap(solver);
 
  IVector u(2);
  u[0] = 1.;
  u[1] = 2.;
  interval initTime = 4;
 
  C0HORect2Set set(u,initTime);
  timeMap.stopAfterStep(true);
  interval finalTime = 5.;
  do{
    timeMap(finalTime,set);
    cout << "currentTime: " << set.getCurrentTime() << endl;
    cout << "x(t) = " << (IVector)set << endl;
    cout << "enclosure over last time step: " << set.getLastEnclosure() << endl << endl;
  }while(!timeMap.completed());
 
  return 0;
}
/* output:
currentTime: [4.140625, 4.140625]
x(t) = {[1.2647227942438677, 1.2647227942438684],[1.7603531149933569, 1.7603531149933582]}
enclosure over last time step: {[0.99999999999999989, 1.2647227942438684],[1.7603531149933569, 2.0000000000000004]}
 
currentTime: [4.3046875, 4.3046875]
x(t) = {[1.5287156355947062, 1.528715635594708],[1.4545060500284652, 1.4545060500284688]}
enclosure over last time step: {[1.2647227942438675, 1.528715635594708],[1.4545060500284654, 1.7603531149933584]}
 
currentTime: [4.46875, 4.46875]
x(t) = {[1.7412935160292706, 1.741293516029274],[1.1358729403345047, 1.1358729403345105]}
enclosure over last time step: {[1.528715635594706, 1.741293516029274],[1.1358729403345049, 1.454506050028469]}
 
currentTime: [4.6484375, 4.6484375]
x(t) = {[1.9138452882401127, 1.913845288240118],[0.78520666257286698, 0.78520666257287419]}
enclosure over last time step: {[1.7412935160292704, 1.913845288240118],[0.78520666257286686, 1.1358729403345107]}
 
currentTime: [4.8359375, 4.8359375]
x(t) = {[2.0272354552591754, 2.0272354552591838],[0.42585209942088714, 0.42585209942089625]}
enclosure over last time step: {[1.9138452882401125, 2.0272354552591838],[0.42585209942088709, 0.7852066625728743]}
 
currentTime: [5, 5]
x(t) = {[2.0719023157504193, 2.0719023157504308],[0.12013839737212174, 0.12013839737213276]}
enclosure over last time step: {[2.0272354552591749, 2.0719023157504308],[0.12013839737212177, 0.4258520994208963]}
 
*/

Solutions to IVPs as functions

Class ITimeMap provides methods for computing solutions to IVPs that behave like regular functions. This mechanism is written in the spirit of functional programming. Using this method you will get a function solution that represents solution to IVP and you will be able to evaluate this function at any intermediate time

IVector u = solution(t);

without additional integration of ODE.

This method is recommended if you

• do not know exact time t at which you will need the solution
• you need solution values for many intermediate times t

In order to obtain solution curve you have to

• define an instance of SolutionCurve - this is type defined in class ITimeMap

  interval initialTime = ....;
  ITimeMap::SolutionCurve solution(initialTime); 

• specify initial condition and call your timeMap integrator

  interval finalTime = ...;
  C0Rect2Set set(...);
   
  timeMap(finalTime,set,solution);

  After the last line the object solution becomes a function that you can evaluate at given time. We will describe details in the sequel.

The type ITimeMap::SolutionCurve provides (in particular) methods for

• obtaining the domain at which this curve is defined - this should be interval [initialTime,finalTime]

  double L = solution.getLeftDomain();
  double R = solution.getRightDomain();

• evaluating the function. For one can compute

  IVector u = solution(t); 

  which is guaranteed bound for the solution to the IVP at time t (it can be interval).

  Note

  An exception is thrown when the argument t in call to solution(t) is out of the domain
  

Complete example (from examples/odesrig/ITimeMapSolutionCurveExample.cpp):

#include <iostream>
#include "capd/capdlib.h"
 
using namespace capd;
using namespace std;
 
int main(){
  cout.precision(17);
  // define an instance of class IMap that describes the vector field
  IMap pendulum("time:t;par:omega;var:x,dx;fun:dx,sin(omega*t)-sin(x);");
  pendulum.setParameter("omega",1.);
 
  int order=20;
  IOdeSolver solver(pendulum,order);
  ITimeMap timeMap(solver);
 
  // specify initial condition
  IVector u(2);
  u[0] = 1.;
  u[1] = 2.;
  double initTime = 4;
  double finalTime = 10;
  C0Rect2Set set(u,initTime);
 
  // define functional object
  ITimeMap::SolutionCurve solution(initTime);
 
  // and integrate
  timeMap(finalTime,set,solution);
 
  // then you can ask for any intermediate point on the trajectory
  cout << "domain = [" << solution.getLeftDomain() << "," << solution.getRightDomain() << "]\n";
  cout << "solution(4) =" << solution(4) << endl;
  cout << "solution(5.,5.125) = " << solution(interval(5,5.125)) << endl;
  cout << "solution(7,8) =" << solution(interval(7,8)) << endl;
  cout << "solution(10)=" << solution(10) << endl;
 
  return 0;
}
 
/* output:
domain = [4,10]
solution(4) ={[1, 1],[2, 2]}
solution(5.,5.125) = {[2.061890364489789, 2.0922917471924896],[-0.10693311935888665, 0.12198962338481042]}
solution(7,8) ={[-1.3710853587244543, -0.44825760563602185],[-1.6756032292973486, 0.056820499411927523]}
solution(10)={[1.2048370333482143, 1.2048370333483036],[1.1704148266177143, 1.170414826617781]}
*/