Диссертация (1144294), страница 17
Текст из файла (страница 17)
fig ’] , ’compact ’) ;4849505152535455hFig2 = figure ; hold on ;plot3 ( x (: ,1) , x (: ,3) , x (: ,4) , ’ Color ’ , [0 , 0 , 1]) ;axis on ; grid on ;xlabel ( ’ x_1 ’) , ylabel ( ’ x_3 ’) , zlabel ( ’ x_4 ’) ;view (3) ; hold off ; box off ;savefig ( hFig2 ,[ ’ fitts_cont_ ’ int2str ( iStep -1) ’ _x1x3x4 . fig ’] , ’compact ’) ;end11756endБ.3Локализация хаотического решения в системе Фиттса с кусочно-гладкой и гладкой правой частью сиспользованием обратного сценария разрывной аппроксимации и метода продолжения по параметруЛистинг Б.13:fittsSatSyst.m – функция, задающая систему (1.21) скусочно-гладкой нелинейностью типа ”насыщение”.1function out = fittsSatSyst (t ,x , a0 , a1 , a2 , a3 , N )2A = [0 , 1, 0 , 0;0 , 0 , 1 , 0;0 , 0 , 0 , 1;-a0 , -a1 , -a2 , - a3 ];34567b = [0; 0; 0; 1];89c = [0; 0; -1; 0];1011sat = 1 / (2* N ) * ( abs (c ’ * x + N ) - abs (c ’ * x - N ) );1213out = A * x + b * sat ;1415endЛистинг Б.14:runFittsSat.m – функция, для построения аттрактора всистеме (1.21) с нелинейностью типа ”насыщение”.1function runFittsSat2345acc = 1e -8; RelTol = acc ; AbsTol = acc ; InitialStep = acc /10;solverOptions = odeset ( ’ RelTol ’ , RelTol , ’ AbsTol ’ , AbsTol , ...’ InitialStep ’ , InitialStep , ’ NormControl ’ ,’ on ’) ;678m1 = 0.9; m2 = 1.1; beta = 0.1;1189101112abcd====4* beta ;m1 ^2+ m2 ^2+6* beta ^2;(2*( m1 ^2+ beta ^2) ) * beta +2* beta *( m2 ^2+ beta ^2) ;( m1 ^2+ beta ^2) *( m2 ^2+ beta ^2) ;131415traj0 = [1.18318278566033; -0.595564919059928; ...-0.190135252095761; 0.453183915321858];1617tEnd = 1 e3 ; N = 5e -3;18192021% Integration of the trajectory :[t , traj ] = ode45 ( @ (t , x ) fittsSatSyst (t , x , d , c , b , a , N) , ...[0 , tEnd ] , traj0 , solverOptions ) ;222324252627figure (1) ; hold on ;plot3 ( traj (: ,1) , traj (: ,4) , traj (: ,3) , ’ Color ’ , [0.5 , 0 , 0.5])axis on ; grid on ;xlabel ( ’ x_1 ’) ; ylabel ( ’ x_4 ’) ; zlabel ( ’ x_3 ’) ;view ([38 ,10]) ; hold off ; box off ;282930313233figure (2)plot (t , [0; 0; -1; 0] ’ * traj ’ , ’ Color ’ , ’ blue ’, ’ LineWidth ’ , 1) ;xlabel ( ’t ’) ; ylabel ( ’ x_3 ’) ;axis ([100 , 200 , -0.6 , 0.6]) ;endЛистинг Б.15:fittsTanhSyst.m – функция, задающая правую часть системы(1.21) при переходе от кусочно-гладкой нелинейности типа ”насыщение” (при = 0) гладкой нелинейности типа гиперболический тангенс (при = 1).1function out = fittsTanhSatSyst (t , x , a0 , a1 , a2 , a3 , N , eps )23456A = [0 , 1, 0 , 0;0 , 0 , 1 , 0;0 , 0 , 0 , 1;-a0 , -a1 , -a2 , - a3 ];78b = [0; 0; 0; 1];910c = [0; 0; -1; 0];11911sat = 1 / (2* N ) * ( abs (c ’ * x + N ) - abs (c ’ * x - N ) );1213out = A * x + b * ( sat + eps * ( tanh (c ’ * x / N ) - sat ) ) ;1415endЛистинг Б.16:runFittsTanh.m – функция, реализующая процедурупродолжения по параметру для системы (1.21) при переходе откусочно-гладкой нелинейности типа ”насыщение” (при = 0) гладкойнелинейности типа гиперболический тангенс (при = 1)..1function runFittsTanh2345acc = 1e -8; RelTol = acc ; AbsTol = acc ; InitialStep = acc /10;solverOptions = odeset ( ’ RelTol ’ , RelTol , ’ AbsTol ’ , AbsTol , ...’ InitialStep ’ , InitialStep , ’ NormControl ’ ,’ on ’) ;67m1 = 0.9; m2 = 1.1; beta = 0.1;89101112a3a2a1a0====4* beta ;m1 ^2+ m2 ^2+6* beta ^2;(2*( m1 ^2+ beta ^2) ) * beta +2* beta *( m2 ^2+ beta ^2) ;( m1 ^2+ beta ^2) *( m2 ^2+ beta ^2) ;131415traj0 = [1.18318278566033; -0.595564919059928; ...-0.190135252095761; 0.453183915321858];1617tEnd = 1 e3 ; N = 1e -2; epsSpan = 0 : 0.5 : 1;181920for iEps = 1 : numel ( epsSpan )eps = epsSpan ( iEps ) ;21222324% Integration of the trajectory :[t , traj ] = ode45 ( @ (t , x) fittsTanhSatSyst (t , x , a0 , a1 , a2 , a3, N , eps ) , [0 , tEnd ] , traj0 , solverOptions ) ;traj0 = traj ( end ,:) ;252627figure ; hold on ;plot3 ( traj (: ,1) , traj (: ,4) , traj (: ,3) , ’ Color ’ , [0.5 , 0 , 0.5]) ;120282930axis on ; grid on ;xlabel ( ’ x_1 ’) ; ylabel ( ’ x_4 ’) ; zlabel ( ’ x_3 ’) ;view ([38 ,10]) ; hold off ; box off ;31323334353637figure ;plot (t , [0; 0; -1; 0] ’ * traj ’ , ’ Color ’ , ’ blue ’ , ’ LineWidth ’ ,1) ;xlabel ( ’t ’) ; ylabel ( ’ x_3 ’) ;axis ([100 , 200 , -0.6 , 0.8]) ;endend121Приложение ВЧисленное моделирование системы лоренцевского типаЛистинг В.1:glorenzSyst.m – функция, задающая обобщенную системуЛоренца (2.1).12function out = glorenzSyst (t , X , sigma , r , d , b )out = zeros (3 ,1) ;3% Праваяout (1) =out (2) =out (3) =45678ч а с т ь о б о б ще н н о й с и с т е м ы Л о р е н ца :- sigma * X (1) + sigma * X (2) ;r * X (1) - d * X (2) - X (1) * X (3) ;- b* X (3) + X (1) * X (2) ;endЛистинг В.2:lorenzLikeSyst.m – функция, задающая систему лоренцевскоготипа (2.2).12function out = lorenzLikeSyst (t , X , alpha , beta , lambda )out = zeros (3 ,1) ;3% Праваяout (1) =out (2) =out (3) =45678ч а с т ь с и с т ем ы ло р е н ц е в с к о го т и па :X (2) ;- lambda * X (2) - X (1) * X (3) + X (1) - X (1) ^3;- alpha * X (3) - beta *X (1) * X (2) ;endЛистинг В.3:J.m – функция, задающая матрицу линейной части системы(2.2).12212345function OUT = J (X , alpha , beta , lambda )OUT = [0 1 0;-3* X (1) ^2 - X (3) +1 - lambda -X (1) ;- beta * X (2) - beta * X (1) - alpha ];endВ.1Реализация алгоритма численного определения границ областей неустойчивости системы лоренцевскоготипаЛистинг В.4:1simulateTrajectory.m – функция.function simulateTrajectory (T , T1 , T2 , b0 , sigma0 , r0 , d0 , saveDir )23global bsigma r d ;45678b = b0 ; sigma = sigma0 ; r = r0 ; d = d0 ;function OUT = J ( X )OUT = [ - sigma sigma 0 ; r - X (3) -d -X (1) ; X (2) X (1) -b ];end910S0 = [0 , 0 , 0];111213S12_XY = sqrt ( b *( r - d ) ) ;S12_Z = r - d ;141516[V , D ] = eig ( J ( S0 ) ) ;[ ~ , IX ] = sort ( diag (D ) , ’ descend ’) ;171819S1 = [ S12_XY , S12_XY , S12_Z ];S2 = [ - S12_XY , - S12_XY , S12_Z ];202122[ V1 , D1 ] = eig ( J ( S1 )) ;[ ~ , IX1 ] = sort ( real ( diag ( D1 ) ) , ’ descend ’) ;232425delta = 0.01;x_unst_01 = S0 + delta *V (: , IX (1) ) ’;123262728x0 = 1000; y0 = x0 ; z0 = b *0.9* x0 ^2 + x0 ^2; x_unst_02 = [ x0 , y0 , z0 ];x0 = 10000; y0 = x0 ; z0 = 0; x_unst_03 = [ x0 , y0 , z0 ];293031rel_tol = 1e -8; abs_tol = 1e -8;options = odeset ( ’ RelTol ’ , rel_tol , ’ AbsTol ’ , abs_tol );3233[t , x_ust1 ] = ode45 ( @glorenzSyst , [0 , T ] , x_unst_01 , options ) ;3435hFig = figure ( ’ Visible ’ , ’ off ’) ;3637383940xStart = plot3 ( x_unst_01 (1) , x_unst_01 (2) , x_unst_01 (3) , ’.
’ , ’markersize ’ , 15 , ’ Color ’ , ’ black ’) ;hold on ;hS0 = plot3 ( S0 (1) , S0 (2) , S0 (3) , ’. ’ , ’ markersize ’ , 15 , ’ Color ’ , ’red ’) ;text ( S0 (1) , S0 (2) , S0 (3) , ’ S_0 ’ , ’ fontsize ’ , 18) ;41424344c = 1: numel ( t ) ;% # colorsh = surf ([ x_ust1 (: ,1) , x_ust1 (: ,1) ] , [ x_ust1 (: ,2) , x_ust1 (: ,2) ] , [x_ust1 (: ,3) , x_ust1 (: ,3) ] , [c (:) , c (:) ] , ’ EdgeColor ’ , ’ flat ’ , ’FaceColor ’ , ’ none ’) ;colormap ( jet ( numel ( t) ) );45464748495051525354grid on ; axis auto ;xlabel ( ’x ’) ; ylabel ( ’y ’) ; zlabel ( ’z ’) ;view ([30 ,15]) ;saveas (h , [ saveDir ’ /0 ’ ’ - [ T = ’ num2str (T , ’ %.6 g ’) ’ _b = ’ num2str (b ,’ %.2 f ’) ’ _r = ’ num2str (r , ’ %.6 g ’) ’ _d = ’ num2str (d , ’ %.6 g ’) ’ _sigma =’ num2str ( sigma , ’ %.6 g ’) ’]. png ’ ]) ;view (0 ,90) , title ( ’X - Y ’) ;saveas (h , [ saveDir ’ /0 ’ ’ - [ T = ’ num2str (T , ’ %.6 g ’) ’ _b = ’ num2str (b ,’ %.2 f ’) ’ _r = ’ num2str (r , ’ %.6 g ’) ’ _d = ’ num2str (d , ’ %.6 g ’) ’ _sigma =’ num2str ( sigma , ’ %.6 g ’) ’] _XY .
png ’ ]) ;view (0 ,0) , title ( ’X - Z ’);saveas (h , [ saveDir ’ /0 ’ ’ - [ T = ’ num2str (T , ’ %.6 g ’) ’ _b = ’ num2str (b ,’ %.2 f ’) ’ _r = ’ num2str (r , ’ %.6 g ’) ’ _d = ’ num2str (d , ’ %.6 g ’) ’ _sigma =’ num2str ( sigma , ’ %.6 g ’) ’] _XZ . png ’ ]) ;view (90 ,0) , title ( ’Y - Z ’) ;1245556saveas (h , [ saveDir ’ /0 ’ ’ - [ T = ’ num2str (T , ’ %.6 g ’) ’ _b = ’ num2str (b ,’ %.2 f ’) ’ _r = ’ num2str (r , ’ %.6 g ’) ’ _d = ’ num2str (d , ’ %.6 g ’) ’ _sigma =’ num2str ( sigma , ’ %.6 g ’) ’] _YZ . png ’ ]) ;close ( hFig ) ;5758endЛистинг В.5:123mainVarEpsB.m – функция.function mainVarEpsB (T , T1 , T2 )outDir = ’ ./ OUT / ’;sigma = 36; epsillon_min = 0.1; epsillon_max = 1.0; e_step = 0.1;4567891011% Routing the paths for plot directoriescurrProjDir = [ outDir ’/ epsilonProject / Line5 - ’ num2str (epsillon_min , ’ %.6 g ’) ’ - ’ num2str ( e_step , ’ %.6 g ’) ...’ - ’ num2str ( epsillon_max , ’ %.6 g ’) ...’/ proj ’ ’ -’ num2str (T , ’ %.6 g ’) ’ - ’ num2str ( T1 , ’ %.6 g ’) ’ - ’num2str ( T2 , ’ %.6 g ’) ];if ~ exist ( currProjDir , ’ dir ’)mkdir ( currProjDir ) ;end12for i = epsillon_min : e_step : epsillon_maxcurrDir0 = [ currProjDir ’/ epsillon = ’ num2str (i , ’ %.2 g ’) ];if ~ exist ( currDir0 , ’ dir ’)mkdir ( currDir0 ) ;endbmin = 4 * i - 0.15; step = 0.05; bmax = 4 * i - 0.05;for j = bmin : step : bmaxcurrDir = [ currDir0 ’/ b = ’ num2str (j , ’ %.2 g ’) ];if ~ exist ( currDir , ’ dir ’)mkdir ( currDir ) ;endd = i - sigma ; r = - sigma - d ;plot_and_save_single_trajectory (T , T1 , T2 , j , sigma , r ,d , currDir ) ;endend13141516171819202122232425262728end125Листинг В.6:plotBorder.m – функция отрисовки численно посчитаннойграницы области неустойчивости системы (2.1).1function plotBorder23outDir = ’ ./ OUT / ’;45678currProjDir = [ outDir ’/ epsilonProject / Divider / ’ ];if ~ exist ( currProjDir , ’ dir ’)mkdir ( currProjDir ) ;end9101112b_val = [0.5 , 0.9 , 1.3 , 1.7 , 2 , 2.4 , 2.8 , 3.3 , 3.8 , 4.2];b_val2 = [0.39 , 0.795 , 1.19 , 1.59 , 2, 2.4 , 2.8 , 3.25 , 3.7 , 4.15];epsilon_val = 0.1:0.1:1;13141516171819fig = figure ( ’ vis ’ , ’ off ’) ;plot ( b_val2 (:) , epsilon_val (:) , ’x - ’ , ’ LineWidth ’ ,2 , ’ markers ’ , 12 ,’ MarkerEdgeColor ’ , ’ red ’) ;grid on ;axis ([0 , 4.5 , 0 , 1]) ;xlabel ( ’b ’) ;ylabel ( ’\ sigma + d ’) ;202122saveas ( fig , [ currProjDir ’/ ’ ’pic - divider .