Computations of higher order derivatives

computations

In the following examples we show how to use CAPD to integrate ODE along with variational equations to obtain not only image of a set but also normalized derivatives up to given order r with respect to initial conditions.

In general classes that we use have ICn prefix, which means that they rigorously (I) compute and store derivatives to given order (Cn).

Integration of ODE

We integrate a harmonic oscilator for a half of period with fixed time step. So we expect fist derivative to be and the other derivatives to be equal 0.

We need to include only one file

#include "capd/capdlib.h"
  // We define a vector field
  IMap f("var:x,y;fun:-y,x;",degree);
 
  int order = 10;  // order of numerical method
  int numberOfSteps = 128;
  double timeStep = PI/numberOfSteps;
 
  // We define ICnSolver - numerical method for ODE integration
  ICnOdeSolver solver(f, order);
  solver.setStep(timeStep); // fix time step and turn off step control
 
  // We set an initial condition
  IVector v(2);
  double radius = 0.1;
  v[0] = interval(1-radius, 1+radius);
  v[1] = interval(-radius, radius);
 
  CnMultiMatrixRect2Set set(v, degree);
 
  // We move set with given dynamical system (we make numberOfSteps with given timeStep)
  for(int i=0; i<numberOfSteps; ++i){
    set.move(solver);
  }
 
  // We extract information from the CnRect2Set.
  std::cout << "\n initial condition : " << v
            << "\n time of integration : " << numberOfSteps*timeStep
            << "\n value : \n  " << IVector(set)
            << "\n C^1 derivative :\n" << IMatrix(set)
            << "\n all derivatives" << set.currentSet().toString();

The next example shows how to compute higher order derivatives with respect to initial condition with automatic step control along the trajectory.

Integration of ODE with automatic time step control (TSC)

In this example we changed the ODE and initial condition but in fact the only important difference is that instead of moving set "manually" we will use TimeMap class with automatic Time Step Control (TSC).

  ICnTimeMap timeMap(solver);

This class computes image of a given set after a given time.

  timeMap(time, set);

Full source of this example

#include <iostream>
#include "capd/capdlib.h"
 
using namespace capd;
using namespace std;
 
int main(){
  const int degree = 3;
  std::cout << "\n Michelson equation \n----------------------\n";
  // vector field
  IMap f("par:c;var:x,y,z;fun:y,z,c^2-y-0.5*x*x;",degree);
  f.setParameter("c",interval(1.0));
 
  // ICnSolver - numerical method for ODE integration
  int order = 20;                 // order of numerical method
  ICnOdeSolver solver(f, order);
 
  // initial condition
  IVector v(3);
  double sizeOfSet = 0.005;
  v[0]=0.0;    v[1]= interval(1.563-sizeOfSet, 1.563+sizeOfSet);  v[2]=0.0;
  CnRect2Set set(v, degree);
 
  ICnTimeMap timeMap(solver);
  interval time = 2;
  timeMap(time, set);
 
  cout << set.currentSet().toString();
 
  return 0;
}

How to compute derivatives of order bigger than 3?

We need to specify the maximal degre of the vector field IMap in the last argument of the contructor.

  const int degree = 5;
 
  std::cout << "\n Michelson equation \n----------------------\n";
  // vector field
  IMap f("par:c;var:x,y,z;fun:y,z,1-y-0.5*x*x;", degree);

Now we can use our types instead of built-in and the rest remains the same:

#include <iostream>
#include "capd/capdlib.h"
using namespace capd;
using namespace std;
 
int main(){
  const int degree = 5;
 
  std::cout << "\n Michelson equation \n----------------------\n";
  // vector field
  IMap f("par:c;var:x,y,z;fun:y,z,1-y-0.5*x*x;", degree);
 
  // ICnSolver - numerical method for ODE integration
  int order = 20;                 // order of numerical method
  ICnOdeSolver solver(f, order);
 
  // initial condition
  IVector v(3);
  double sizeOfSet = 0.005;
  v[0]=0.0;    v[1]= interval(1.563-sizeOfSet, 1.563+sizeOfSet);  v[2]=0.0;
  CnMultiMatrixRect2Set set(v, degree);
 
  // Time map with automatic Time Step Control (TSC)
  ICnTimeMap timeMap(solver);
  interval time = 2;
 
  timeMap(time, set);
  cout << set.currentSet().toString();
 
  return 0;
}