Example

This is a slightly modified example based on capdDynSys4/examples/encloseTrajectory/encloseTrajectoryBetweenTimeSteps.cpp file. The main difference is the ODE that it computes the enclosed trajectory: "par:k_0;var:t_0, v_0;fun:1.0, (k_0*t_0)" and the time of integration, 17.

Modified lines are commented with // <----- MODIFIED.

Code

/////////////////////////////////////////////////////////////////////////////
//
/// @file encloseTrajectoryBetweenTimeSteps.cpp
///
/// @author Daniel Wilczak
//
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) CAPD group
//
// This file constitutes a part of the CAPD library,
// distributed under the terms of the GNU General Public License.
// Consult  http://capd.ii.uj.edu.pl/ for details.


#include "capd/capdlib.h"

using namespace capd;
using namespace std;

int main()
{
  try{
    cout.precision(12);
    IMap vectorField("par:k_0;var:t_0, v_0;fun:1.0, (k_0*t_0)");   // <----- MODIFIED
    interval k_0=interval(1.9, 2.1);                               // <----- MODIFIED

    // set parameter values
    vectorField.setParameter("k_0",mu);

    // The solver uses high order enclosure method to verify the existence of the solution. 
    // The order will be set to 20.
    IOdeSolver solver(vectorField,20);

    ITimeMap timeMap(solver);
    // This is our initial condition
    IVector x(2);                         // <----- MODIFIED
    x[0]=0.0;                             // <----- MODIFIED
    x[1]=0.0;                             // <----- MODIFIED

    // define a doubleton representation of the interval vector x
    C0Rect2Set s(x);

    // Here we start to integrate. The time of integration is set to T=17.

    double T=17;                          // <----- MODIFIED
    timeMap.stopAfterStep(true);
    interval prevTime(0.);
    do 
    {
      timeMap(T,s);
      interval stepMade = solver.getStep();
      cout << "\nstep made: " << stepMade;

      // This is how we can extract an information 
      // about the trajectory between time steps.
      // The type CurveType is a function defined
      // on the interval [0,stepMade].
      // It can be evaluated at a point (or interval).
      // The curve can be also differentiated wrt to time.
      // We can also extract from it the 1-st order derivatives wrt.
      const IOdeSolver::SolutionCurve& curve = solver.getCurve();
      interval domain = interval(0,1)*stepMade;
   
      // Here we use a uniform grid of last time step made
      // to enclose the trajectory between time steps.
      // You can use your own favorite subdivision, perhaps nonuniform,
      // depending on the problem you want to solve.
      int grid=2;
      for(int i=0;i<grid;++i)
      {
        interval subsetOfDomain = interval(i,i+1)*stepMade/grid;
        // The above interval does not need to be a subset of domain.
        // This is due to rounding to floating point numbers.
        // We take the intersection with the domain.
        intersection(domain,subsetOfDomain,subsetOfDomain);

        // Here we evaluated curve at the interval subsetOfDomain.
        // v will contain rigorous bound for the trajectory for this time interval.
        IVector v = curve(subsetOfDomain);
        std::cout << "\nenclosure for t=" << prevTime + subsetOfDomain << ":  " << v;      
        std::cout << "\ndiam(enclosure): " << diam(v);
      }
      prevTime = timeMap.getCurrentTime();
      cout << "\ncurrent time: " << prevTime << endl << endl;

    }while(!timeMap.completed());
  }catch(exception& e)
  {
    cout << "\n\nException caught!\n" << e.what() << endl << endl;
  }
}  // END

Output

step made: [1, 1]
enclosure for t=[0, 0.5]:  {[-0, 0.5],[-0, 0.2625]}
diam(enclosure): {[0.5, 0.5],[0.2625, 0.2625]}
enclosure for t=[0.5, 1]:  {[0.5, 1],[0.2375, 1.05]}
diam(enclosure): {[0.5, 0.5],[0.8125, 0.8125]}
current time: [17, 17]

It picks the step size [1,1]. To my understanding, it should repeat the inner loop 17 times to reach the time [17, 17]. However, the inner loop is run only once because timeMap.iscompleted() is true after the first iteration.

Environment

Mac OSX 10.10.1
gcc-4.9.2

soonhokong/capddynsys-4.0-timestep_example.md

Example

Code

Output

Environment