CAPD4-PoincareSection-Example

capd-poincare.md

Code from CAPD Author

#include <iostream>
#include "capd/capdlib.h"
using namespace std;
using namespace capd;

int main(){
  cout.precision(17);
  // two oscillators
  // vf1 will operarate for x>0
  IMap vf1("var:x,y;fun:-y,x;");

  // vf2 will operarate for x<0
  IMap vf2("var:x,y;fun:-2*y,2*x;");
  
  // The Poincare section is x=0
  
  IOdeSolver s1(vf1,20);
  IOdeSolver s2(vf2,20);
  
  
  ICoordinateSection section(2,0);
  IPoincareMap p1(s1,section);
  IPoincareMap p2(s2,section);
  
  interval returnTime;
  
  IVector x(2);  
  x[0] = 1;
  x[1] = 1;

  // we are in x>0 region. So we are using vf1
  C0Rect2Set set1(x);
  IVector y = p1(set1,returnTime);
  
  cout << "over the time: [0," << returnTime.leftBound() << "] the solution remains in region x>=0.\n";
  cout << "intersection with x=0: " << y[1] << endl;

  // set y[0]=0 because we are on the section
  y[0] = 0;
  // and integrate second system
  C0Rect2Set set2(y);
  y = p2(set2,returnTime);
  cout << "over the time: [0," << returnTime.leftBound() << "] the solution remains in region x<=0.\n";
  cout << "intersection with x=0: " << y[1] << endl;
}

Output

over the time: [0,0.78539816339744684] the solution remains in region x>=0.
intersection with x=0: [1.4142135623730929, 1.4142135623731089]
over the time: [0,1.5707963267948957] the solution remains in region x<=0.
intersection with x=0: [-1.4142135623731129, -1.4142135623730909]