Permalink
Browse files

Added the capability of finding if the orbit is periodic or not, don'…

…t know if it works correctly
  • Loading branch information...
1 parent c41697a commit 33de1d4fca4561b4b28fb40517fc34a28662a8ea muzgash committed Jun 13, 2011
Showing with 78 additions and 22 deletions.
  1. +3 −1 include/rossler.h
  2. +63 −12 src/rossler.cpp
  3. +12 −9 test/find.cpp
View
@@ -11,7 +11,9 @@ class rossler{
double sstep, double ttf);
vector<double>* pyragas();
- int isPeriodic(double a[4][30], double b[4][30]);
+ bool attractor();
+ int areEqual(double a[4][30], double b[4][30]);
+ int isPeriodic;
//mutators
void setParameters(double *params);
View
@@ -20,13 +20,15 @@ int equations(double t, const double y[], double f[],
double F=par[5];//0.0;
f[0] = -y[1]-y[2];
f[1] = y[0]+a*y[1]+y[3];
- f[2] = b+y[0]*y[2]+C*(F-y[2]);
- f[3] = -c*y[2]+d*y[3];
+ f[2] = b+y[0]*y[2];
+ f[3] = -c*y[2]+d*y[3]+C*(F-y[1]);
return GSL_SUCCESS;
}
rossler::rossler(double aa, double bb,double cc, double dd,
double sstep, double ttf){
+
+ isPeriodic=2;
timeSeries=new vector<double>;
timeSeries2=new vector<double>;
a=aa;
@@ -42,7 +44,7 @@ rossler::rossler(double aa, double bb,double cc, double dd,
pos[2]=0.0;//0.1;
pos[3]=10.0;
err=new double[4];
- transientTime=100;
+ transientTime=200;
period=2.85;
sampleRating=finalTime/step;
}
@@ -90,7 +92,7 @@ void rossler::setStep(double s){
vector<double>* rossler::pyragas(){
vector<double> feedback;//2.456;
- double params[6]={a,b,c,d,C,0.0};//crear atributo de clase
+ double params[6]={a,b,c,d,0.0,0.0};//crear atributo de clase
const gsl_odeiv_step_type *T=gsl_odeiv_step_rk4;
gsl_odeiv_step *WT=gsl_odeiv_step_alloc(T,4);
@@ -106,12 +108,15 @@ vector<double>* rossler::pyragas(){
//filling of the feedback vector
for(t=transientTime;t<transientTime+period;t+=step){
if(gsl_odeiv_step_apply(WT,t,step,pos,err,NULL,NULL,&sys)!=GSL_SUCCESS) break;
- feedback.push_back(pos[2]);
+ feedback.push_back(pos[3]);
+// cout<<pos[0]<<" "<<pos[1]<<endl;
}
+ params[4]=C; //se agrega el término de control diferente de cero
// cout<<"t= "<<t<<" period="<<20*period<<" FB size="<<feedback.size()<<endl;
// cout.flush();
//now the control comes
- double periodicity0[30],periodicity1[30]; //variables to know if there's periodicity
+ double periodicity0[4][30],periodicity1[4][30]; //variables to know if there's periodicity
+ int o=0,p=0;//para avanzar en los arreglos de periodicidad
for(i=0;i<(int)((finalTime-transientTime-period)/period);i++){
int n=0;
while(t<transientTime+(i+1)*period){
@@ -120,19 +125,43 @@ vector<double>* rossler::pyragas(){
if(gsl_odeiv_step_apply(WT,t,step,pos,err,NULL,NULL,&sys)!=GSL_SUCCESS) break;
if(t>transientTime+200*period){
- timeSeries->push_back(pos[2]);
+ timeSeries->push_back(pos[0]);
timeSeries2->push_back(pos[1]);
- cout<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<endl;
-// cout.flush();
+ //relleno los arreglos a comparar para determinar periodicidad
+// cout<<t<<" "<<transientTime+200*period<<" "<<transientTime+202*period<<endl;
+ if(o<30 && t>transientTime+201*period && t<transientTime+203*period){
+ periodicity0[0][o]=pos[0];
+ periodicity0[1][o]=pos[1];
+ periodicity0[2][o]=pos[2];
+ periodicity0[3][o]=pos[3];
+ o++;
+ cout<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<" "<<pos[3]<<endl;
+ cout.flush();
+ }
+// cout<<endl;
+ if(p<30 && t<transientTime+205*period && t>transientTime+203*period){
+ periodicity1[0][p]=pos[0];
+ periodicity1[1][p]=pos[1];
+ periodicity1[2][p]=pos[2];
+ periodicity1[3][p]=pos[3];
+ p++;
+ cout<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<" "<<pos[3]<<endl;
+ cout.flush();
+ }
+// cout<<endl;
+ if(areEqual(periodicity0,periodicity1)) isPeriodic=0;
+ else isPeriodic=1;
}
+// else cout<<pos[0]<<" "<<pos[1]<<endl;
// cout<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<endl;
// cout<<pos[1];
- feedback[n]=pos[2];
+ feedback[n]=pos[3];
// cout<<n;
// cout.flush();
n++;
t=t+step;
}
+ if(isPeriodic==false) break;
}
@@ -141,7 +170,28 @@ vector<double>* rossler::pyragas(){
return timeSeries;
}
-int rossler::isPeriodic(double a[4][30], double b[4][30]){
+bool rossler::attractor(){
+ double params[6]={a,b,c,d,0.0,0.0};//crear atributo de clase
+
+ const gsl_odeiv_step_type *T=gsl_odeiv_step_rk4;
+ gsl_odeiv_step *WT=gsl_odeiv_step_alloc(T,4);
+ gsl_odeiv_system sys={equations,NULL,4,&params};
+
+ double t=0.0;
+
+ //transient
+ for(t=0.0;t<transientTime;t+=step){
+ if(gsl_odeiv_step_apply(WT,t,step,pos,err,NULL,NULL,&sys)!=GSL_SUCCESS) break;
+ }
+ //filling of the feedback vector
+ for(t=transientTime;t<finalTime;t+=step){
+ if(gsl_odeiv_step_apply(WT,t,step,pos,err,NULL,NULL,&sys)!=GSL_SUCCESS) break;
+// cout<<pos[0]<<" "<<pos[3]<<endl;
+ }
+ return true;
+}
+
+int rossler::areEqual(double a[4][30], double b[4][30]){
int count=0;
double tol=0.02;
for(int i=0;i<30;i++){
@@ -151,6 +201,7 @@ int rossler::isPeriodic(double a[4][30], double b[4][30]){
}
if(mod<=tol) count++;
}
- if(count>25) return true;
+
+ if(count>10) {cout<<count;return true;}
else return false;
}
View
@@ -9,41 +9,44 @@ using namespace std;
int main(int argc, char **argv){
ofstream tSData;
ofstream atractor;
- double p[6]={0.5,1.3,2.25,2.3,2.85,2.9};
- for(int i=0;i<6;i++){
+// double p[6]={0.5,1.3,2.25,2.3,2.85,2.9};
+// for(double i=0.4;i<8.5;i+=0.01){
rossler cc(0.25,3.0,0.5,0.05,
- 0.001,500.0);
+ 0.01,1000.0);
vector<double> *y,*y2;
y=new vector<double>;
y2=new vector<double>;
// for(double c=0.01;c<=1;c+=0.001){
- cc.setControlWeight(0.4);
+ cc.setControlWeight(0.3);
- cc.setPeriod(p[i]);
+ cc.setPeriod(2.8);
y=cc.pyragas();
+ if(cc.isPeriodic==0){
y2=cc.get_ts2();
stringstream name;
stringstream name2d;
// name=new stringstream(stringstream::in | stringstream::out);
- name<<"../etc/timeSeries_"<<p[i]<<".dat";
+ name<<"../etc/timeSeries_"<<2.8<<".dat";
tSData.open(name.str());
- name2d<<"../etc/timeSeries_"<<p[i]<<"_2d.dat";
+ name2d<<"../etc/timeSeries_"<<2.8<<"_2d.dat";
atractor.open(name2d.str());
-
+//
for(int j=0;j<y->size();j++){
atractor<<y->at(j)<<" "<<y2->at(j)<<endl;
tSData<<y->at(j)<<endl;
}
tSData.close();
atractor.close();
+// cc.attractor();
}
+// }
// if(argv[1]=="--plot"){
// cout<<"Plotting"<<endl;
// system("python ../etc/plot.py");
// }
-// return 0;
+ return 0;
}

0 comments on commit 33de1d4

Please sign in to comment.