CAPD::DynSys Library  6.0.0
ODEs - variational equations

Table of Contents

In this section we will give an overview on simultaneous integration of ODEs and associated first order variational equations.

One step Taylor method

Class [L]DOdeSolver can also integrate first order variational equations of ODEs. Creating an instance of solver is exactly the same as in the case of computing the trajectory only (see section One-step Taylor method).

In order to compute derivatives with respect to initial conditions one has to

  • specify initial conditions for the trajectory and for variational equations
    double t = ...;
    DVector initPoint(...);
    DMatrix initMatrix(...)
    capd::vectalg::Matrix
    Definition: Matrix.h:65
    capd::vectalg::Vector
    Definition: Vector.h:54
  • define matrix that will store solution to variational equations (monodromy matrix)
    DMatrix D(...);
  • call operator that simultaneously computes solution to IVP for ODE and first order variational equations
    DVector y = solver(t,initPoint,initMatrix,D);

After the last instruction the vector y will contain approximate solution after specified time step and the matrix D will contain solution to the variational equations at the same time.

Note
Two arguments initMatrix and D can be the same object instance. This is used to obtain monodromy matrix after many steps of integration.
int dimension = ...;
DMatrix D = DMatrix::Identity(dimension);
for(int i=0;i<numberOfSteps;++i)
u = solver(t,u,D,D);
cout << "t=" << t << ", u=" << u << endl;
cout << "Monodromy matrix: " << D << endl;
}
capd::vectalg::Matrix::Identity
static Matrix Identity(size_type dim)
Definition: Matrix.hpp:109

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

#include <cstdio>
#include "capd/capdlib.h"
using namespace capd;
int main(){
// Define an instance of class DMap that describes the vector field
DMap pendulum("time:t;par:omega;var:x,y;fun:y,sin(omega*t)-sin(x);");
pendulum.setParameter("omega",1.);
// define an instance of ODE solver
int order = 20;
DOdeSolver solver(pendulum,order);
// We set a fixed time step - this disables automatic step control.
// If one really wants to use fixed time step, it is recommended to put a floating point number with many tailing zeros in the mantissa.
solver.setStep(0.125);
// specify initial condition
DVector u(2);
u[0] = 1.;
u[1] = 2.;
double t = 4;
DMatrix D = DMatrix::Identity(2);
// use one step Taylor method to integrate
do{
u = solver(t,u,D,D);
printf("t=%6f, x=%6f, y=%6f, du1/dx=%6f, du1/dy=%f6, du2/dx=%6f, du2/dy=%6f\n",t,u[0],u[1],D[0][0],D[0][1],D[1][0],D[1][1]);
}while(t<5);
return 0;
}
/* output:
t=4.125000, x=1.236999, y=1.788258, du1/dx=0.996334, du1/dy=0.1248596, du2/dx=-0.054193, du2/dy=0.996888
t=4.250000, x=1.446319, y=1.558578, du1/dx=0.987553, du1/dy=0.2491456, du2/dx=-0.082043, du2/dy=0.991905
t=4.375000, x=1.626223, y=1.318748, du1/dx=0.976821, du1/dy=0.3730036, du2/dx=-0.085989, du2/dy=0.990894
t=4.500000, x=1.775834, y=1.074735, du1/dx=0.966890, du1/dy=0.4972216, du2/dx=-0.069940, du2/dy=0.998278
t=4.625000, x=1.894909, y=0.830721, du1/dx=0.959990, du1/dy=0.6230426, du2/dx=-0.038190, du2/dy=1.016892
t=4.750000, x=1.983628, y=0.589383, du1/dx=0.957828, du1/dy=0.7519716, du2/dx=0.005228, du2/dy=1.048134
t=4.875000, x=2.042431, y=0.352264, du1/dx=0.961636, du1/dy=0.8856116, du2/dx=0.056802, du2/dy=1.092206
t=5.000000, x=2.071902, y=0.120138, du1/dx=0.972242, du1/dy=1.0255256, du2/dx=0.113523, du2/dy=1.148294
}
*/
capd::map::Map
This class is used to represent a map .
Definition: Map.h:125
main
int main(int argc, char *argv[])
Definition: argdemo.cpp:81
order
int order
Definition: tayltst.cpp:31
capd
Definition: ApplicationDesc.h:23
capd::DOdeSolver
capd::dynsys::BasicOdeSolver< capd::DMap > DOdeSolver
Definition: typedefs.h:26

First order variational equations

The class [L]DTimeMap (see Long-time integration) provides interface for long-time integration of ODEs with associated variational equations (i.e. derivatives of solution with respect to initial condition). One has to specify initial conditions for both main equations and for variational equations

int dimension = ...;
double initTime = ...;
DVector x(...);
DMatrix initMatrix(...);

and call suitable operator

double finalTime = ...;
DMatrix monodromyMatrix(dimension,dimension);
DVector y = timeMap(finalTime,x,initMatrix,monodromyMatrix,initTime);
Note
When solving IVPs with variational equations the built in step control mechanism predicts time steps to fix requested error tolerances both for main equations and for variational equations.
Note
The argument initMatrix can be skipped. Then the initial condition for variational equations is set to identity.
DVector y = timeMap(finalTime,x,monodromyMatrix,initTime);
Note
The last argument initTime can be skipped. Then it is set to zero by default.
DVector y = timeMap(finalTime,x,monodromyMatrix);

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

#include <cstdio>
#include "capd/capdlib.h"
using namespace capd;
int main(){
// 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);");
pendulum.setParameter("omega",1.);
int order = 20;
DOdeSolver solver(pendulum,order);
DTimeMap timeMap(solver);
// specify initial condition
DVector u(2);
u[0] = 1.;
u[1] = 2.;
double initTime = 4;
double finalTime = 8;
double data[] = {1,1,3,4}; // {{1,1},{3,4}}
DMatrix initMatrix(2,2,data);
DMatrix D(2,2);
// integrate
u = timeMap(finalTime,u,initMatrix,D,initTime);
// print result
printf("t=%6f, x=%6f, y=%6f, du1/dx=%6f, du1/dy=%f6, du2/dx=%6f, du2/dy=%6f\n",finalTime,u[0],u[1],D[0][0],D[0][1],D[1][0],D[1][1]);
return 0;
}
// output: t=8.000000, x=-1.329739, y=0.056758, du1/dx=4.389353, du1/dy=5.7680156, du2/dx=-7.130663, du2/dy=-9.142525
capd::poincare::TimeMap
TimeMap class provides methods for transport of sets (or points) by a given flow over some time inter...
Definition: TimeMap.h:34

Second order variational equations

The class [L]DC2TimeMap (see Long-time integration) provides interface for long-time integration of ODEs with associated first and second order variational equations. One has to specify initial conditions for main ODE and for first and second order variational equations.

int dimension = ...;
double initTime = ...;
DVector x(...);
DMatrix initMatrix(...);
DHessian initHessian(...);
capd::DHessian
capd::diffAlgebra::Hessian< capd::DVector::ScalarType, capd::DVector::csDim, capd::DMatrix::ColumnVectorType::csDim > DHessian
Definition: fdlib.h:40

and call suitable operator

double finalTime = ...;
DMatrix D(dimension,dimension);
DHessian H(dimension,dimension);
DVector y = timeMap(finalTime,x,initMatrix,initHessian,D,H,initTime);
Note
When solving IVPs with variational equations the built in step control mechanism predicts time steps to fix requested error tolerances both for main equations and for first and second order variational equations.
Note
The arguments initMatrix and initHessian can be skipped. Then the initial condition for first order variational equations is set to identity and for the second order variational equations is set to zero.
DVector y = timeMap(finalTime,x,D,H,initTime);
Note
The last argument initTime can be skipped. Then it is set to zero by default.
DVector y = timeMap(finalTime,x,D,H);

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

#include <cstdio>
#include "capd/capdlib.h"
using namespace capd;
int main(){
// 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);");
pendulum.setParameter("omega",1.);
int order = 20;
DC2OdeSolver solver(pendulum,order);
DC2TimeMap timeMap(solver);
// specify initial condition
DVector u(2);
u[0] = 1.;
u[1] = 2.;
double initTime = 4;
double finalTime = 8;
double initDerData[] = {1,1,3,4}; // {{1,1},{3,4}}
double initHessData[] = {1,3,3,0,2,1}; // {dxdx,dxdy,dydy} = {{1,3,3},{0,2,1}}
DMatrix initMatrix(2,2,initDerData);
DHessian initHessian(2,2,initHessData);
DMatrix D(2,2);
DHessian H(2,2);
// integrate
u = timeMap(finalTime,u,initMatrix,initHessian,D,H,initTime);
// print result
printf("t=%6f, x=%6f, y=%6f\n",finalTime,u[0],u[1]);
// print derivative
printf("du1/dx=%6f, du1/dy=%f6, du2/dx=%6f, du2/dy=%6f\n",D[0][0],D[0][1],D[1][0],D[1][1]);
// print normalized hessian
printf(
"d2u1/dxdx=%6f, d2u1/dxdy=%f6, d2u1/dydy=%f6, d2u1/dxdx=%6f, d2u1/dxdy=%f6, d2u1/dydy=%f6\n",
H(0,0,0),H(0,0,1),H(0,1,1),H(1,0,0),H(1,0,1),H(1,1,1)
);
return 0;
}
/* output:
t=8.000000, x=-1.329739, y=0.056758
du1/dx=4.389353, du1/dy=5.7680156, du2/dx=-7.130663, du2/dy=-9.142525
d2u1/dxdx=33.076326, d2u1/dxdy=86.6810126, d2u1/dydy=54.8328166, d2u1/dxdx=-27.228978, d2u1/dxdy=-74.3629156, d2u1/dydy=-48.3478676
*/
capd::DC2OdeSolver
capd::dynsys::BasicC2OdeSolver< capd::DMap > DC2OdeSolver
Definition: typedefs.h:27

Higher order variational equations

The class [L]DCnTimeMap (see Long-time integration) provides interface for long-time integration of ODEs with associated first, second and higher order variational equations. One has to specify initial conditions for main ODE and for all variational equations.

int dimension = ...;
double initTime = ...;
DVector x(...);
DJet initJet(...);
capd::DJet
capd::diffAlgebra::Jet< capd::DMatrix, 0 > DJet
Definition: fdlib.h:44

and call one of overloaded operators

double finalTime = ...;
DJet J(dimension,degree);
DVector y1 = timeMap(finalTime,x,J,initTime);
// or
DVector y2 = timeMap(finalTime,initJet,J,initTime);
Note
When solving IVPs with variational equations the built in step control mechanism predicts time steps to fix requested error tolerances both for main equations and for all variational equations up to requested degree.
Note
After call to the first version of operator
DVector y1 = timeMap(finalTime,x,J,initTime);
the initial condition for first order variational equations is set to identity. Higher order varational equations have initial conditions set to zero.
Note
After call to the second version of operator
DVector y1 = timeMap(finalTime,initJet,J,initTime);
the initial condition for the main equations and all variational equations are copied from the object initJet.
Note
In both cases the last argument initTime can be skipped. Then it is set to zero by default.
DVector y1 = timeMap(finalTime,x,J);
DVector y2 = timeMap(finalTime,initJet,J);

Complete example (from examples/odes/DTimeMapJetExample.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 = 8;
// 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);
// integrate
timeMap(finalTime,initJet,J,initTime);
// print result
std::cout << J.toString() << std::endl;
return 0;
}
/* output:
value :
{0,0} : {-1.32974,0.0567577}
Taylor coefficients of order 1 :
{1,0} : {4.38935,-7.13066}
{0,1} : {5.76801,-9.14253}
Taylor coefficients of order 2 :
{2,0} : {33.0763,-27.229}
{1,1} : {86.681,-74.3629}
{0,2} : {54.8328,-48.3479}
Taylor coefficients of order 3 :
{3,0} : {181.037,-106.984}
{2,1} : {735.217,-438.7}
{1,2} : {975.287,-580.696}
{0,3} : {423.543,-248.531}
Taylor coefficients of order 4 :
{4,0} : {1035.45,-207.342}
{3,1} : {5608.26,-1220.75}
{2,2} : {11214.2,-2581.6}
{1,3} : {9815.9,-2339.67}
{0,4} : {3175.13,-771.844}
*/
capd::DCnOdeSolver
capd::dynsys::BasicCnOdeSolver< capd::DMap > DCnOdeSolver
Definition: typedefs.h:28