// 2-parameter 分岐図 #include "trs5-inc.mm" #define CHOTTO 0.0000001 void dat_output(); void rem_dat_output_16(); List bifudiag_output(); Real outvf(), secf(); Integer is_identical(); Real distp(); List find_unchecked_point(); void see_point(); void output_chaotic_region_to_matfile(); List prmvalue_for(); Matrix flag, statem; Real t_fin_bifu; Func void main() { Matrix x; Real t; Integer USE_GRAPHICS; Integer px, py; Integer meshsize_x, meshsize_y; settimer(); USE_GRAPHICS= no; h=0.0; t=0.0; t_fin=0.0; x = Z(5,1); // 状態変数は5次元 // 各ベクトル、行列の領域確保 eps= Z(x); xbase= Z(x); x_ini= Z(x); set_parameters(); // パラメータのデフォルト値を設定 // 計算する次元 comp_dimにしたがって resize x = x(1:comp_dim, 1); eps = eps(1:comp_dim, 1); xbase = xbase(1:comp_dim, 1); x_ini = x_ini(1:comp_dim, 1); // ********** 一時的なパラメータ値変更はここで。***************** // UlV1= 1.95; t_fin_bifu= 1500.0; // ****** 分岐図の可変パラメータ。...変更するときは,下の distp()も変更。*** // PRM1 #define PRM1 TdoD1 #define PRM1STR "TdoD1" #define PRM1MIN 0.7 #define PRM1MAX 12.0 #define PRM1STEP 0.05 // PRM2 #define PRM2 TV1 #define PRM2STR "TV1" #define PRM2MIN 0.8 #define PRM2MAX 5.0 #define PRM2STEP 0.01 // flag の意味 #define CHAOTIC_CHECKED 4 #define STEPOUT 3 #define CHAOTIC 2 /* see_point()しただけ */ #define NOTCHAOTIC 1 /* see_point()しただけ */ #define UNCHECKED 0 meshsize_x= Integer( (PRM1MAX-PRM1MIN)/PRM1STEP ) + 1; meshsize_y= Integer( (PRM2MAX-PRM2MIN)/PRM2STEP ) + 1; printf("# meshsize %d x %d.\n",meshsize_x, meshsize_y); flag= Z(meshsize_x, meshsize_y); statem= Z(meshsize_x, meshsize_y, x_ini); // 出発点 px= Integer( (1.1 - PRM1MIN)/PRM1STEP ) + 1; py= Integer( (1.1 - PRM2MIN)/PRM2STEP ) + 1; printf("# start point px=%d py=%d.\n",px, py); read x_ini << "xini.mat"; // 出発点でのアトラクタ上の点 printf("# read x_ini from xini.mat.\n"); statem(px, py, x_ini) = x_ini; flag(px, py)= CHAOTIC; printf("# free parameters are %s and %s.\n",PRM1STR,PRM2STR); ftrs5(t,x); rem_prm_output(); // globalパラメータの値を出す while(1){ // まだチェックしていない点を見つける。 {px, py} = find_unchecked_point(px, py); printf("\npx=%d py=%d ",px,py); if (px<0 || py<0){ break; } // 調べる点がもうなければ終了。 x_ini = statem(px, py, x_ini); // px, py でのアトラクタ上の点 // 周り８方の点のアトラクタを確かめて，flag と statem を設定 see_point(px+1, py-1, x_ini); see_point(px+1, py, x_ini); see_point(px+1, py+1, x_ini); see_point(px-1, py-1, x_ini); see_point(px-1, py, x_ini); see_point(px-1, py+1, x_ini); see_point(px, py-1, x_ini); see_point(px, py+1, x_ini); flag(px, py)= CHAOTIC_CHECKED; // チェック済の印 } output_chaotic_region_to_matfile(); // 結果を出力 printf("# CPU time= %g(%gh)\n",gettimer(),gettimer()/3600.0); } Func List find_unchecked_point(px, py) Integer px, py; { Integer j, k; // まずは手近なところを探す for(j=-1; j<=1; j++){ for(k=-1; k<=1; k++){ if (j==0 && k==0){ continue; } px= px+j; py= py+k; if ( px<1 || py<1 || px>Rows(flag) || py>Cols(flag) ){ continue; } if ( Integer(flag(px,py))==CHAOTIC ){ return {px, py}; } } } // 全部調べる for(px=1; px<=Rows(flag); px++){ for(py=1; py<=Cols(flag); py++){ if ( Integer(flag(px,py))==CHAOTIC ){ return {px, py}; } } } return {-1, -1}; // 見つからなかった場合 } Func void see_point(px, py, x_ini) Integer px, py; Matrix x_ini; { Matrix x, tmpx; Real t; Integer conv, p; Integer fd; // 範囲を外れた場合は何もしない if ( px<1 || py<1 || px>Rows(flag) || py>Cols(flag) ){ return; } // カオスであると分かっている場合は再度は調べない if ( Integer(flag(px,py)) == CHAOTIC ){ return; } if ( Integer(flag(px,py)) == CHAOTIC_CHECKED ){ return; } fd= fopen("/dev/null","w"); // 分岐図は出力しない {PRM1, PRM2} = prmvalue_for(px, py); x = x_ini; for(t=0.0; ;){ // ------------------------------------------- x = rngkut4(t,x,ftrs5,h); t = t + h; // {t, x} = rk4(ftrs5,t,x,h); // １ステップ進める。 tmpx = x; // xが間違って変更されないように tmpx に複製 // 定常状態に収束 または 既定の時間に達した場合 conv==1. {conv, p}= bifudiag_output(PRM1, secf(x), outvf(x), t, t_fin_bifu, is_identical, tmpx, fd); if ( x(Ndlt1)>1000.0 ){ p= -2; conv=1; } // 脱調 if ( conv==1 ){ break; } } // -------------------------------------------------------- fclose(fd); printf("t=%g ",t); if (p>12){ flag(px,py)= CHAOTIC; statem(px,py,x) = x; }else if (p==-2){ flag(px,py)= STEPOUT; // 脱調の時 statem(px,py,x) = x; }else{ flag(px,py)= NOTCHAOTIC; statem(px,py,x) = x; } } // 結果を出力 Func void output_chaotic_region_to_matfile() { Integer px, py; Real p1, p2; Matrix Rchaos, Rchaos2, Rnotchaos, Rstepout; Rchaos= Z(0); Rchaos2= Z(0); Rnotchaos= Z(0); Rstepout= Z(0); for(py=1; py<=Cols(flag); py++){ for(px=1; px<=Rows(flag); px++){ {p1, p2} = prmvalue_for(px, py); switch( Integer(flag(px,py)) ){ case CHAOTIC_CHECKED: Rchaos= [ Rchaos, [[p1][p2]] ]; break; case CHAOTIC: Rchaos2= [ Rchaos2, [[p1][p2]] ]; break; case NOTCHAOTIC: Rnotchaos= [ Rnotchaos, [[p1][p2]] ]; break; case STEPOUT: Rstepout= [ Rstepout, [[p1][p2]] ]; break; default: break; } } } print Rchaos >> "Rchaos.mat"; print Rchaos2 >> "Rchaos2.mat"; print Rnotchaos >> "Rnotchaos.mat"; print Rstepout >> "Rstepout.mat"; } // 格子 (x, y) でのパラメータ値を返す Func List prmvalue_for(px, py) Integer px, py; { return {PRM1MIN + (px-1)*PRM1STEP, PRM2MIN + (py-1)*PRM2STEP}; } // secf(x)=0 がポアンカレ断面となる Func Real secf(x) Matrix x; { return x(Nw1); } // 分岐図の縦軸になる値 Func Real outvf(x) Matrix x; { return x(Ndlt1)*180.0/PI; } // 同じ点とみなしていいかどうかの判定 Func Integer is_identical(x1, x2, mul) Matrix x1, x2; Real mul; // 普通は1.0。甘くしたいときに大きくする。 { Real distance, threshold; // UlV によるアトラクタの伸縮を考慮 threshold= mul * 1.0e-4*(UlV1-1.90929); distance= norm( Array(x2-x1)/Array(xbase) ); if ( distance < threshold ){ // printf("- %g %g\n",distance, threshold); return 1; }else{ return 0; } } Func void dat_output(x) Matrix x; { switch(systype){ case FirstOrderAVR: printf("%g %g %g %g\n",x(Ndlt1)/PI*180, x(Nw1), x(NeqD1), x(NgV1)); break; case WithVMdelay: default: printf("%g %g\n",Vt1, Pe1); break; } } Func void rem_dat_output_16(x, s) Matrix x; String s; { switch(systype){ case FirstOrderAVR: printf("#%.16g %.16g %.16g %.16g %s",x(1)/PI*180, x(2), x(3), x(4), s); break; case WithVMdelay: printf("#%.16g %.16g %.16g %.16g %.16g %s",x(1), x(2), x(3), x(4), x(5), s); default: break; } printf("\n"); }