dynIbex.cpp

#include <iostream>
#include "ibex.h"
#include <list>
#include <math.h>       /* fabs */

using namespace std;
using namespace ibex;

#define __PREC__ 1e-11
#define __METH__ RK4
#define __DURATION__ 10.0

bool isEqual (IntervalVector first, IntervalVector second);
bool isPossible(simulation simu);
bool isAccepted(simulation simu);

int main(){
//Declarer les variables
  Variable y(2);
	
	Interval ref(10);
	IntervalVector K(2);
	K[0] = Interval(0.0,1000.0);
	K[1] = Interval(0.0,100.0);
	Interval m(990.0,1010.0);
	Interval b(50.0);
	Interval c(0.4);

	// Initial conditions
	IntervalVector yinit(2);
	yinit[0] = Interval(0.,0.);
	yinit[1] = Interval(0.,0.);

  std::list<IntervalVector> stack,stack_sol;
  stack.push_back(K);
  double tol = 1e-10;

  while(!stack.empty()){
		IntervalVector K_current = stack.front();
		stack.pop_front();
  	Function ydot(y,Return(
	 		(K_current[0]*(y[0]-ref)+K_current[1]*y[1]-b*y[0]-c*y[0]*y[0])/m,
			y[0]-10)
		);
		// Ivp contruction (initial time is 0.0)
		ivp_ode problem = ivp_ode (ydot, 0.0, yinit);
		simulation simu = simulation (&problem, __DURATION__, __METH__, __PREC__);
		simu.run_simulation();
		//cout<<simu.get(10.)<<endl;
    //Si cette interval contient une interval acceptable
		if(isPossible(simu)){
    //S'il est une interval acceptable
      			cout<<"This interval is possible"<<endl;
      			if(isAccepted(simu)){
				cout<<"This interval is accepted"<<endl;
  			  	stack_sol.push_back(K_current);
		  	}
			else{
			    LargestFirst bbb(0);
			    std::pair<IntervalVector,IntervalVector> p = bbb.bisect(K_current);
			    stack.push_back(p.first);
			    stack.push_back(p.second);
			}
		}
  	}

	std::list<IntervalVector> result;
  	std::list<IntervalVector>::const_iterator it = stack_sol.begin();
  	std::list<IntervalVector>::const_iterator resIt;
    	while(it != stack_sol.end()){
  		resIt = result.begin();
			bool isNew = true;
			while(resIt != result.end()){
				if(isEqual(*resIt,*it)){
					isNew = false;
				}
				resIt++;
			}
			if(isNew){
				result.push_back(*it);
				std::cout<<*it<<endl;
			}
			it++;
		}

  	return 0;
}


// a binary predicate implemented as a function:
bool isEqual (IntervalVector first, IntervalVector second)
{ 
  bool isEq = true;
  for(int i = 0; i <5; i++){
		double _1u = first[i].ub();
		double _2u = second[i].ub();
		double _1l = first[i].lb();
		double _2l = second[i].lb();
		double detau = fabs(_1u-_2u);
		double detal = fabs(_1l-_2l);
    if(!((detau<1e-1)&&(detal<1e-1))){
      isEq=false;
    }else{
		}
  }
  return isEq; 
}

bool isPossible(simulation simu){
	return true;
}
bool isAccepted(simulation simu){
	return true;
}