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
- define an instance of ITimeMap ITimeMap timeMap(solver);
- specify initial condition interval initTime = ...; // needed for nonautonomous systems onlyIVector initialCondition(...);C0Rect2Set set(initialCondition,initTime);// or C0HORect2Set set(initialCondition,initTime);// or C0TripletonSet set(initialCondition,initTime);// or C0HOTripletonSet set(initialCondition,initTime);
- and integrate
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;
// 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
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 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
![$ t\in[L,R]$](form_87.png)
one can compute which is guaranteed bound for the solution to the IVP at time t (it can be interval).IVector u = solution(t);- 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]}
*/
