Solutions to IVPs as functions

Class [L]DTimeMap 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, say solution, that represents solution to IVP and you will be able to evaluate this function at any intermediate time

DVector u = solution(t);

capd::vectalg::Vector

Definition: Vector.h:54

without additional integration of ODE.

This method is recommended if you

do not know exact time t at which you will need the solution or
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 [L]DTimeMap
double initialTime = ....;

DTimeMap::SolutionCurve solution(initialTime);

capd::diffAlgebra::SolutionCurve
Definition: SolutionCurve.h:131
specify initial condition and call your timeMap integrator
double finalTime = ...;

DVector u(...);

timeMap(finalTime,u,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 DTimeMap::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]$ you can compute
DVector u = solution(t);

which is an approximate solution to IVP at time t.
Note
An exception is thrown when the argument t in call to solution(t) is out of the domain

Complete example (from examples/odes/DTimeMapSolutionCurveExample.cpp):

#include <iostream>
#include "capd/capdlib.h"
 
using namespace capd;
using namespace std;
 
int main(){
  cout.precision(21);
 
  // define an instance of class DMap that describes the vector field
  LDMap pendulum("time:t;par:omega;var:x,dx;fun:dx,sin(omega*t)-sin(x);");
  pendulum.setParameter("omega",1.);
 
  int order=20;
  LDOdeSolver solver(pendulum,order);
  LDTimeMap timeMap(solver);
 
  // specify initial condition
  LDVector u(2);
  u[0] = 1.;
  u[1] = 2.;
  long double initTime = 4;
  long double finalTime = 10;
 
  // define functional object
  LDTimeMap::SolutionCurve solution(initTime);
 
  // and integrate
  timeMap(finalTime,u,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(7) =" << solution(7) << endl;
  cout << "solution(10)=" << solution(10) << endl;
 
  return 0;
}
 
/* output:
domain = [4,10]
solution(4) ={1,2}
solution(7) ={-0.451828248033488868443,-1.67522823683216929268}
solution(10)={1.20483703334825757551,1.17041482661774720001}
*/

First order variational equations - functional approach

As in the case of solving IVPs (see Solutions to IVPs as functions) one can obtain solution to variational equation as a functional object that can be evaluated at any intermediate time. In principle there are two differences:

We can specify initial condition for variational equations. This argument (initMatrix in the example below) is optional. If it is skipped then the Identity matrix will be chosen as an initial condition for variational equations.
We have to send an additional argument to operator that receives monodromy matrix at the end of trajectory and indicates that we intend to integrate variational equations.

Complete example (from examples/odes/DTimeMapMonodromyMatrixCurveExample.cpp):

#include <iostream>
#include "capd/capdlib.h"
 
using namespace capd;
using namespace std;
 
int main(){
  // define an instance of class DMap that describes the vector field
  LDMap pendulum("time:t;par:omega;var:x,dx;fun:dx,sin(omega*t)-sin(x);");
  pendulum.setParameter("omega",1.);
 
  int order=20;
  LDOdeSolver solver(pendulum,order);
  LDTimeMap timeMap(solver);
 
  // specify initial condition
  LDVector u(2);
  u[0] = 1.;
  u[1] = 2.;
  long double initTime = 4;
  long double finalTime = 10;
  long double data[] = {1,1,3,4};
  LDMatrix initMatrix(2,2,data); // {{1,1},{3,4}}
  LDMatrix monodromyMatrix(2,2);
 
  // define functional object
  LDTimeMap::SolutionCurve solution(initTime);
 
  // and integrate
  timeMap(finalTime,u,initMatrix,monodromyMatrix,solution);
 
  // then you can ask for any intermediate point on the trajectory
  cout.precision(21);
  cout << "domain = [" << solution.getLeftDomain() << "," << solution.getRightDomain() << "]\n";
  cout << "solution(4) =" << solution(4) << endl;
  cout << "solution(7) =" << solution(7) << endl;
  cout << "solution(10)=" << solution(10) << endl;
 
  cout << "monodromyMatrix(4)=" << solution.derivative(4) << endl;
  cout << "monodromyMatrix(7)=" << solution.derivative(7) << endl;
  cout << "monodromyMatrix(10)=" << solution.derivative(10) << endl;
  cout << monodromyMatrix << endl;
  return 0;
}
 
/* output:
domain = [4,10]
solution(4) ={1,2}
solution(7) ={-0.451828248033488868443,-1.67522823683216929268}
solution(10)={1.20483703334825757551,1.17041482661774720001}
monodromyMatrix(4)={
{1,1},
{3,4}
}
monodromyMatrix(7)={
{10.1969455920725050385,13.1893571903452595699},
{-3.06072164458348302597,-3.8608572219154975219}
}
monodromyMatrix(10)={
{-7.80565667723112838715,-10.0010040066793847611},
{-2.61890487701602830734,-3.48358623656441465785}
}
{
{-7.80565667723112838715,-10.0010040066793847611},
{-2.61890487701602830734,-3.48358623656441465785}
}
*/

Second order variational equations - functional approach

As in the case of solving IVPs (see Solutions to IVPs as functions) one can obtain solution to first and second order variational equation as a functional object that can be evaluated at any intermediate time. In principle there are two differences:

We can specify initial condition for first and second order variational equations. These arguments (initMatrix and initHessian in the example below) are optional. If both of them are skipped then the Identity matrix will be chosen for first order variational equations and zero for the second order variational equations.
We have to send two additional arguments to the operator that computes trajectory. These arguments receive monodromy matrix and normalized hassian at the end of trajectory, respectively. Moreover, they indicate that we intend to integrate second order variational equations.

Complete example (from examples/odes/DTimeMapHessianCurveExample.cpp):

#include <iostream>
#include "capd/capdlib.h"
 
using namespace capd;
using namespace std;
 
void print(const LDHessian& h, long double t){
  cout << "hessian(" << t << ")={";
  for(LDHessian::size_type i=0;i<h.imageDimension();++i)
    for(LDHessian::size_type j=0;j<h.dimension();++j)
      for(LDHessian::size_type c=j;c<h.dimension();++c)
        cout << " " << h(i,j,c);
 
  cout << " }\n";
}
 
int main(){
  // define an instance of class DMap that describes the vector field
  LDMap pendulum("time:t;par:omega;var:x,dx;fun:dx,sin(omega*t)-sin(x);");
  pendulum.setParameter("omega",1.);
 
  int order=20;
  LDC2OdeSolver solver(pendulum,order);
  LDC2TimeMap timeMap(solver);
 
  // specify initial condition
  LDVector u(2);
  u[0] = 1.;
  u[1] = 2.;
  long double initTime = 4;
  long double finalTime = 10;
  long double initDerData[] = {1,1,3,4};
  long double initHessData[] = {1,3,3,0,2,1}; // {dxdx,dxdy,dydy} = {{1,3,3},{0,2,1}}
  LDMatrix initMatrix(2,2,initDerData); // {{1,1},{3,4}}
  LDHessian initHessian(2,2,initHessData);
 
  LDMatrix D(2,2);
  LDHessian H(2);
 
  // define functional object
  LDC2TimeMap::SolutionCurve solution(initTime);
 
  // and integrate
  timeMap(finalTime,u,initMatrix,initHessian,D,H,solution);
 
  // then you can ask for any intermediate point on the trajectory
  cout.precision(21);
  cout << "domain = [" << solution.getLeftDomain() << "," << solution.getRightDomain() << "]\n";
  cout << "solution(4) =" << solution(4) << endl;
  cout << "solution(7) =" << solution(7) << endl;
  cout << "solution(10)=" << solution(10) << endl;
 
  cout << "monodromyMatrix(4)=" << solution.derivative(4) << endl;
  cout << "monodromyMatrix(7)=" << solution.derivative(7) << endl;
  cout << "monodromyMatrix(10)=" << solution.derivative(10) << endl;
  cout << D << endl;
 
  print(solution.hessian(4),4);
  print(solution.hessian(7),7);
  print(solution.hessian(10),10);
  print(H,10);
  return 0;
}
 
/* output:
domain = [4,10]
solution(4) ={1,2}
solution(7) ={-0.451828248033488868443,-1.67522823683216929268}
solution(10)={1.20483703334825757551,1.17041482661774720001}
monodromyMatrix(4)={
{1,1},
{3,4}
}
monodromyMatrix(7)={
{10.1969455920725050385,13.1893571903452595699},
{-3.06072164458348302597,-3.8608572219154975219}
}
monodromyMatrix(10)={
{-7.80565667723112838715,-10.0010040066793847611},
{-2.61890487701602830734,-3.48358623656441465785}
}
{
{-7.80565667723112838715,-10.0010040066793847611},
{-2.61890487701602830734,-3.48358623656441465785}
}
hessian(4)={ 1 3 3 0 2 1 }
hessian(7)={ 42.6352224242584004026 114.667151261346096119 73.2754844459227481224 17.2114064486085741219 42.4899350405156879182 26.8990700123135739949 }
hessian(10)={ -29.1661902479786202484 -79.9855736889665852526 -52.1961182105122759693 -14.8470974997925887675 -38.9769249314583959565 -24.2494060090330983866 }
hessian(10)={ -29.1661902479786202484 -79.9855736889665852526 -52.1961182105122759693 -14.8470974997925887675 -38.9769249314583959565 -24.2494060090330983866 }
*/

Higher order variational equations - functional approach

The class [L]DCnTimeMap can compute solutions to higher order variational equations as a functional objects that can be evaluated at any intermediate time. In principle there are two differences:

We can specify initial condition for all variational equations. This argument (initJet in the example below) is optional. If it is skipped then the Identity matrix will be chosen for first order variational equations and zero for higher order variational equations.
We have to send one additional arguments to the operator that computes trajectory. This argument receives jet of solution and indicates that we intend to integrate higher order variational equations.

Complete example (from examples/odes/DTimeMapJetCurveExample.cpp):

#include <iostream>
#include "capd/capdlib.h"
 
using namespace capd;
 
int main(){
  // maximal degree of truncated Taylor series
  const int degree = 4;
  // define an instance of class DMap that describes the vector field
  DMap pendulum("time:t;par:omega;var:x,dx;fun:dx,sin(omega*t)-sin(x);",degree);
  pendulum.setParameter("omega",1.);
 
  int order = 20;
  DCnOdeSolver solver(pendulum,order);
  DCnTimeMap timeMap(solver);
 
  double initTime = 4;
  double finalTime = 10;
 
  // specify initial condition
  DJet initJet(2,degree);
  // set initial conditions
  // main ODE
  initJet(0) = 1; initJet(1) = 2;
  // init matrix
  initJet(0,0) = 1; initJet(0,1) = 1; initJet(1,0) = 3; initJet(1,1) = 4;
  // init hessian
  initJet(0,0,0) = 1; initJet(0,0,1) = 3; initJet(0,1,1) = 3; initJet(1,0,0) = 0; initJet(1,0,1) = 2; initJet(1,1,1) = 1;
  // higher order coefficients remain zero
 
  DJet J(2,degree);
  // define functional object
  DCnTimeMap::SolutionCurve solution(initTime);
 
  // integrate
  timeMap(finalTime,initJet,J,solution);
 
  // print result
  std::cout << "solution(7)=" << solution.jet(7).toString() << std::endl;
  std::cout << "solution(10)="<< J.toString() << std::endl;
  return 0;
}
 
/* output:
 
solution(7)=
value :
   {0,0}  : {-0.451828,-1.67523}
Taylor coefficients of order 1 :
   {1,0}  : {10.1969,-3.06072}
   {0,1}  : {13.1894,-3.86086}
Taylor coefficients of order 2 :
   {2,0}  : {42.6352,17.2114}
   {1,1}  : {114.667,42.4899}
   {0,2}  : {73.2755,26.8991}
Taylor coefficients of order 3 :
   {3,0}  : {145.438,187.641}
   {2,1}  : {605.182,743.839}
   {1,2}  : {812.539,969.994}
   {0,3}  : {353.441,416.38}
Taylor coefficients of order 4 :
   {4,0}  : {389.64,1142.97}
   {3,1}  : {2220.14,6172.14}
   {2,2}  : {4601.83,12302}
   {1,3}  : {4123.34,10734.3}
   {0,4}  : {1351.71,3462}
 
solution(10)=
value :
   {0,0}  : {1.20484,1.17041}
Taylor coefficients of order 1 :
   {1,0}  : {-7.80566,-2.6189}
   {0,1}  : {-10.001,-3.48359}
Taylor coefficients of order 2 :
   {2,0}  : {-29.1662,-14.8471}
   {1,1}  : {-79.9856,-38.9769}
   {0,2}  : {-52.1961,-24.2494}
Taylor coefficients of order 3 :
   {3,0}  : {-142.07,-97.0567}
   {2,1}  : {-575.954,-389.605}
   {1,2}  : {-757.59,-513.922}
   {0,3}  : {-323.487,-223.165}
Taylor coefficients of order 4 :
   {4,0}  : {-554.512,-714.868}
   {3,1}  : {-3064.1,-3836.76}
   {2,2}  : {-6216.23,-7630.76}
   {1,3}  : {-5495.57,-6668.73}
   {0,4}  : {-1790.96,-2161.71}
*/

Table of Contents

Solutions to IVPs as functions

First order variational equations - functional approach

Second order variational equations - functional approach

Higher order variational equations - functional approach