Диссертация (1144294), страница 19
Текст из файла (страница 19)
S12 * ( vU1_S1 + vU2_S1 ) / norm (vU1_S1 + vU2_S1 ) ;sol . sepaS1 = feval ( ODE . solverName , @ (t , x ) lorenzLikeSyst (t , x ,alpha , beta , lambda ) , [0 , tEnd . trans ] , initPointS1 , ODE . options ) ;if ~ isempty ( plotOptions )plotLorenzLikeSepa ( fileName , sol . sepaS1 .x , sol . sepaS1 .y ’ ,plotOptions ) ;if ODE . lim13462636465666768attrS1 = deval ( sol . sepaS1 , tLimSpan ) ;plotLorenzLikeSepa ( fileName , tLimSpan , attrS1 ’,plotLimOptions ) ;endendelsesol . sepaS1 = [];end6970warning ( ’ on ’ , ’ MATLAB : odearguments : RelTolIncrease ’) ;7172endЛистинг В.13:plotLorenzLikeSepa.m – функция для построения исохранения картинок фазового пространства с системы (2.1).1function plotLorenzLikeSepa ( fileName , t , traj , plotOptions )23S0 = [0 , 0 , 0]; S1 = [1 , 0 , 0]; S2 = [ -1 , 0 , 0];45hFigMain = figure ; hold on ;67set ( hFigMain , ’ visible ’ , plotOptions .
figVisible ) ;891011plot3 ( S0 (1) , S0 (2) , S0 (3) , ’. ’ , ’ markersize ’ , 15 , ’ Color ’ , ’ red ’) ;plot3 ( S1 (1) , S1 (2) , S1 (3) , ’. ’ , ’ markersize ’ , 15 , ’ Color ’ , ’ red ’) ;plot3 ( S2 (1) , S2 (2) , S2 (3) , ’. ’ , ’ markersize ’ , 15 , ’ Color ’ , ’ red ’) ;121314151617plot3 ( traj (: ,1) , traj (: ,2) , traj (: ,3) ) ;if plotOptions . symmetricplot3 ( - traj (: ,1) , - traj (: ,2) , traj (: ,3) ) ;endgrid on ; xlabel ( ’x ’) ; ylabel ( ’v ’) ; zlabel ( ’u ’) ;181920212223if isempty ( plotOptions .
axis )axis auto ;elseaxis ( plotOptions . axis ) ;end2425if isempty ( plotOptions . view )135view (3) ;2627elseview ( plotOptions . view ) ;2829end30313233if plotOptions . figSavesavefig ( hFigMain , [ fileName , ’. fig ’] , ’ compact ’) ;end34353637if strcmp ( plotOptions . figVisible , ’ off ’)close ( hFigMain ) ;end383940if plotOptions . colorhFigColor = figure ; hold on ;4142set ( hFigColor , ’ visible ’ , plotOptions .
figVisible ) ;43444546plot3 ( S0 (1) , S0 (2) , S0 (3) , ’. ’ , ’ markersize ’ ,15 , ’ Color ’ , ’ red ’) ;plot3 ( S1 (1) , S1 (2) , S1 (3) , ’. ’ , ’ markersize ’ ,15 , ’ Color ’ , ’ red ’) ;plot3 ( S2 (1) , S2 (2) , S2 (3) , ’. ’ , ’ markersize ’ ,15 , ’ Color ’ , ’ red ’) ;474849505152c = 1: numel ( t ) ;% # colorssurf ([ traj (: ,1) , traj (: ,1) ] , [ traj (: ,2) , traj (: ,2) ] , ...[ traj (: ,3) , traj (: ,3) ] , [ c (:) , c (:) ] , ...’ EdgeColor ’ , ’ flat ’ , ’ FaceColor ’ ,’ none ’) ;colormap ( jet ( numel ( t) ) ) ;5354grid on ; xlabel ( ’x ’) ; ylabel ( ’v ’) ; zlabel ( ’u ’) ;555657585960if isempty ( plotOptions . axis )axis auto ;elseaxis ( plotOptions . axis ) ;end6162636465if isempty ( plotOptions .
view )view (3) ;elseview ( plotOptions . view ) ;136end6667if plotOptions . figSavesavefig ( hFigColor , [ fileName , ’ _col . fig ’] , ’ compact ’) ;end68697071if strcmp ( plotOptions . figVisible , ’ off ’)close ( hFigColor ) ;end727374end7576endЛистинг В.14:figs2png.m – функция, для построения и сохраненияPNG-картинок.1function figs2png ( delta , betaList , sList )23for iBeta = 1 : length ( betaList )45beta = betaList ( iBeta ) ;678for iS = 1 : length ( sList )s = sList ( iS ) ;91011DIR = [ ’ ./ GRID_TEST / Delta = ’ num2str ( delta , ’ %.2 g ’) ’/ Beta = ’ ,num2str ( beta , ’ %.2 g ’) ];fileName = [ DIR , ’/ ’ int2str ( iS ) ’ -s = ’ , num2str (s , ’ %.16 g ’) ];12if exist ( fileName , ’ file ’)hFig = openfig ([ fileName ’. fig ’], ’ invisible ’) ;set ( hFig , ’ Position ’ , [1 , 1 , 200 , 200]) ;print ( hFig , ’ - r300 ’ , [ fileName , ’. png ’] , ’ - dpng ’) ;close ( hFig ) ;131415161718fileNameNoTrans = [ fileName ’ _noTrans ’ ];hFig2 = openfig ([ fileNameNoTrans ’.
fig ’] , ’ invisible ’) ;set ( hFig2 , ’ Position ’ , [1 , 1 , 200 , 200]) ;print ( hFig2 , ’ - r300 ’ ,[ fileNameNoTrans , ’. png ’] , ’ - dpng ’) ;close ( hFig2 ) ;1920212223end2425end1372627endendВ.2.2 Алгоритм для численного моделирования отображения ПуанкареЛистинг В.15:bifChaoticMain.m – скрипт для запуска процедурымоделирования поведения сетки точек на сечении Пуанкаре припоследовательном применении отображения Пуанкаре.1function bifChaoticMain234DIR = ’ ./ FIG / ’;delta = 0.9; beta = 2.899;56789sL = 0.7955; sR = 0.7958;numPoincareMaps = 100;isContinue = 0;lorenzLikeSystPoincareSec ( DIR , delta , sL , numPoincareMaps , isContinue ) ;1011endЛистинг В.16:lorenzLikeSystPoincareSec.m – процедура моделированияповедения сетки точек на сечении Пуанкаре при последовательномприменении отображения Пуанкаре.1function lorenzLikeSystPoincareSec ( DIR , delta , beta , s ,numPoincareMaps , isCont )23warning ( ’ off ’ , ’ MATLAB : odearguments : RelTolIncrease ’) ;45T = 1000;67alpha = delta * sqrt (1 - s ) ; lambda = s / sqrt (1 - s ) ;89101112% % Event 1 : Is trajectory in the vicinity of the eq .
S0 :function [ value , isterminal , direction ] = isInVic (~ , y ,epsHomoclin , dir )value = norm ( y ) - epsHomoclin ;isterminal = 1;138direction = dir ;1314end15161718192021% % Event 2 :function [ value , isterminal , direction ] = isOnPS0 (~ , y , p ,vNorm )value = vNorm (1) * ( y (1) - p (1) ) + vNorm (2) * ( y (2) - p (2) )+ vNorm (3) * ( y (3) - p (3) ) ;isterminal = 1;direction = 1;end2223242526272829% % Event 3 :function [ value , isterminal , direction ] = isOnPS1 (~ , y , p ,vNorm )value = [ vNorm (1) * ( y (1) - p (1) ) + vNorm (2) * ( y (2) - p (2) )+ vNorm (3) * ( y (3) - p (3) ) , ...- vNorm (1) * (y (1) + p (1) ) - vNorm (2) * ( y (2) + p (2) ) +vNorm (3) * (y (3) - p (3) ) ];isterminal = [1 , 1];direction = [1 , 1];end30313233343536% % Event 4 :function [ value , isterminal , direction ] = isTurning ( ~ , y , vS ,vSS )value = [ y (1) , y (2) , det ([ y (1) , y (2) , y (3) ; vS (1) , vS (2) , vS(3) ; vSS (1) , vSS (2) , vSS (3) ]) ];isterminal = [1 , 1 , 0];direction = [1 , -1 , 0];end373839%%odeLorenzLike = @ (t , x ) lorenzLikeSyst (t , x , alpha , beta , lambda ) ;404142S0 = [0 , 0 , 0]; S1 = [1 , 0 , 0]; S2 = [ -1 , 0 , 0];equilPoints = [ S0 ; S1 ; S2 ];434445% % Stability of the equilibria :[V , D ] = eig ( J ( S0 , alpha , beta , lambda ) ) ;13946474849505152[ ~ , IX ] = sort ( real ( diag ( D ) ) , ’ descend ’) ;fprintf ( ’ lambda_1 ( S0 ) = %s , lambda_2 ( S0 ) = %s , lambda_3 ( S0 ) = % s \ n ’,D ( IX (1) , IX (1) ) , D ( IX (2) , IX (2) ) , D ( IX (3) , IX (3) ) ) ; %if real ( D( IX (1) , IX (1) ) ) > 0fprintf ( ’ Eq .
S0 is unstable \ n ’) ;elsefprintf ( ’ Eq . S0 is stable \ n ’) ;end535455565758596061[ ~ , D1 ] = eig ( J ( S1 , alpha , beta , lambda ) ) ;[ ~ , IX1 ] = sort ( real ( diag ( D1 ) ) , ’ descend ’) ;fprintf ( ’ lambda_1 ( S1 ) = %s , lambda_2 ( S1 ) = %s , lambda_3 ( S1 ) = % s \ n ’, D1 ( IX1 (1) , IX1 (1) ) , D1 ( IX1 (2) , IX1 (2) ) , D1 ( IX1 (3) , IX1 (3) )) ; %if real ( D1 ( IX1 (1) , IX1 (1) ) ) > 0fprintf ( ’ Eq . S1 and S2 are unstable \ n ’) ;elsefprintf ( ’ Eq . S1 and S2 are stable \n ’) ;end6263646566% % Eigenvectors of J in S0vU = V (: , IX (1) ) ’ / norm ( V (: , IX (1) ) ’) ;vS = V (: , IX (2) ) ’ / norm ( V (: , IX (2) ) ’) ;vSS = V (: , IX (3) ) ’ / norm ( V (: , IX (3) ) ’) ;676869% % Parameters of ODE solver :acc = 1e -15; rel_tol = acc ; abs_tol = acc ; InitialStep = acc /10;70717273% % Choosing ’ epsHomoclin ’ for init .
point in the vicinity of the eq. s S0epsHomoclinInit = 1e -16;x_unst_01 = S0 + epsHomoclinInit * vU ;747576% % Choosing ’ epsVicinity ’ for placing the Poincare sections Pi0epsVicPS0 = 5e -1; % 3e -5;777879808182% % Define P0 - point on the Poincare section Pi_0P0 = S0 + epsVicPS0 * vS ;optionsOnPS0 = odeset ( ’ RelTol ’ , rel_tol , ’ AbsTol ’ , abs_tol , ...’ InitialStep ’ , InitialStep , ’ NormControl ’ , ’ on ’ , ...’ Events ’ , @ (t , x ) isOnPS0 (t , x , P0 , - vS ) ) ;140838485[ ~ , sepa ] = ode45 ( odeLorenzLike , [0 T ] , x_unst_01 , optionsOnPS0 ) ;sepaEnd = sepa ( end , :) ;86878889909192% % Construct a rectangular grid of points on Pi_0wNorm = [ det ([ vS (2) , vSS (2) ; vS (3) , vSS (3) ]) , ...- det ([ vS (1) , vSS (1) ; vS (3) , vSS (3) ]) , ...det ([ vS (1) , vSS (1) ; vS (2) , vSS (2) ]) ];vPS = cross ( vS , wNorm ) / norm ( cross ( vS , wNorm ) );vPSN = - cross ( vS , vPS ) / norm ( cross ( vS , vPS ) ) ;939495epsPoincareWidth = 2.5 e -2; % norm ( sepaEnd - P0 ) + epsVicPS0 /10;disp ( epsPoincareWidth ) ;96979899A = det ([ vS (2) , vS (3) ; vSS (2) , vSS (3) ]) ;B = - det ([ vS (1) , vS (3) ; vSS (1) , vSS (3) ]) ;C = det ([ vS (1) , vS (2) ; vSS (1) , vSS (2) ]) ;100101102distSepaStMan = abs ( A * sepaEnd (1) + B * sepaEnd (2) + C * sepaEnd(3) ) / ( sqrt ( A ^2 + B ^2 + C ^2) ) ;disp ( distSepaStMan ) ;103104105horNumSteps = 5;frameThikness = epsPoincareWidth / 2; % epsVicPS0 /20;106107108epsPoincareHeight = epsPoincareWidth + frameThikness ;disp ( epsPoincareHeight ) ;109110111112% % Vertex of the framefrVertex = frameVertex ( P0 , vPS , vPSN , frameThikness , ...epsPoincareWidth , epsPoincareHeight ) ;113114115116117% % Rectangular on Poincare section Pi0rectHorSize = 6e -2;rectVertSize = 6e -2;rect3D = rectangle3d ( P0 , vPS , vPSN , rectHorSize , rectVertSize ) ;118119120121% % Distance from saddle to Poincare section Pi1epsVicPS1 = 1.5 * rectVertSize ;P1 = S0 + epsVicPS1 * vU ;141122optionsOnPS1 = odeset ( ’ RelTol ’ , rel_tol , ’ AbsTol ’ , abs_tol , ’InitialStep ’ , InitialStep , ’ NormControl ’ , ’ on ’ , ’ Events ’ , @ (t , x) isOnPS1 (t , x , P1 , vU ) ) ;123124125axisSize3D = [];126127128% % Check continuation of num .
procedureif isCont == 0129130131pointsPS0 = frameGrid ( P0 , vPS , vPSN , frameThikness , ...epsPoincareWidth , epsPoincareHeight ,horNumSteps ) ;132133134[ pointsBackPS0 , pointsPS1_1 , pointsPS1_2 ] = ...execPoincareMap ( pointsPS0 , odeLorenzLike , [0 T ] ,optionsOnPS0 , optionsOnPS1 ) ;135136fileName = [ DIR ’ 00 ’ ];137138plotPoincareSec ( equilPoints , vS , vSS , rect3D , sepa , frVertex ,pointsPS0 , pointsPS1_1 , pointsPS1_2 , pointsBackPS0 , axisSize3D ,fileName ) ;139140141142143144145146147148149150lastIndex = 0;currPointsPS0 = pointsBackPS0 ;save ([ DIR ’ currPointsPS0 .