ODEs - variational equations

Table of Contents

• Representation of initial condition
• Long-time integration of variational equations
• Variational equations - functional approach

Representation of initial condition

The rigorous solver requires that the initial condition for variational equations is given in one of the acceptable representation. We have implemented several algorithms and representations but C1Rect2Set and C1HORect2Set proved to be most efficient in most typical applications.

These classes represent a subset of in the form of doubleton where

• are interval vectors, contain zero
• are interval matrices, with close to orthogonal
  

They represent monodromy matrix also in the form of doubleton where

• are interval matrices, contain zero
• are interval matrices, with close to orthogonal
  

One can define an instance of C1Rect2Set or C1HORect2Set by constructor call

IVector initialCondition(...);
interval initialTime = ...;
C1Rect2Set set(initialCondition,initialTime);
// or C1HORect2Set set(initialCondition,initialTime);

Note

initialTime can be skipped and then it is set to zero

Note

the above constructor sets Identity as the initial condition for variational equations

If the initial condition is not an interval vector but an affine set

we recommend to use one of the following overloaded constructors

IVector x(...);
IMatrix C(...);
IVector r0(...);
interval initialTime = ...;
 
C1Rect2Set set(x,C,r0,initialTime);
// or
C1Rect2Set set(x,r0,initialTime);

Long-time integration of variational equations

The class TimeMap (see Long-time integration) provides interface for long-time integration of ODEs with associated variational equations. One has to specify initial condition

C1Rect2Set set(...);

and call suitable operator

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

After successful integration one can extract monodromy matrix from the object set

IMatrix monodromyMatrix = (IMatrix)set;

Note

When solving ODEs with variational equations step control mechanism predicts time steps to fix requested error tolerances both for main equations and for variational equations.

Below is an example of integration of the forced pendulum equation.

See also an example of integration of the Rossler system with variational equations: Integration of ODE along with a variational equations

Complete example (from examples/odes/ITimeMapVariationalEquationsExample.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.;
  interval initTime = 4;
  C1Rect2Set set(u,initTime);
 
  // integrate
  interval finalTime = 8;
  u = timeMap(finalTime,set);
 
  // print result
  cout << "u=" << u << endl;
  cout << "monodromyMatrix=" << (IMatrix)set << endl;
  return 0;
}
 
/*output:
u={[-1.3297388770241241, -1.3297388770240786],[0.056757688397800897, 0.05675768839784965]}
monodromyMatrix={
{[0.25336614973537314, 0.25336614973554084],[1.3786621618302277, 1.3786621618305588]},
{[-1.0950745933673036, -1.0950745933671608],[-2.0118627006388627, -2.011862700638559]}
}
*/

Variational equations - functional approach

As in the case of computing trajectories (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.

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

#include <iostream>
#include "capd/capdlib.h"
 
using namespace capd;
using namespace std;
 
int main(){
  // 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.;
  interval initTime = 4;
  C1Rect2Set set(u,initTime);
 
  // define functional object
  ITimeMap::SolutionCurve solution(initTime);
 
  // and integrate
  interval finalTime = 10;
  timeMap(finalTime,set,solution);
 
  // then you can ask for any intermediate point on the trajectory
  cout.precision(17);
  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(10)=" << solution(10) << endl;
 
  cout << "monodromyMatrix(4)=" << solution.derivative(4) << endl;
  cout << "monodromyMatrix(5,5.125)=" << solution.derivative(interval(5,5.125)) << endl;
  cout << "monodromyMatrix(10)=" << solution.derivative(10) << endl;
  return 0;
}
 
/* output:
domain = [4,10]
solution(4) ={[1, 1],[2, 2]}
solution(5,5.125) ={[2.0659018710648156, 2.0826100235006182],[-0.10794601246272083, 0.12123454798834073]}
solution(10)={[1.2048370333482177, 1.2048370333482994],[1.1704148266177106, 1.1704148266177841]}
monodromyMatrix(4)={
{[1, 1],[-0, 0]},
{[-0, 0],[1, 1]}
}
monodromyMatrix(5,5.125)={
{[0.97223163741647423, 0.99012219196297169],[1.0255143836882481, 1.1731157012509985]},
{[0.11326407278286743, 0.17284240065817016],[1.1480297669989168, 1.2147644133774944]}
}
monodromyMatrix(10)={
{[-1.2196146888865109, -1.2196146888862074],[-2.1953473294485577, -2.1953473294479529]},
{[-0.024860798371007034, -0.024860798370730811],[-0.86468135954864922, -0.86468135954812031]}
}
*/